Pitfalls in Modeling Walls and Neutrals Physics in Gas Discharges Using Parallel Particle-in-Cell Monte Carlo Collision Algorithms

Owing to their ability to model the physics of low-pressure plasmas away from thermodynamical equilibrium, particle-in-cell (PIC) techniques have become one of the tools of choice to simulate the operation of many plasma devices. This trend is reinforced by the growing access to parallel computing ressources which enables tackling problems that were previously intractable with PIC techniques. However, accurate modelling of these plasmas often depends critically on the detailed description of a variety of physical phenomena ranging from microscopic to macroscopic scale and from electrons’ to neutral particles’ timescale. Among those are coupling phenomena between charged particles and neutrals. We illustrate here how the implementation of simplified models for scattering kinematics, neutrals dynamics and particle-wall interaction can affect simulation results. Until the full breadth of these effects can be captured in models, these results underline the importance of using extensive parametric scans to assess the importance of these effects.


INTRODUCTION
Low-pressure gas discharges and plasma sources are used in a great variety of applications [1,2], such as space propulsion [3,4], neutral beam injection [5] and plasma separation [6][7][8]. To optimize these devices and expand the range of application of low-pressure plasmas, numerical modeling capabilities are highly desirable. Yet, experimental results indicate that such plasmas are often far from thermal equilibrium [2]. Furthermore, the mean free path of charge particles in low-pressure discharges is often comparable to, if not greater than, the vacuum vessel dimensions. Under such conditions standard fluid models do not hold [9,10] and fluid simulation tools can not fully capture the physics at play. In these plasma regimes, kinetic modeling tools such as particle-in-cell (PIC) techniques [11][12][13] and Vlasov's codes [14,15] are required for high-fidelity simulations. Kinetic modeling however comes at a significantly larger computing cost which sets an upper limit to the plasma density and plasma volume one can model. This in turn limits the applicability of PIC techniques for whole device simulations. Yet, owing to the expansion of parallel computing, PIC techniques are increasingly used to model a large range of low-pressure plasma sources [16][17][18][19].
PIC methods solve numerically Vlasov's equation which describes a fully ionized and collisionless plasma. PIC techniques [11,13] consists in following the dynamics of a very large number of macro-particles, up to 10 13 [20], interacting with the self-generated electromagnetic (EM) fields (and also possibly superimposed externally applied fields).In the explicit PIC implementations that are used in most PIC simulations, macroparticles are pushed at each time-step by solving Lorentz equation using the EM fields at the beginning of the step. Fields are then updated by solving Maxwell's equations, or Poisson's equation in electrostatic (ES) PIC, using the new macro-particle positions and velocities. These methods are particularly suitable for parallelization and have been shown to scale very well up to millions of computing cores [20].
Yet, since plasma sources and laboratory plasmas are neither fully ionized or collisionless, they can not be modeled by standard PIC methods alone. Additional models are indeed required alongside the PIC scheme to handle neutral dynamics and also neutral interaction with charged particles in partially ionized plasmas. For dense plasmas, Coulomb collisions can also be important and require their own dedicated models. Furthermore, laboratory plasmas are inherently of finite volume. Physical boundaries are therefore facing the plasma. PIC schemes hence need to be modified to adequately model the interaction of charged and neutral particles with the walls. Implementation within PIC techniques of detailed models for each of these phenomena is key to ensure that simulation results do reproduce the physics at play in actual devices. On the other hand, modeling of these phenomena leads to additional computational cost and can have a detrimental effect on parallelization performances. Efficient implementation of these models is made even more challenging by the fact that these effects occur on vastly different timescales, ranging from electron to neutral timescale (10 −10 to 10 −4 s).
In this paper, we review how the modeling of a few specific physical phenomena can play, under some conditions, an important role on PIC simulation results, and offer insights into implementation strategies with an eye to parallelization performances. In section 2, the role played by neutrals dynamics and neutrals interaction with charged particles in low-pressure discharges is highlighted, and possible modeling options are discussed. In section 3, the physics of particle-wall interaction and its importance on the properties of low-pressure plasmas is reviewed. In section 4, the main findings are summarized. Not that this paper is not aimed at covering all aspects of particlein-cell methods, but rather at underlining the importance that the choice of certain simplifications in physical models can have on simulation results. For an in-depth and detailed description of particle-in-cell methods, the reader is referred to the many excellent reviews available [see 11,13,21,22].

