Coupling plasma physics and chemistry in the PIC model of electric propulsion: Application to an air-breathing, low-power Hall thruster

This work represents a first attempt to include the complex variety of electron-molecule processes in a full kinetic particle-in-cell/test particle Monte Carlo model for the plasma and neutral gas phase in a Hall thruster. Particular emphasis has been placed on Earth’s atmosphere species for the air-breathing concept. The coupling between the plasma and the gas phase is self-consistently captured by assuming the cold gas approximation and considering gas-wall and gas recycling from the walls due to ion neutralization. The results showed that, with air molecular propellants, all the most relevant thruster performance figures degraded relative to the nominal case using Xe propellant. The main reasons can be ascribed to a reduced ionization cross-section, a larger gas ionization mean free path due to lighter mass air species, and additional electron collisional power losses. While vibrational excitations power losses are negligible, dissociation and electronic excitations compete with the ionization channel. In addition, for molecular oxygen, the large dissociation leads to even faster atoms, further reducing their transit time inside the discharge channel. Future studies are needed to investigate the role of non-equilibrium vibrational kinetics and metastable states for stepwise ionization.

so-called atmospheric-breathing electric propulsion (ABEP) [4,5] for Earth (N 2 and O 2 ) or Venus and Mars (CO 2 ). The ABEP system includes an intake to collect and compress the molecular gas [6][7][8] and an efficient electric thruster to compensate for the atmospheric drag operating on the satellite. Due to its similarity to the ramjet gas chemical propulsion used in aerospace applications, the ABEP concept is also known as RAMEP. Different electrothermal, electrostatic, and electromagnetic configurations (Hall, helicon, microwave, ECR, etc.) are being considered as possible candidates for ABEP systems [4,5].
While the experimental characterization of different thruster configurations using air species as propellant is quite robust [9,10], their numerical simulation remains in its early stages and currently lacks a detailed representation of plasma chemistry and plasma-gas-surface coupling in different global and fluid models [11][12][13]. A fully kinetic multi-dimensional simulation is a fundamental tool to study possible changes in size-geometry, electrode arrangement, and magnetic field topology to optimize the efficiency of the mass propellant utilization of the ABEP thruster.
This study addresses this question and develops a novel sophisticated Monte Carlo module in a fully kinetic particlebased model to include electron-induced processes typical of air molecular propellants such as N 2 and O 2 , and a gas mixture N 2 -O (typical composition of the atmospheric air at 200 km altitude). The model has been applied to an SPT20 Hall thruster configuration [14] and the plasma parameters and thruster performances have been compared to those for the nominal Xe propellant case.
This article is organized as follows. Section 2 describes in detail the numerical model developed in this study. The plasma and gas simulations are described, with particular attention to the Monte Carlo Collision module used to process all the electronneutral reactions for atoms and molecules of air composition. Section 3 presents and analyzes the results by comparing the plasma parameters and thruster performances for different propellants and with nominal working thruster parameters. Finally, Section 4 reports the conclusions and future development.