Motivations for Modeling Neutral Particles
Over the years, the interplay between charged particles and neutral particles has been shown to lead to a variety of phenomena across a wide range of plasma parameters. In magnetically confined fusion plasmas, neutrals play an important role on edge physics [23] and plasma confinement [24,25]. In high-density plasma sources, neutral depletion [26] can be the cause of enhanced plasma transport [27], oscillations [28] and neutral-heating [29]. Neutrals can also control transport and separation capabilities in rotating plasmas [30]. In astrophysical plasmas, neutrals can lead to instabilities in interplanetary space [31,32], contribute to the tail structure of comets [33] and modify the properties of planetary atmospheres [34,35]. The interaction between charged particles and neutrals has also been proposed to explain the formation of the solar system [36,37] through the critical ionization velocity phenomena [38].
Particular attention must therefore be taken when handling collisions to ensure that simulation results reproduce accurately the coupling mechanisms between charged particles and neutrals. In PIC models, collisions between charged particles and neutrals are typically handled through the addition of a Monte-Carlo collision (MCC) module [39][40][41]. Between each iteration of the PIC cycle, the MCC module computes the probability for each charged particle in the simulation to undergo a collision of kind k during the time increment t.
Here v r = |v − v n | is the relative velocity, x and v are respectively the position and the velocity of the charged particle, v n and n n are respectively the neutral velocity and density, and is the cross-section for the k th collisional process, with χ the angle of rotation of the velocity after collision. For ionization processes, the angle differential cross-section is itself the result of the integration of a doubly differential cross section which contains information about the energy partitioning between the scattered particle and the secondary electron. Note that in practice, a null-collisions technique [42] is implemented to speed-up the Monte-Carlo collision step.
Equations (1,2) show that two conditions need to be satisfied to ensure an accurate description of collisions between charged particles and neutrals. First, the differential cross-section assumed in the numerical model must be adequate to ensure a physical determination of post-collision velocities. This point is illustrated in section 2.2. Second, the accurate prediction of collisions is conditioned by an accurate modeling of the distribution functions of neutral particles. This point is addressed in section 2.3. These two conditions also highlight the multi-scale nature of charged particle -neutral collision. Electrons evolve and collide with neutrals on the fastest timescales (100 ps or less for n e = 10 12 cm −3 ), while neutral dynamics is much slower (100 µs or more for T n = 300 K).

Effect of Scattering Kinematics
A useful parameter to discuss the effects of collisions between charged particles and neutrals is the ionization fraction ι = n i /(n n + n i ), where n i and n n are the ion and neutral density, respectively. For low ionization fraction, say ι ≤ 10 −3 , the effects of charged particles on neutrals can be neglected to first order. This is for example the case in radio-frequency (RF) capacitive discharges, where 10 −6 ≤ ι ≤ 10 −4 [43]. In this regime, collisions may have a strong influence on the electron energy distribution function [44], but the distribution function of neutral particles is not modified by charged particles. Eq. (1) can then be simplified by assuming a uniform and steadystate density n n of neutrals with a maxwellian distribution of temperature T n . In this regime, proper handling of collisions then hinges only on the information about scattering kinematics contained in the cross-section.