Model description
This section describes the model used to characterize the physics and chemistry of plasma-gas coupling in a typical lowpower Hall thruster (HT) configuration.
The model is based on a fully kinetic particle representation; namely, the particle-in-cell/Monte Carlo collision (PIC-MCC) and test particle Monte Carlo (TPMC) methods for the plasma and gas phases, respectively [15][16][17][18][19][20]. The simulated region (Figure 1) includes the internal channel discharge and the near-field plume region (up to a distance from the exit plane equal to the external diameter of the channel). The coordinate system is cylindrical and the simulation region is reduced to a two-dimensional axisymmetric plane in the radial-axial (r − z) directions. This includes both the axis of symmetry and the dielectric material of the channel walls, which is considered with its own relative permittivity, thus leading to a more rigorous representation of the plasma-dielectric interface. Particle velocities are tracked in three dimensions (v r ; v θ ; v z ). Within this representation, uniformity is assumed along the azimuthal direction: no self-consistent azimuthal electric field is generated by space charges and no azimuthal component of the external magnetic field is considered. The azimuthal position of the particles is not updated and every cell of the (plasma and gas) mesh corresponds to a ring of cross-sectional area ΔA ΔrΔz (the grid is uniform) and volume V j 2πr j ΔrΔz, where j represents the mesh index along the radial direction. Furthermore, the cells are squares with a radial cell size equal to their axial size Δr Δz. Due to the axisymmetric nature of the simulation, effects related to the off-axis location of the cathode are ignored. Figure 2 shows the typical PIC/TPMC flowchart used for the HT simulation with typical input data exchanged between the PIC and TPMC modules. Table 1 summarizes the simulation parameters used for the cases analyzed in Section 3.
The plasma and gas modules inside the code workflow ( Figure 2) are separately and sequentially solved due to the very different time scales characteristic of the plasma and gas dynamics [21]. Inside the plasma module (Section 2.1), the typical PIC tasks (free flight of the charged particles, field solver, and bulk/surface collisional processes) are iterated using a small time step (Δt PIC < 0.2/ω pe , where ω pe is the largest Langmuir frequency detected during the simulation) for 1500 time steps on a plasma grid, with a cell size smaller than the minimum Debye length detected during the simulation, Δz PIC < λ De . During all PIC iterations, the gas background is considered fixed, and the non-uniform atomic/molecular densities obtained at the end of the previous gas TPMC iterations are used in the current MCC task (Section 2.1.3) of the plasma module. At the end of the plasma module, the gas TPMC module (Section 2.2) is called: the neutral particles (atoms and molecules) are moved with their characteristic time step (Δt TPMC 10 4 Δt PIC ) for 1500 cycles on a coarser spatial mesh (Δz TPMC 10Δz PIC ). The number of cycles corresponds to the time required by the heaviest neutrals (Xe atoms) to cover the entire computational domain from the anode to the outflow boundary, thus filtering breathing mode oscillations. The collisions between neutral particles and electrons are then processed, considering the electron background to be frozen, using the electron density n e (r, z) and the electron-neutral collisional rate coefficient k eg (r, z) maps calculated during the previous plasma iteration (averaged over the last 1500 PIC cycles).
The two different PIC plasma and TPMC gas modules with their unique tasks are detailed in the following sub-sections.

PIC plasma module
The model presented in this work simulates the full plasma discharge evolution starting from the gas breakdown to plasma sustainment in a steady-state condition. Initially, the simulation domain is empty with only gas particles from the first TPMC module call (Section 2.2). A fixed prescribed number of macroelectrons is initially injected into the plume region and the electron multiplication cascade occurs, leading to gas breakdown. When a macroscopic ion current is detected at the outflow boundary (the cathode line in Figure 1), the electron current injected into the plume is automatically set as equal to the total current collected at the anode I e,inj I anode I e,anode − I i,anode .
( 1 ) This guarantees that the total current computed at the outflow boundary is globally zero, leading to a perfectly neutralized flow at all timesteps. In the steady-state condition, the anode (or discharge) current I d must be equal to the cathode current I e,inj , that partially enters the discharge channel to sustain it and partially enters the plume to neutralize the ion plume current (see Figure 3). The latter condition is not forced but verified during the simulation.  Therefore, the discharge current is not a priori prescribed in the simulation but is instead computed self-consistently, resulting in a function of the imposed working conditions: discharge voltage, magnetic field, mass flow rate, and chemical composition of the gas propellant.
Due to the axisymmetric nature of the simulation, the electron emission from the out-of-axis cathode location cannot be reproduced and the electrons are injected uniformly in position (by considering the cylindrical metrics) in the plume region, with a full Maxwellian distribution in the velocity space, which is characterized by a temperature T e0 = 2 eV. Electron macro-particles, just like ions, are eliminated from their respective particle lists, whenever they cross a simulation boundary, including an external boundary of the simulation domain or a material boundary. The current version of the code does not perform selective reflection (based on the electron kinetic energy) at the downstream boundaries, as done in other reports [22,23]. Consequently, a numerical plasma sheath develops close to the downstream boundary to guarantee a perfect local ambipolar current condition.
At each PIC time step, the electron charge is deposited on the four cell nodes through an area weighting interpolation function. When the electric field is solved on the mesh (Section 2.1.1), it is interpolated back to the particle location by using the same interpolation function. The Boris algorithm is used to integrate the Lorentz force and the particles are moved according to the leapfrog method [15]. Due to the large ion mass (according to the propellant used in this study), a subcycling method is used to move ions and deposit their charge to the nodes only every 15 electron PIC cycles (Δt ion 15Δt PIC ). The electric field used to push the ions is calculated as the average over the last electron cycles. The related speed-up can reach up to 30%.
Anomalous electron transport contributions due to azimuthal fluctuations are considered as additional scattering with a prescribed frequency where the fitting coefficient k is set to 0.2 by matching the results of experimental measurements obtained using Xe as a propellant [24]. The anomalous collision output velocity is computed as for an elastic isotropic one, although a recent study demonstrated how the isotropic and anisotropic characters of anomalous scattering can lead to different discharge characteristics [25]. The transport due to electron-wall interaction (so-called nearwall conductivity) is self-consistently considered due to the radial nature of the simulation. Finally, the model implemented no numerical tricks (reduced ion mass, increased vacuum permittivity [26], or geometrical scaling [27]) to accelerate the execution time of the simulation.