Electron -Neutral Collisions
Further simplifications can be made for electron -neutral collisions. Indeed, since T e ≫ T n and m e ≪ m n , one finds v n ≪ v e , with v e the typical velocity of an electron. The electron velocity v e can then be substituted in lieu of the relative velocity v r to recast the differential cross-section as σ k (ǫ, χ). Here ǫ = m e v e 2 /2 is the kinetic energy of an electron.
While experimental data clearly shows [see 45 and references therein] that the angular scattering of electrons rapidly becomes forward as the incident kinetic increases (ǫ ≥ 20 eV), it is common in the literature to see studies in which electron-neutral scattering kinematics is treated as isotropic. Different possible explanations can be brought forward to explain this choice. First, experimental differential cross sections σ (ǫ, χ) are not as broadly available as integrated cross sections σ (ǫ). This is because measuring σ (ǫ, χ) is much more demanding than measuring σ (ǫ). Isotropic treatment can hence be used by default. Second, compared to an isotropic treatment, implementing an anisotropic treatment in a Monte Carlo collision (MCC) module requires generating one additional random number and doing one additional table look-up (interpolation) per collision event. This essentially doubles the computing cost of the collision treatment. While collisions typically only make for a small fraction of the total computing cost of PIC-MCC simulations of gas discharges, this additional computing cost can be canceled by treating electron-neutral collisions as isotropic scattering events. A legitimate question to ask is therefore whether this simplification can affect simulation results.
For elastic collisions, results from Janssen et al. [46] suggest that using an isotropic model instead of an anisotropic model does not lead to significant differences, even at higher energies. Note though that, in order to obtain this agreement, it is essential for the momentum transfer cross-section to be held constant across the isotropic and anisotropic models. In other words, the collision frequency has to be scaled up when modeling elastic collisions as anisotropic (forward) compared to the isotropic case to preserve momentum transfer. However, isotropic modeling of elastic collision may not be accurate in nearly collisionless regimes. To illustrate this last point, consider an electron packet initially directed alongx propagating in a lowpressure neutral background for different energies ǫ b . A measure of the isotropy of this electron packet is the angle ϑ 95 which is defined as with (ϑ, φ) the velocity space coordinates (0 ≤ φ ≤ 2π the azimuth and −π/2 ≤ ϑ ≤ π/2 the inclination) and f the normalized electron packet distribution function. The time evolution of ϑ 95 obtained both for isotropic and anisotropic elastic collisions is plotted in Figure 1 for ǫ b = 5, 50 and 500 eV. Anisotropic scattering is assumed to follow the empirical angular distribution cross-section suggested by Surendra et al. [47].
Results after a few collisions (tν c ≥ 2, with ν c the momentumtransfer collision frequency) confirm minimal deviation between isotropic and anisotropic models when maintaining σ m k (ǫ) constant, in agreement with Janssen et al. [46]. On the other hand, results for tν c ∼ 0.07 and ǫ b = 500 eV show that ϑ 95 is over 50% larger when assuming isotropic collisions. This over-prediction of fast electron deflection could have implications in certain lowpressure discharges with nearly collisionless trapped electrons, such as hollow cathodes [48] and wire plasma sources [49].
For inelastic collisions, anisotropic modeling is required at higher kinetic energy (ǫ ≥ 20 eV). For instance, Surendra et al. showed that isotropic scattering leads to a large overprediction of the ionization rate in a DC glow discharge [47]. Besides DC discharges, the effects of anisotropy may generally be important in discharges featuring high-energy electrons, such as runaway electrons in atmospheric pressure discharges [45,50] and electron beams formed by sheaths in DC biased RF discharges [51].
Another feature of electron-neutral scattering often neglected in plasma modeling is how kinetic energy is partitioned between scattered and secondary electrons in electron impact ionization processes. This information is described by the doubly differential cross-section for ionization σ ion (ǫ p , ǫ s , χ) [52,53], which is related to the total ionization cross section σ ion (ǫ p ) through where ǫ iz is the ionization energy and ǫ s and ǫ p are the kinetic energy of secondary and scattered electrons, respectively. Data suggest that the post-collision kinetic energy ǫ p − ǫ iz is randomly distributed between the scattered and secondary electrons at lowenergy (≤ 50 eV), but increasingly passed preferentially onto the scattered electron as ǫ p increases [52]. Not accounting for this high-energy behavior leads to underestimating the highenergy electron population. Specifically, Tzeng and Kunhardt showed how the choice of a particular energy partitioning model affects the electron energy distribution function (EEDF) [54]. Since the energy partitioning increasingly deviates from random distribution as ǫ p increases, the accurate description of these effects is, similarly to angular scattering effects, expected to become important in plasmas featuring high-energy electrons.