Poisson equation solver
The Poisson equation is solved in both the plasma and the dielectric material regions (green regions in Figure 1) to consider the effect of the relative permittivity ϵ ε/ε 0 and more realistic boundary conditions instead of imposing the normal electric field at the plasma-dielectric interface. The 2D(r, z) axial-symmetric generalized Poisson equation in strong (differential) form is: If integrated in a volume Ω j,k around the node (j, k) (with j and k representing the node index along the radial and axial directions, respectively), this leads to the weak (variational) form of the Poisson equation: where C j,k is the enclosing contour surface around the node (j, k) and dn ndS is the differential unit normal vector ( Figure 4). The total charge Q tot j,k corresponds to the sum of all charges present inside the volume Ω j,k , including both volumetric free charges Q V j,k and static surface charges Q S j,k deposited on the plasma-dielectric interface, and is updated when an electron/ ion hits the dielectric surface. Expanding the gradient operator, breaking the integral in the four different contributions corresponding to the four different faces of the volume Ω j,k in 2D, approximating the derivatives with central differences schemes, and evaluating the permittivity on the contour surface element as an average over the two cells (the permittivity is defined at the cell centers; i.e., in a staggered mesh), the discretized version of Poisson Eq. 3 becomes where (assuming Δr Δz) the different coefficients are: On the axis of symmetry (r 0), the volume Ω 0k reduces to a cylindrical volume, and, by imposing zϕ/zr 0, the Poisson coefficients reduce to At the remaining boundaries (the cathode and anode lines), simple Dirichlet conditions are applied. The linear algebraic system derived from Eq. 4 (and including the boundary conditions) is finally solved using the Petsc package [28].

Plasma-wall interaction
When an electron strikes the dielectric wall, the number of secondary electrons emitted is selected from the Vaughan formula [29].
using the following fitting parameters corresponding to BNSiO 2 : Y max 2.016, E max 299, k 0.563. For simplicity, all secondary emitted electrons are considered true secondaries: their energy is sampled from a half-Maxwellian flux distribution with temperature T see 2 eV, corresponding to a cosine distribution in the emission polar angle and a random distribution in the azimuthal angle. All ions hitting the dielectric surfaces are deleted from the particle list since they are assumed to be neutralized and reinjected from the dielectric surfaces as neutrals. A counter is activated to update the flux of the different neutral species g (henceforth with g or G we will denote the generic neutral species used in the present study, N, N 2 , O, O 2 , or Xe) at the inner iw and outer ow walls as a function of the axial location Γ g,iw (z) and Γ g,ow (z). These quantities are used in the next neutral module call to consider ion recycling from the walls. In the cases of nitrogen and oxygen atomic ion N + and O + , the recombination coefficients (associative recombination and emission as neutral molecules) γ N + 0.07 and γ O + 0.17 [30,31] are used, respectively.

Monte Carlo collision plasma model
For simplicity, the present model simulates only electron-gas collisions. Previous HT studies showed the negligible effect of Coulomb processes involving electron and ion species [18], while ion-neutral collisions have a relevant impact on the plume expansion and plume-spacecraft interactions. For electronatom collisions, elastic scattering, electronic excitation, and single ionization are considered, while vibrational excitations, dissociation, and dissociative ionization are additionally considered for electron-molecule collisions. Collisions between electron and molecular oxygen O 2 also provide the possibility to include the production of negative ions O − through the dissociative attachment process. The production of multiply charged ions of atoms and molecules is neglected due to the low power range investigated.
An important assumption, valid for low-pressure plasmas, is that the neutral particle target G is always found by the electron projectile in the electronic ground states since the spontaneous relaxation from an electronically excited state is considered to be much faster than the electron collisional time [32]. This is true for radiative excited states, but not for Frontiers in Physics frontiersin.org 05 metastable states. The possible implications of two-step ionization processes (ionization of metastable states) will be investigated in the future by adding metastable states to the ground state in the gas module. In addition, the vibrational kinetics of molecules are deactivated and all the molecules are considered to be in their vibrational ground level ] = 0.
The different cross-sections used in the model are reported and discussed in Section 2.1.3.1, while the Monte Carlo algorithm used for the electron-neutral collisions is presented in Section 2.1.3.2.