Ion -Neutral Collisions
An important difference between ion-neutral and electronneutral collisions lies in the fact that the ion speed v i may be comparable to neutral speed v n . For this reason, the collision kinematics must be considered in the center-of-mass frame.  (4). Both isotropic and anisotropic models converge toward an isotropic distribution after a few collisions. However, for tν c ≤ 0.2, with ν c the momentum-transfer collision frequency, the isotropic model is observed to over-predict electron deflection in the high energy cases.
Transposing methods used for modeling neutral-neutral collisions, ion-neutral collisions are often treated as hard spheres collisions in the literature, that is to say that ion-neutral scattering is assumed isotropic [see 21]. However, experimental measurements [55] and analytic derivations [56,57] indicate that elastic scattering in symmetric systems is not isotropic. On the contrary, these findings suggest that the angular differential cross-section can be approximated by where the first term on the right-hand-side corresponds to the mostly small angle momentum-exchange scattering and the second term on the right-hand-side corresponds to the mostly large angle charge-exchange scattering [58]. Here, α > 1, β > 1, A > 0 and B > 0 are function of the incident ion kinetic energy ǫ. However, care must be taken to ensure that, no matter what the chosen form for σ (ǫ, χ) is, simulations reproduce swarm data [see [59][60][61][62]] in regimes for which the local field approximation is valid. Failing to account for anisotropic scattering thus leads to an overestimation of ions scattered at intermediate angles, and an underestimation of ions scattered at small and large angles. The influence of these effects on the properties of a xenon ion beam propagating in a rarefied gas was highlighted by Giuliano and Boyd [63]. Wang et al. further showed that the choice of a particular differential cross-section for ion-neutral elastic collisions affects the ion distribution function in a helium discharge at intermediate pressure [64]. These results suggest that particular attention should be given to these phenomena, especially for devices where ion dynamics plays a key role such as plasma thrusters [65] and ion sources [66,67]. Another motivation for capturing these effects lies in the role of fast neutral particles on secondary emission at the walls, as discussed in section 3.