Electron-gas cross-section set for air species and Xe
This section presents the different electron-induced processes involving the air species (N, N 2 , O, and O 2 ) and Xe included in the model. The corresponding cross-sections are also discussed.  Detailed overviews of the relevant cross-sections with nitrogen species are given in Brunger and Buckman [33], Itikawa [34], Tabata et al. [35,36], Thorsteinsson and Gudmundsson [30], and Kawaguchi et al. [37]. The complete reaction set used in the present model is reported in Table 2.
The cross-section for elastic scattering by the molecular ground state N 2 (X 1 Σ g + ) is taken from [34][35][36]. The vibrational excitation cross-section corresponds to the sum of the transitions from ] = 0 to all possible final states [38]. In our model, we have considered the electron impact excitation of the following 14 electronic excited states of the nitrogen molecule [34][35][36]: Furthermore, we assumed that the total molecular dissociation cross-section includes three channels: 1) a resonant dissociation through the intermediation of the unstable nitrogen negative ion N 2 + e → (N 2 − )* → 2N + e leading to a sharp narrow peak around 10 eV; 2) predissociation; i.e., excitation to certain electronic and vibrational states that automatically dissociate into two nitrogen atoms by an internal conversion from an excited state towards a repulsive dissociative state [40]; and 3) the remaining contributions for E > 14 eV. The cross-sections of molecular and dissociative ionization with energy thresholds of 15.58 eV and 24.34 eV, respectively, are taken from [34][35][36].
Dissociative double ionization with the production of two atomic ions is also considered, with an energy threshold of 34 eV and cross-section taken from [34][35][36].
The cross-section for elastic scattering by the atomic nitrogen ground state N( 4 S) is taken from [39]. In addition to atomic excitations to the N( 2 D) state at 2.39 eV and to the N( 2 P) state at 3.57 eV, 24 other electronic excitations with thresholds >10 eV are considered [39]. The cross-section for atomic ionization (at 14.54 eV) is taken from [39].
The cross-sections for the collisions between electron and the oxygen species (O and O 2 ) were presented by Brunger and Buckman [33], Itikawa [41], Vahedi and Surendra [42], Gudmundsson et al. [31], and Vass et al. [43]. The complete reaction set is reported in Table 3.
The cross-section used for the elastic scattering of electrons by the ground state of the oxygen molecule O 2 (X 3 Σ g − ) is that recommended by Itikawa [41]. The vibrational cross-section of O 2 shows completely different behavior in the energy regions above and below 1 eV. At energies >1 eV, the cross-section shows a broad peak at about 10 eV. At 0.2-1 eV, the cross-section consists of a set of very sharp peaks due to a temporary electron capture of O to form a negative ion state O − ( 2 Π g ). In between, the cross-section value is very small. Here, we have included the contributions of four different transitions ] = 0 → 1-4 and neglected the resonant peaks due to their small width and energy range compared to the typical electron energies in HT.   Frontiers in Physics frontiersin.org 07 threshold energies (≈4 eV) are very close. Finally, the electron energy loss spectrum for O 2 electronic excitation shows a broad peak ranging from 7 to 9.5 eV. This is called the Schumann-Runge (SR) continuum and is caused by the excitation of the O 2 (B 3 Σ u − ) state. Its energy threshold is 6.12 eV. The excitation of the SR continuum also contributes to a neutral dissociation   [41] and the corresponding ionization potential are 12.06 eV and 18.73 eV, respectively. For the dissociative attachment from the ground state oxygen molecule, the cross-section is taken from Itikawa [41]. The threshold energy is 4.2 eV. The incident electron loses its energy, which is absorbed by the oxygen molecule to form O 2 − , which subsequently dissociates to form the fragments O and O − . The potential energy for the O+ O − pair is 3.63 eV above the ground state potential for O 2 . The remaining incident electron energy (E inc -3.63) eV is divided between the two fragments.
The cross-section for elastic collisions of electrons with oxygen atoms is taken from Itikawa and Ichimura [45]. The cross-sections for electron-impact excitation of the atomic oxygen ground 3 P state to the 1 D, 1 S, and 3 P 0 excited states and to the Rydberg 5 S 0 , 3 S 0 , 5 P, 3 P, 5 D 0 , and 3 D 0 states are taken from Laher-Gilmore [46]. The cross-section for electronimpact ionization is taken from [47]. The ionization potential is 13.62 V.
Finally, the electron-Xe reaction set (Table 4) includes the elastic momentum transfer of the Xe ground state 5p 6