Assessing the Effect of Anisotropy
For a subset of species and processes, differential cross sections have been experimentally determined and compiled in the literature [see [68][69][70][71][72]. When this information is available, the influence of anisotropy on simulation results can be directly assessed by comparing results obtained with and without anisotropic handling. In the many cases where experimental data for differential cross sections is not available, or limited to a range of energies and angles, it is sometimes possible to make use of computed cross sections. For instance, differential elastic scattering cross sections calculated using the relativistic Dirac partial-wave method for elements up to 96 amu [73]. As a last recourse, different analytical models with varying accuracy have been proposed (see [46] for a discussion of these models).
Irrespective of whether differential cross section are obtained experimentally, numerically or analytically, it is critical to ensure proper normalization of the differential cross section when comparing isotropic and anisotropic models. For example, the momentum transfer cross-section should be conserved across models for elastic electron neutral collisions. Furthermore, parametric scans are highly advisable since cross-section data can vary between studies. Ideally, one would compare the different differential cross section data sets available and run simulation cases reflecting the uncertainty of the scattering data.

Effect of the Neutral Distribution Function
Up to this point, the effects of charged particles on the neutral distribution function have been neglected on the ground that the ionization fraction is small in many plasma devices. This allowed us to consider a steady-state maxwellian distribution function for neutrals when modeling collisions.
However, certain low-pressure gas discharges feature much larger ionization fractions, with for example ι = O(1) in electron cyclotron resonance (ECR) sources [43] and helicon sources [74] operating in the mTorr range. For such conditions, coupling mechanisms with charged particles have been shown to modify both the neutral density, for example in the form of neutral depletion [26], and the neutral temperature [29]. Accurate modeling of these effects hence requires accounting for the spatial and temporal evolution of the neutral distribution function together with those of charged particles.
Even at lower ionization fractions (ι ∼ 10 −2 ), the distribution function of neutrals may be affected by charged particles under some conditions. For instance, neutral depletion has also been predicted [18] in the negative ion source prototype developed for the Neutral Beam Injector (NBI) of the International Thermonuclear Experimental Reactor (ITER) and of the Demonstration power station (DEMO) [5,75,76]. For a given background pressure, the plasma density growth with input power is in this case observed to saturate as the neutral population begins being depleted, at which point further increase in input power leads to an increase of the electron temperature T e . This variation in T e then affects the potential drop across the sheath and hence the kinetic energy of ions impacting the walls. This in turn modifies the energy distribution of neutrals formed by ion recombination at the wall. Although depletion occurs in this case for a plasma density about an hundred times smaller than the neutral density, accounting for neutral dynamics is essential to capture these effects. Furthermore, the role of neutral-wall interaction on the neutral dynamics underlines the importance of wall effects discussed in Sec. 3.
Yet another configuration where neutrals modeling may be called for are devices featuring a neutral flow. For instance, the neutral density in plasma thrusters [77,78] and ion sources [67] typically varies by two order of magnitude or more between the gas injection region and the beam propagation region. In such configurations, neutral flow properties have to be modeled a priori if charged particles are assumed to not affect neutrals, or concurrently if they are. An example of the latter is plasmas where charged particles are expected to impart momentum to neutrals [79].

Neutral Modeling
Various numerical models have been proposed and implemented to capture the dynamics of neutrals concurrently with a particlein-cell scheme for charged particles. These models can be split into two groups depending on whether or not neutrals can be described as a continuous medium. Introducing the Knudsen number Kn, defined as the ratio of the neutral mean-free-path to the characteristic length of the system, continuous media and rarefied gas correspond to Kn ≤ 1 and Kn ≥ 1, respectively. In either regime, a suite of models exists, ranging from fast but very simplified models to computationally intensive rigorous models handling the full multi-scale nature of the problem.
In regimes where neutrals can be treated as a continuous medium, i. e Kn ≤ 1, the simplest option consists in assuming a constant neutral velocity. The neutral density map is then obtained by solving a continuity equation [80]. For configurations where neutral -wall interaction is expected to play a role, Barral et al. [81] proposed to supplement the continuity equation with a momentum equation in which an adjustable coefficient is used to model the diffusion of neutrals at the wall. Although this method is arbitrary and lacks scientific background, it has proven useful to investigate certain properties of plasma thrusters [see 77,78]. In coupled models, fluid and PIC-MCC numerical schemes are concurrently advanced in time. The neutral density computed by the fluid model at a given timestep is used as input for the MCC module, while the predictions of the MCC module at a given time-step are used to derive the collision rates needed in the neutrals' continuity and momentum equations.
In gas discharges and low-temperature plasmas, the Knudsen number is often larger than 1. In this regime, Navier-Stokes equations can be inaccurate and one has to turn to kinetic methods such as Direct Simulation Monte Carlo (DSMC) method [41,82]. A significant challenge in coupling DSMC models with PIC models lies is the very large computational cost of advancing concurrently in time charged particles and neutrals. Serikov et al. [83] suggested that, by using a specific weighted Monte-Carlo scheme to couple these models, PIC and DSMC and their different timescales can be treated separately. However, this choice leads to additional constraints with respect to acceptable particle weights and time-steps. This scheme is also limited to the modeling of quasi-stationary plasmas. To address these limitations and guarantee self-consistent results, one has to simulate both charged particles and neutrals through the entire range of timescales, from electron's to neutral's [84]. In allowing for longer time-steps, implicit methods can in principle relax the constraints on the time dynamics to be simulated, but at the expense of a more complex and computationally expensive field update [85].
Even in the absence of "sub-cycling" scheme such as the one proposed by Serikov et al. [83], particular care is needed when choosing the weight of neutrals and charged species macro-particles, i. e., the number of physical particles that a macro-particle represents. Indeed, performing collisions between macro-particles requires a large number of particles per cell to obtain good statistics since pairs are randomly selected within a given cell. As discussed next, this can become an issue depending on the parallelization strategy implemented, namely particle decomposition or domain decomposition. Furthermore, the different weights for neutral and charged species macroparticles have to be accounted for in the MCC module. Different strategies, including variable weights, rejection methods and accumulation, have been suggested with their own advantages and inconvenients [see [86][87][88]].

Collision Statistics and Parallelization Strategy
Domain decomposition (DD) is generally considered the most advantageous parallelization option for PIC methods [89]. In this approach, each MPI (Message Passing Interface) thread operates on a subset of the whole domain called a subdomain. DD is often supplemented by particle decomposition where OpenMP (Open Multi-Processing) threads are attached to a given MPI thread. OpenMP threads then have access to the whole sub-domain and share an equal amount of particles. This hybrid approach is appealing since it alleviates the load balancing issues which arise in pure DD implementations when density gradients form.
On the other hand, particles of a given OpenMP process can be located anywhere in the subdomain. This means that the number of macro-particles in a given cell and in a given OpenMP thread can be limited, which degrades the collision description. Note that this issue is also found with pure particle decomposition strategies. One option to address this shortcoming and improve the resolution is to use a "sorting" algorithm which, within each subdomain, groups in memory particles that are close together in physical space [90]. Indeed, by running this "sorting" algorithm regularly, the probability to find both incident and target macro-particle in the same cell increases significantly, which is advantageous both for Monte-Carlo and DSMC algorithms [41,82]. A side benefit of particles sorting is that it limits the access to computer memory and, as a consequence, speed up the calculation of both particle trajectories and source terms. This result is illustrated in Figure 2, which shows a speed up by a factor 5 when sorting particles in 3D. Details about this numerical model can be found in Fubiani et al. [18].

IMPORTANCE OF WALLS PHYSICS
Since laboratory plasmas are inherently of finite spatial extension, some fraction of particles in the plasma necessarily interact with the walls. This includes both charged and neutral particles. In certain configurations, these wall effects have been shown to govern the plasma dynamics. For instance, secondary electron emission resulting from ion impact at the cathode is essential to the maintenance of a DC discharge [91]. Accurate simulation of laboratory plasmas hence requires appropriate numerical models to describe the interaction of particles with the wall. However, the description of these phenomena in PIC models is made particularly challenging by the range of scale-lengths and timescales involved.

Electron-Wall Interaction
The physics of electron emission from electron impact on materials is complex. Depending on its properties (energy and angle of incidence) and on the target (material and surface condition), an incident electron can be lost, elastically or inelastically backscattered, or can induce the emission of one or multiple secondary electrons [92,93]. Furthermore, the emission of an electron upon impact of an incident electron at the wall is not only characterized by a scalar probability, but also by energy and angular distributions for backscattered and emitted electrons. Yet, often in the literature, electron emission is much simplified and assumed to follow a constant probability for reemission (independent of the incident electron's energy and angle), with electrons re-emitted at a given (low) temperature [see 94,95].
The electron emission coefficient δ is defined as the number of emitted electrons per incident electron. The typical evolution of δ with incident energy ǫ observed in experiments is illustrated in Figure 3. At low energy, the emission coefficient increases with ǫ and δ = 1 for an electron energy ǫ I of a few tens of eV. The secondary emission coefficient reaches a maximum δ max for an electron energy ǫ max of a few hundreds of eV [96]. Experimental measurements indicate that ǫ I , δ max and ǫ max are function of the material and surface condition [97]. For modeling purposes, empirical analytical models for δ(ǫ) can be derived from this set of experimental measurements. In low-temperature plasma simulations, an often used model is Vaughan's model [98], which fits well experimental results at high electron energy when δ approaches δ max . This model has for example been incorporated in a modified version of the XPDC1 code [9] and used to study RF discharges [99,100] and Hall thruster [101][102][103][104]. In RF discharges, Horvath et al. have shown the importance of electron emission on sheath dynamics. In particular, they highlighted how the predicted plasma density can be strongly affected by the electron emission model for large discharge voltage [100]. In Hall thrusters, secondary electron emission from the walls can play a significant role on the electron transport perpendicular to the magnetic field [102]. Recently, another analytical model has been proposed [105] to better fit the electron emission near the first cross-over energy (ǫ ∼ ǫ I ) where multipactor effects can be important, which may be useful to simulate undesirable RF discharges on satellites. Yet, a limitation of Vaughan's model remains its inability to properly account for secondary emission from low-energy incident electrons since it assumes zero emission below a threshold incident electron energy. Scholtz and al. [93] proposed a fix to address this limitation, but at the cost of a degraded accuracy in a different energy range.
A refined model addressing these issues is that of Furman and Pivi (FP) [106]. Indeed, while it is still obtained by fitting experimental data, FP's model captures separately the three components of electron emission (true-secondaries, elastically and inelastically backscattered) as well as the energy distribution of emitted electrons. Comparison with Vaughan's model highlighted the importance of emission from low-energy electrons in the simulation of multipactor effects [107]. Another advantage of FP's model is that it captures the energy distribution of emitted electrons, which can affect the plasma parameters since it contributes to the electron energy balance in the discharge. A limit of FP's model is its large number of parameters (about 45), which limits in turn its applicability. Nevertheless, while input parameters have thus far only been determined for a few materials (copper and stainless-steel), recent experimental measurements obtained for additional materials [108] may be used to extend FP's model. In addition, parameters for metals may also be derived from a simple and validated analytical model [109].
Finally, electron emission from the wall is generally assumed to be isotropic. This assumption is justified for secondary electrons emitted from the target since the properties of the emitted electron have no relation with the incident electron. On the other hand, this justification does not hold true for elastically backscattered electrons. On the contrary, detailed modeling based on Monte Carlo simulations reveal that the angular distribution of elastically backscattered electrons is strongly dependent of the angular distribution of incident electrons. As shown in Figure 4, the angular distribution presents two lobes of emission, one in the incident direction and the other specular to the incident direction [110]. To our knowledge, this property has not yet been included in low-temperature plasma numerical models.
In summary, while existing models describing secondary electron emission by electron impact capture many important facets of this phenomenon (detailed emission process, energy distribution of secondary electrons), new features such as anisotropic emission continue to be uncovered. Consequently, it remains presumptuous to use these models for predictive simulations. Nevertheless, when compared with simpler models or use parametrically, these models can provide valuable insights into the effects of secondary emission.

Heavy Particle-Wall Interaction
Similarly to the physics of electron impact at the wall discussed in the previous paragraph, the physics of heavy particle impact at the wall is intricate and can have a significant impact on the plasma properties. Electron emission by heavy particle impact stems from two contributions, namely kinetic and potential emission. Potential emission occurs when a charged particle (ion) impact a surface with a kinetic energy ǫ greater than twice the work function of the target material W. Potential emission decreases with W and is relatively independent of ǫ [111]. On the other hand, kinetic emission depends strongly on the energy of the incident heavy particle (ion or neutral), and becomes comparable or larger than potential emission when the incident energy reaches a fraction of a keV. As a result, the secondary emission coefficient γ , defined as the number of emitted electrons per incident ion, is relatively constant for ǫ ≤ 500 eV and increases with ǫ past this energy. Here again, empirical models based on experimental data fits have been derived (see [112]), but their use is limited to heavy particle -target pairs for which detailed experimental measurements are available. Even when data is available, the experimental observation that surface conditions, such as roughness and pollution, can have a strong effect on emission [113,114] questions the validity of these models to describe experiments since surface conditions are intrinsically hard to determined. The influence of electron secondary emission by heavy particle impact can nevertheless be assessed in numerical simulations through parametric scans, as shown for example by Daksha et al. in capacitively coupled plasmas [115].
Besides the electron emission coefficient, another parameter of importance in the description of heavy particle -wall interaction is the thermal accommodation coefficient α. This coefficient relates the energy of the heavy particle leaving the wall to both the energy of the incident heavy particle and the wall temperature. This coefficient varies from α = 0 when a heavy particle is backscattered off the surface with its incident energy to α = 1 when the heavy particle looses its kinetic energy and leaves with an energy corresponding to the surface temperature. Accurate quantification of this coefficient in numerical models is important since it may, under some conditions, affect simulation results. Such dependence has for example been observed by Fubiani et al. [18] when modeling the operation of the BATMAN ion source [75]. In this source, changing α from 0.5 to 1 has been predicted to lead to a decrease of the hydrogen atom temperature by a factor of 5 and an increase of the neutral hydrogen density by a factor of 2. As seen in Figure 5, the neutral hydrogen temperature scales nearly linearly with α over the 0.5 − 1 range. Unfortunately, the limited experimental data available combined with a strong dependence of α on surface conditions makes the determination of this coefficient challenging. This cast doubts on the ability of numerical methods to reproduce accurately experimental results from first principles, especially in light of the important role of the neutrals dynamics underlined in section 2. As a baseline for a parametric scan, one can consider that the incoming particle only interacts with the atoms of the wall (i. e., a clean surface), in which case the accommodation coefficient may be estimated numerically. A database for various materials and incident particle masses has been derived using the TRIM code [116].

Wall Physics Modeling
In light of the high complexity of particle -wall interaction phenomena and the limited experimental data available to date to support models, it seems presumptuous to expect numerical simulations to predict, or even fully capture, these effects. The difficulty of accurately describing these phenomena is made even more acute in PIC codes by the various scale-lengths (atomic spacing to surface roughness of a few µm) and time-scales (electron time-scale for backscattering to neutral time-scale for accommodation) involved. Despite these shortcomings, there are important benefits in implementing wall physics models in PIC codes. First, it allows assessing the relative importance of the various processes at play. Second, parametric scans within a range of realistic experimental values can be used to refine the analysis and examine the influence of idealized experimental conditions on plasma parameters. Practically, since the impact of charged particles at the wall can lead to the emission of neutrals and reciprocally, modeling wall physics in PIC codes requires dealing with the interaction (creation, removal) of macro-particles with different weights. This situation can be handled similarly to collisions in volume between macro-particles with different weights, as discussed in section 2.4.

Wall Physics and Parallelization
In low-pressure plasma PIC modeling, the numerical integration of particle trajectories is typically responsible for most of the computing time. As a result, the implementation of particlewall models should generally have minimal effects on the overall performances. In particular, the implementation of wall models is not expected to modify the pros and cons of the different parallelization strategies discussed in section 2.4.
Yet, for configurations where wall effects dominates bulk effects, or configurations with large particle fluxes to the walls, wall effects handling may in principle lead to load unbalance between subdomains with walls and subdomains without walls. One option to alleviate this limitation is to choose a domain decomposition scheme such that the number of boundary cells is, as much as possible, equally distributed over the subdomains. However, the associated benefits might be offset by a greater number of particles crossing from one subdomain to another at each time-step. Another option for such rare cases might be to use a particle decomposition strategy.

SUMMARY
Particle-in-cell (PIC) techniques, coupled with Monte-Carlo Collision (MCC) modules, are increasingly used to simulate a variety of low-pressure plasmas. However, and besides proper implementation, the ability of PIC-MCC techniques to reproduce and simulate accurately these plasmas hinges upon the description of various phenomena. A particular challenge here lies in the fact that these phenomena cover a wide range of timescales and scale-lengths. In this paper, we discuss more specifically the importance of both neutral dynamics and wall physics.
Properties of low-pressure plasmas are often strongly affected by collisions between charged-particles and neutrals. Accounting for these effects in numerical models requires both a detailed description of the kinematics of charged particle-neutral collisions and a proper description of the neutrals distribution function f n . This is made difficult by the fact that electronneutral collisions have to be modeled on electron timescales while f n evolves on the much slower neutral particle timescale. One important feature in charged particle-neutral collision is the anisotropy of these processes, in particular for energetic particles. Failure to account for this property may, under some conditions, affect simulation results and lead to inaccurate or even unphysical results.
Many low-pressure laboratory plasmas are also strongly affected by the walls. This makes capturing in numerical models the detailed interaction of electrons, ions and neutrals with the wall a requirement. However, describing these phenomena is particularly challenging since they cover various scale-lengths (atomic spacing to surface roughness of a few µm) and timescales (electron time-scale for backscattering to neutral timescale for accommodation). Special care must therefore be taken to ensure that all potentially important processes are accounted for.
Unfortunately, while models for particles scattering and wall physics continue to be refined and the experimental database supporting these models grows, it remains ill-fated to assume that implementing these models as is in a PIC-MCC code ensures simulation results accuracy. Until global models describing reliably these complex phenomena become available, the recommended approach to build confidence in simulation results is to start from well established simplified models, and from there to include an additional effect and rely on parametric scans around realistic experimental conditions to quantify the importance of this effect.
Although the detailed scattering physics and wall physics may only be of limited importance in many low-pressure discharges, it should be emphasized that they could play a significant role in some. While, as intended in this paper, some generic guidelines can be put forward to identify conditions where these effects can be important, it most certainly does not cover all cases. This makes it even more important to assess the importance of these effects in numerical simulations.

AUTHOR CONTRIBUTIONS
All three authors RG, GF, and LG have contributed to the bibliographic survey, to the analysis and discussion of the results and to the writing of the manuscript.