MCC algorithm for electron-gas collisions
The collision probability of the ith electron can be written as [16].
where n g is the density of the gas species g in the electron home cell (calculated during the previous TPMC module) and σ tot,g is the total cross-section considering all possible collisional processes l between an electron and the neutral species g. In the case of molecular oxygen O 2 , for example, the total crosssection is σ tot,O2 l σ l,O2 σ el,O2 + σ vib,O2 + σ exc,O2 + σ diss,O2 + σ ion,O2 The sum of the vibrational and electronic excitations and for the dissociation in Eq. 11 includes all possible channels (Table 3) for these reactions.
The null Monte Carlo collision method [42] is used to save computational time and the Nanbu method [48] simultaneously determines: 1) whether an electron collides (henceforth, R 01 is a generic random number uniformly distributed in the range [0,1]); 2) which gas species g is selected as the gas target partner among the total N g neutral species The Monte Carlo module is called every 60 PIC cycles (Δt coll 60Δt PIC ) to guarantee that the collisional operator is correctly resolved; i.e., P i < 0.1 [16, 42,48].
Under the cold gas approximation used (the gas-particle velocity is much smaller than the electron velocity), the relative and electron velocities are equal. Neglecting the electron mass compared to the neutral mass m ≪ M, the kinematic collision equations [42,48] are greatly simplified. No kinematic information related to the neutral target particles is required and the electron post-collisional velocity is generally calculated as follows in all elastic and inelastic processes: whereṽ v(1 − E th E ) 1/2 and h [ṽ ⊥ cos η, −ṽ rṽθ cos η+ṽṽz sin η v⊥ , −ṽ rṽz cos η+ṽṽ θ sin η v⊥ ]. Here, v and E are the electron pre-collisional velocity and energy, respectively, E th is the threshold energy of the inelastic process (third column of Tables 2-4), while χ and η are the scattering and azimuthal angles, respectively. All collisions are assumed to be isotropic; i.e., the collisional angles are sampled as [16, 42,48]: In the case of an ionization event (single ionization for an atom and molecule or dissociative ionization for a molecule) the Frontiers in Physics frontiersin.org 09 total post-collisional energy is equally distributed between the scattered (parent) and ejected (progeny) electrons: where E ion is the ionization potential (Tables 2-4). The scattered electron is deflected at an angle χ assumed to be isotropic [Eq. 16a], while the ejected electron is deflected at an angle χ′, assumed to be orthogonal to the χ direction (χ′ χ + π/2). The velocity components of the ion byproduct are sampled from a full Maxwellian distribution with the local temperature T g (r, z) computed in the previous gas module call. In the case of dissociative ionization or attachment, the ionic fragment randomly shares the post-collisional energy and is assumed to be isotropically scattered.
In the last MCC call of the PIC module iteration, the reaction rate constant of the different collisional events l of the electron-g collision k l,g (r, z) 〈σ l,g (E)v〉 (18) is calculated [here the brackets refer to the average over all electrons belonging to the home cell centered around the location (r, z)] and passed together with the electron density n e (r, z) to the following TPMC module. These quantities are fixed and used for all TPMC iterations in the next gas module call to calculate the collision probability of the neutral particles with the electron target, as explained in Section 2.2.

Gas-phase test particle Monte Carlo model
During the neutral module iteration, gas particles are allowed to enter the simulation domain from the anode location and from the inner and outer wall (with an axial-dependent flux Γ g,iw (z) and Γ g,ow (z) from the ion recycling calculated in the previous plasma module call) with a half-Maxwellian flux distribution with temperature T g0 500 K. Neutral-wall interactions are simulated by assuming a reflection coefficient R n 1 for all species with an energy accommodation coefficient α 1 (neutral emitted from the surface with wall temperature T gw 500 K). The same recombination coefficients γ N γ N + 0.07 and γ O γ O + 0.17 [30,31] used for ion-wall neutralization are applied for nitrogen and oxygen atoms, respectively. For example, given the recombination coefficient of nitrogen atoms γ N 0.07, a nitrogen atom is always reflected when it hits the wall in some form: as an atom N with probability p 1-0.07 93% and as a molecule N 2 with probability p = 0.07/ 2 = 3.5% (with the factor of two because two N atoms are needed to produce an N 2 molecule).
Once injected, neutral macro-particles move in straight lines until they experience a collision (see Section 2.2.1), hit the walls, or exit the simulation domain.

Collisional algorithm for the TPMC module
For every ith gas particle of the g species, the neutral-electron collisions probability is calculated as P i n e (r, z)Δt TPMC k tot,g (r, z) n e (r, z)Δt TPMC Nc l 1 k l,g (r, z), where the electron density n e (r, z) and the total reaction rate coefficient k tot,g (r, z) Nc l 1 k l,g (r, z) in the neutral particle home cell come from the previous PIC iteration, Eq. 18. Here, the cold gas approximation makes the reaction rate coefficient dependent only on the local electron energy distribution function. The acceptance-rejection method, Eq. 12, is used to decide if the neutral collides, and Eq. 14 to decide the individual collisional events c. Regarding the outcome, should a collision actually occur, the following applies: • elastic and excitation events: the neutral particle velocity and internal energy level (vibrational and electronic) remain unchanged; • dissociation (dissociative ionization) event: two atoms (one atom) are (is) created and the molecule is deleted from the particle list; the total energy available (depending on the actual dissociation channel) is randomly distributed between the two by-products and they are assumed to be isotropically scattered. • ionization event: the neutral particle is simply removed from the neutral particle list.

FIGURE 6
Magnetic field topology (streamlines in black and B field magnitude in Tesla) used for all analyzed cases.
Frontiers in Physics frontiersin.org  Frontiers in Physics frontiersin.org 11 Due to the low-pressure regime, the possible gas-gas processes are not considered.

Results
The HT configuration selected in the present study corresponds to a low-power SPT-20 HT type [14]. This decision is suggested by the need to minimize the size of the full satellite system to reduce the drag to be compensated and the propellant mass flow rate required to operate. The channel length is z ch 1 cm, while the inner and outer radii are r in 0.5 cm and r out 1 cm, respectively. The power level is P d 80 W, corresponding to a discharge voltage V d 200 V and current I d~0 .4 A. For the different propellants analyzed, the mass flow rate is fixed at _ m 0.1 mg/s. The imposed magnetic field has a classical dipolar topology, as shown in Figure 6: a convex magnetic lens with a dominant radial component, and a bellshaped behavior along the axial direction with a maximum B max 180 Gauss at the exit plane along the channel centerline. The engineering parameters used in the simulations are shown in Table 5. Figure 7 shows the two-dimensional plots of the most relevant HT plasma and gas parameters for the Xe propellant case (Case A): 1) electric potential ϕ(V), 2) electron temperature T e (eV), 3) electron density n e (m −3 ), 4) atomic density n Xe (m −3 ), 5)-6) electronic excitation and ionization rate coefficients, k eE,Xe and k ion,Xe (m 3 s −1 ) (Eq. 18), 7) ionization source term S ion,Xe k ion,Xe n Xe n e (m −3 s −1 ), and 8) ion flux Γ Xe + + (m −3 s −1 ) obtained by averaging over the last 1500 PIC cycles for plasma quantities (a-c, e-h) and over the last 1500 TPMC iterations for gas quantity (d).

Case A: Xenon propellant test case
The maps demonstrate the main features of the typical HT discharge and plume. The results are strongly correlated with the magnetic field topology used: most of the electric field is concentrated near the exhaust, where the neutral density is so low (n Xe < 5·10 19 m −3 ) that the electron mobility is completely determined by the anomalous collision frequency prescribed by Eq. 2. Outside the channel, the potential drop is about 160 V, which is two-thirds of the total discharge potential. There, the electron temperature reaches its maximum of 30 eV. The plasma density peaks (n e~1 .5·10 18 m −3 ) in a region shifted inward relative to the acceleration region: this peak corresponds to the exit plane while the axial electric field peaks 5 mm downstream from the exit plane. The ionization rate S ion,Xe ( Figure 7G) highlights the role of gas recycling at the exit plane by ion neutralization at the inner and outer walls. The average value (along the inner and outer walls) of the secondary electron emission coefficient results in < Y > 0.79, while the thrust is slightly larger than 1 mN. All these quantities agree fairly well with experimental results [14].

Comparisons of different air species propellants
This section compared cases using N 2 (case B), O 2 (case C) and N 2 -O mixture (case D, with 50%-50% in terms of mass flow rate and corresponding to the air composition at an altitude of approximately 200 km) as propellants to the nominal case using Xe (case A), focusing on plasma and performance parameters, which are summarized in Table 6. Case D does not consider air composition changes due to the gas dynamics inside the intake, which can lead to the possible formation of N, O 2 , and NO species entering the thruster channel discharge from the anode [49].   Figure 7 for cases B, C, and D, respectively. The dissociation coefficient for cases B and C is defined as the ratio between the atomic and total (atomic G and molecular parent G 2 ) gas density averaged over the entire computational domain while the mass utilization efficiency η m is defined as the ratio between the ion mass flow rate exhausted into the plume and the neutral propellant mass flow rate injected from the anode with m i and Γ i the mass and particle flow of the i-th ionic species, leaving the computational domain across the cathode line ( Figure 1). The sum extends over all the ion plume species (last column of Table 6).
The most evident result is that the thrust values for all air species propellants (cases B, C, and D) are lower than the corresponding values relative to the Xe propellant (case A). The main reason is ascribed to the lower mass utilization efficiency η m , which is also confirmed by the lower electron density peak (in all the cases always located at the exit plane). Through an in-depth analysis, one can attribute the low mass utilization efficiency of molecular air species to three concomitant effects, which are analyzed below.
1) Lower ionization reaction rate coefficient k ion,g (m 3 s −1 ) [Eq. 18]. Xe is characterized by a smaller ionization potential; at its peak (corresponding to an electron energy E of~100 eV), the Xe ionization cross-section is about three times larger (σ ion,Xe max 6 × 10 −20 m 2 ) than the corresponding cross- Frontiers in Physics frontiersin.org section for atomic and molecular air species, as shown in the full red curves in Figure 5. Figure 11 corroborates this hypothesis by reporting the axial profile along the thruster channel centerline of the ionization reaction rate coefficient of the dominant gas species for the different cases (molecular nitrogen N 2 for case B and atomic oxygen O for cases C and D). This behavior reflects the axial profile of the electron temperature; the highest value after the Xe propellant case A (red curve) corresponds to the ionization of molecular nitrogen (Case B, blue curve). After Xe, N 2 shows the largest ionization cross-section peak σ ion,N2 max 2 × 10 −20 m 2 ( Figure 5B).
2) Nitrogen N 2 and oxygen O 2 molecules are lighter (approximately five-fold) than Xe atoms, which directly affect the thrust (for a given total thrust efficiency and discharge power, ions that are exhausted faster produce a lower thrust) but also indirectly affect the ionization efficiency. The faster the neutral propellant particles, the smaller their transit time through the channel, reducing Frontiers in Physics frontiersin.org 15 their chance of being ionized by the electron cloud. The corresponding ionization mean free path λ ion,g v g n e k ion,g is larger, thus reducing the propellant mass utilization, for a given thruster chamber length. The ion mean free path is further increased by a smaller electron density n e ( Figure 12A) and ionization rate k ion,g ( Figure 11). Additionally, for case C (molecular oxygen O 2 propellant), the high dissociation rate (α 0.9) leads to the production of even lighter atoms (their initial energy after molecular dissociation is of the order of eV, as discussed in 2.1.3.1) thus further worsening the mass utilization efficiency and the thrust compared to the hypothetical case of lower O 2 dissociation degree.
To better understand this point, Figure 12 shows the axial profiles along the channel thruster centerline of the 1) electron density n e and 2) the generalized propellant ionization mean free path for the different cases

FIGURE 11
Axial profiles along the channel thruster centerline of the ionization reaction rate coefficient k ion,G (m 3 s −1 ) (Eq. 18) for the dominant gas species for the different cases (molecular nitrogen N 2 for case B and atomic oxygen O for cases C and D). The black vertical line corresponds to the exit plane location.

FIGURE 12
Axial profiles along the channel thruster centerline with (A) electron density n e (m −3 ) and (B) propellant ionization mean free path λ ion (cm) (Eq. 23). The black vertical line corresponds to the exit plane location. with χ g representing the fraction of the g-th neutral species. By following the Melikov-Morozov criterion [50] (in which a mass utilization efficiency >0.8 requires a minimum channel length z ch twice as large as the ionization length λ ion ), this study shows that the use of air propellant would require longer channel lengths relative to the Xe propellant case A by factors of 4 for the case of molecular oxygen O 2 propellant (case C) and 10 for the cases of molecular nitrogen N 2 and air mixture N 2 -O propellant (cases B and D, respectively).
3) Additional reaction channels lead to extra electron power losses (alternative to ionization). Table 7 and Figure 13 report the repartition of the total power spent by electrons in the main gas propellant collisional channels for the different cases. Although the power losses in molecular vibrational MeV excitations are negligible (P MeV~1 0 −2 W) due to the relatively high electron temperature range in HT (T e > 5 eV), the electronic excitations (atomic AeE and molecular MeE) and molecular dissociations (Mdiss) compete with the propellant ionization. As a figure of merit, one can take the ratio between the power dissipated in ionization (atomic and molecular) and that dissipated in dissociations and electronic excitations. The best performance is obtained with Xe with a value of PAion PAeE = 1.26 with >55% of the total collisional power budget channeled towards the propellant ionization; for the molecular propellants and the air mixture, this value is close or even <1. In the case of molecular nitrogen N 2 propellant (case B), considerable power is dissipated and equally distributed between molecular electronic excitation MeE and dissociation Mdiss. In the molecular oxygen case C, a significant power is dissipated in molecular dissociation Mdis, which is compensated by the large atomic ionization Aion. Unfortunately, the energy level structure and cross-sections of the atomic oxygen also lead to a favorable electronic excitation AeE. O + ions are almost totally produced by the stepwise ionization of atoms produced by dissociation. In contrast, for the molecular nitrogen propellant (case B), N + ions are mostly produced by direct molecular dissociative ionization. The role of negative ions O − for the molecular oxygen case C is negligible since they comprise <1% of total positive ions due to the relatively high electron temperature in HT with respect to usual gas discharges. Finally, the air mixture in case D is characterized by a larger contribution of the atomic oxygen ionization over the molecular nitrogen ionization: the power dissipated is 1.5 times larger, despite the smaller channel transit time. This occurs due to the smaller O atom ionization potential with respect to the corresponding N 2 molecular ionization (almost 2 eV smaller).

Conclusion
This study has developed a PIC-TPMC model to simulate the coupling between plasma and gas in HTs when using molecular propellants that are relevant for air-breathing applications. The model features two different alternating modules for plasma species (PIC) and neutral gas species (TPMC), which are iterated and coupled until convergence. The PIC module features standard algorithms, except for the Poisson solver, which computes the electric potential solution within the dielectric material through a generalized Poisson equation formulation. This allows the use of more realistic boundary conditions. Collisional events include electron-neutral collisions and have been modeled through the well-known null-collision method (against a background of frozen non-uniform neutrals) within the PIC module and using pre-computed (from the PIC) collisional rates for the TPMC module (against a background of frozen non-uniform electrons). A large variety of collisional processes are included; among these, elastic scattering, electronic excitation, and single ionization are considered for collisions with atoms, while vibrational excitations, dissociation, and dissociative ionization are additionally simulated for collisions with molecules (by adding the dissociative electron attachment for O 2 ).
The model has demonstrated how the degradation of the thruster performance compared to the nominal case using Xe can be ascribed mainly to the lower elementary mass of the molecular propellants. This produces a lower propellant mass utilization (and, hence, thrust efficiencies) for a fixed thruster chamber length. The results show that the latter must be elongated by factors of 4 and 10 for O 2 and N 2 /air propellants, respectively, to provide more efficient ionization of molecular and atomic by-products.
The Xe propellant case shows a total power dissipated towards ionization larger than 55%. No other air-species propellants show such a large value due to additional electron energy collisional loss channels. The main losses are due to dissociation (only partially compensated by the following atomic ionization) for O 2 and molecular electronic excitation and dissociation (equally distributed) for N 2 .

FIGURE 13
Repartition of the electron power losses (normalized to the total collisional power loss) among the main electron-gas collisional channels (MeE: molecular electronic excitation, Mdis: molecular dissociation, Mion: molecular ionization, AeE: atomic electronic excitation, and Aion: atomic ionization) for the different cases.
Frontiers in Physics frontiersin.org Finally, this work has assumed that the molecules are in their vibrational ground states and has not considered all the possible metastable states of the atoms and molecules. Additional studies are needed to investigate the role of non-equilibrium vibrational kinetics for molecules and metastable electronic states for atoms and molecules, which may be relevant for the selected propellants and increase the overall efficiency of the molecular propellant cases, through stepwise ionization [51]. Further research is also needed to address unexplored physical phenomena such as the role of ioninduced secondary electron emission or ion/neutral collisions, which have been neglected in the present work.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors without undue reservation.