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.

1 Introduction

The different concepts of electric propulsion [1] often use heavy rare gas species (Ar, Kr, Xe) as propellants due to their atomic nature and relatively high mass-to-ionization energy ratios. Studies in recent years have investigated the use of alternative molecular propellants for economic and environmental reasons (I₂, H₂O) [2, 3] and the potential for harvesting in-situ residual gas for long-lived low-orbit missions. The latter has led to the so-called atmospheric-breathing electric propulsion (ABEP) [4, 5] for Earth (N₂ and O₂) or Venus and Mars (CO₂). The ABEP system includes an intake to collect and compress the molecular gas [6–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–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 particle-based model to include electron-induced processes typical of air molecular propellants such as N₂ and O₂, and a gas mixture N₂-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 electron-neutral 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.

2 Model description

This section describes the model used to characterize the physics and chemistry of plasma-gas coupling in a typical low-power 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–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_{\mathrm{r}};v_{\theta };v_{\mathrm{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 $\Delta A=\Delta r\Delta z$ (the grid is uniform) and volume $V_j=2\pi r_j\Delta r\Delta 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 $\Delta r=\Delta z$. Due to the axisymmetric nature of the simulation, effects related to the off-axis location of the cathode are ignored.

Figure 1

Figure 1. Computational domain with boundary conditions.

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.

Figure 2

Figure 2. Flowchart of the present Hall thruster PIC-TPMC model.

Table 1

Table 1. Simulation parameters used for the present PIC-DSMC model.

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 ($\Delta t_{\mathrm{P}\mathrm{I}\mathrm{C}}<0.2/{\omega }_{\mathrm{p}\mathrm{e}},$ where ${\omega }_{\mathrm{p}\mathrm{e}}$ 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, $\Delta z_{\mathrm{P}\mathrm{I}\mathrm{C}}<{\lambda }_{\mathrm{D}\mathrm{e}}$. 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 ($\Delta t_{\mathrm{T}\mathrm{P}\mathrm{M}\mathrm{C}}={10}^4\Delta t_{\mathrm{P}\mathrm{I}\mathrm{C}}$) for 1500 cycles on a coarser spatial mesh ($\Delta z_{\mathrm{T}\mathrm{P}\mathrm{M}\mathrm{C}}=10\Delta z_{\mathrm{P}\mathrm{I}\mathrm{C}}$). 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_{\mathrm{e}}\left(r,z\right)$ and the electron-neutral collisional rate coefficient $k_{\mathrm{e}g}\left(r,z\right)$ 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.

2.1 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 macro-electrons 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_{\mathrm{e},\mathrm{i}\mathrm{n}\mathrm{j}}=I_{\mathrm{a}\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}}=I_{\mathrm{e},\mathrm{a}\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}}-I_{\mathrm{i},\mathrm{a}\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}}.(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.

Figure 3

Figure 3. Scheme representing the electron injection scheme used in the model.

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 ($\Delta t_{\mathrm{i}\mathrm{o}\mathrm{n}}=15\Delta t_{\mathrm{P}\mathrm{I}\mathrm{C}})$. 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

${\nu }_{\mathrm{a}\mathrm{n}\mathrm{o}}=k\frac{{\omega }_{\mathrm{c}\mathrm{e}}}{16},(2)$

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 near-wall 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.

2.1.1 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 $ϵ=\epsilon /{\epsilon }_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:

$\nabla \bullet \left[ϵ\left(r,z\right)\nabla \varphi \left(r,z\right)\right]=-\frac{\rho \left(r,z\right)}{{\epsilon }_0}.(3)$

If integrated in a volume ${\Omega }_{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:

${\oint }_{C_{j,k}}\left[ϵ\left(r,z\right)\nabla \varphi \left(r,z\right)\right]\bullet \mathrm{d}\mathit{n}=-Q_{j,k}^{\mathrm{t}\mathrm{o}\mathrm{t}}/{\epsilon }_0(4)$

where $C_{j,k}$ is the enclosing contour surface around the node ($j,k$) and $\mathrm{d}\mathit{n}=\mathit{n}\mathrm{d}S$ is the differential unit normal vector (Figure 4). The total charge $Q_{j,k}^{\mathrm{t}\mathrm{o}\mathrm{t}}$ corresponds to the sum of all charges present inside the volume ${\Omega }_{j,k}$, including both volumetric free charges $Q_{j,k}^V$and static surface charges $Q_{j,k}^S$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 ${\Omega }_{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

$a_{j,k}{\varphi }_{j-1,k}+b_{j,k}{\varphi }_{j+1,k}+c_{j,k}{\varphi }_{j,k-1}+d_{j,k}{\varphi }_{j,k+1}+e_{j,k}{\varphi }_{j,k}+Q_{j,k}^{\mathrm{t}\mathrm{o}\mathrm{t}}/{\epsilon }_0=0,(5)$

where (assuming $\Delta r=\Delta z$) the different coefficients are:

$a_{j,k}=\pi \left(r_j-\frac{\Delta r}{2}\right)\left[{ϵ}_{j-\frac{1}{2},k-\frac{1}{2}}+{ϵ}_{j-\frac{1}{2},k+\frac{1}{2}}\right](\mathrm{6a})$

$b_{j,k}=\pi \left(r_j+\frac{\Delta r}{2}\right)\left[{ϵ}_{j+\frac{1}{2},k-\frac{1}{2}}+{ϵ}_{j+\frac{1}{2},k+\frac{1}{2}}\right](\mathrm{6b})$

$c_{j,k}=\pi r_j\left[{ϵ}_{j+\frac{1}{2},k-\frac{1}{2}}+{ϵ}_{j-\frac{1}{2},k-\frac{1}{2}}\right](\mathrm{6c})$

$d_{j,k}=\pi r_j\left[{ϵ}_{j+\frac{1}{2},k+\frac{1}{2}}+{ϵ}_{j-\frac{1}{2},k+\frac{1}{2}}\right](\mathrm{6d})$

$e_{j,k}=-a_{j,k}-b_{j,k}-c_{j,k}-d_{j,k}.(6e)$

Figure 4

Figure 4. Scheme representing the Poisson equation discretization used in the model.

On the axis of symmetry ($r=0$), the volume ${\Omega }_{0k}$ reduces to a cylindrical volume, and, by imposing $\partial \varphi /\partial r=0$, the Poisson coefficients reduce to

$a_{0,k}=0(\mathrm{7a})$

$b_{0,k}=\pi \frac{\Delta r}{2}\left[{ϵ}_{\frac{1}{2},k-\frac{1}{2}}+{ϵ}_{\frac{1}{2},k+\frac{1}{2}}\right](\mathrm{7b})$

$c_{0,k}=\pi \frac{\Delta r}{4}{ϵ}_{0,k-\frac{1}{2}}(\mathrm{7c})$

$d_{0,k}=\pi \frac{\Delta r}{4}{ϵ}_{0,k+\frac{1}{2}}(\mathrm{7d})$

$e_{0,k}=-b_{0,k}-c_{0,k}-d_{0,k}.(7e)$

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].

2.1.2 Plasma-wall interaction

When an electron strikes the dielectric wall, the number of secondary electrons emitted is selected from the Vaughan formula [29].

$Y=Y_{\max }{\left(\frac{E}{E_{\max }}{e}^{1-\frac{E}{E_{\max }}}\right)}^k,(8)$

using the following fitting parameters corresponding to BNSiO₂: $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_{\mathrm{s}\mathrm{e}\mathrm{e}}=$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 re-injected 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₂, O, O₂, or Xe) at the inner iw and outer ow walls as a function of the axial location ${\Gamma }_{g,iw}\left(z\right)$ and ${\Gamma }_{g,ow}\left(z\right)$. 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) ${\gamma }_{N^{+}}=$0.07 and ${\gamma }_{O^{+}}=$ 0.17 [30, 31] are used, respectively.

2.1.3 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 electron-atom 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₂ 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

$G^{*}\to G+h\nu (9)$

is considered to be much faster than the electron collisional time [32]. This is true for radiative excited states, but not for 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.

2.1.3.1 Electron-gas cross-section set for air species and Xe

This section presents the different electron-induced processes involving the air species (N, N₂, O, and O₂) 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.

Table 2

Table 2. Reactions between electron and the neutral nitrogen species (N₂ and N) considered in the proposed PIC-TPMC model.

The cross-section for elastic scattering by the molecular ground state N₂(X¹Σ_g⁺) is taken from [34–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–36]: A³Σ_u⁺, B³Π_g, W³Δ_u, B′³Σ_u⁻, a′¹Σ_u⁻, a¹Π_g, w¹Δ_u, b¹Π_u, b'¹Σ_u⁺, c¹Π_u, c'¹Σ_u⁺, C³Π_u, E³Σ_g^+, and a′′¹Σ_g⁺. 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₂ + e → (N₂⁻)^* → 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–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–36].

The cross-section for elastic scattering by the atomic nitrogen ground state N(⁴S) is taken from [39]. In addition to atomic excitations to the N(²D) state at 2.39 eV and to the N(²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₂) 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.

Table 3

Table 3. Reactions between electron and the neutral oxygen species (O₂ and O) considered in the proposed PIC-TPMC model.

The cross-section used for the elastic scattering of electrons by the ground state of the oxygen molecule O₂(X³Σ_g⁻) is that recommended by Itikawa [41]. The vibrational cross-section of O₂ 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⁻(²Π_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. The cross-sections for the excitations [44] to the two metastable O₂(a¹Δ_g) (energy threshold 0.98 eV) and O₂(b¹Σ_g⁺) (energy threshold 1.63 eV) are distinguished while the excitations to the three states O₂(A³Σ_u⁺), O₂(A′³Δ_u), and O₂(c¹Σ_u⁻) are represented as a unique cross-section since the corresponding threshold energies (≈4 eV) are very close. Finally, the electron energy loss spectrum for O₂ 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₂(B³Σ_u⁻) state. Its energy threshold is 6.12 eV. The excitation of the SR continuum also contributes to a neutral dissociation O₂(B³Σ_u⁻) → O(³P) + O(¹D). The cross-sections for the electron-impact dissociation of the oxygen molecule are taken from [44]. The excitation to the 6.12 eV level leads to dissociation into O(³P) + O(³P), excitation to the 8.4 eV level leads to dissociation into O(³P) + O(¹D), and excitation to the 9.97 eV level leads to dissociation into O(¹D) + O(¹D). The corresponding released energy magnitudes to each pair of heavy fragments are 1.03 eV, 1.27 eV, and 0.88 eV, respectively. The cross-sections for electron-impact and dissociative ionization of the oxygen molecule ground state are taken from Itikawa [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₂⁻, 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₂. 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 ³P state to the ¹D, ¹S, and ³P⁰ excited states and to the Rydberg ⁵S⁰, ³S⁰, ⁵P, ³P, ⁵D⁰, and ³D⁰ states are taken from Laher-Gilmore [46]. The cross-section for electron-impact 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 1}S₀ showing the typical Ramsauer minimum at 0.62 eV and 74 electronic excitations from 8.32 eV to 12.58 eV taken from [39]. Finally, the single ionization of the Xe ground state with an energy threshold of 12.13 eV is from Hayashi [44].

Table 4

Table 4. Reactions between electrons and Xe in the proposed PIC-TPMC model.

All cross-sections used in the model are reported in Figures 5A–E.

Figure 5

Figure 5. Electron-neutral collision cross-sections used in the model for (A) atomic nitrogen N, (B) molecular nitrogen N₂, (C) atomic oxygen, (D) molecular oxygen O₂, and (E) xenon.

2.1.3.2 MCC algorithm for electron-gas collisions

The collision probability of the ith electron can be written as [16].

$P_i\cong {\sum }_g{\sigma }_{\mathrm{t}\mathrm{o}\mathrm{t},g}\left(E_i\right)v_i{n}_g\left(\mathrm{r},\mathrm{z}\right){\Delta t}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}},(10)$

where $n_g$is the density of the gas species $g$ in the electron home cell (calculated during the previous TPMC module) and ${\sigma }_{\mathrm{t}\mathrm{o}\mathrm{t},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₂, for example, the total cross-section is

${\sigma }_{\mathrm{t}\mathrm{o}\mathrm{t},{\mathrm{O}}_2}={\sum }_l{\sigma }_{l,O_2}={\sigma }_{\mathrm{e}\mathrm{l},{\mathrm{O}}_2}+\sum {\sigma }_{\mathrm{v}\mathrm{i}\mathrm{b},{\mathrm{O}}_2}+\sum {\sigma }_{\mathrm{e}\mathrm{x}\mathrm{c},{\mathrm{O}}_2}+\sum {\sigma }_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{s},{\mathrm{O}}_2}+{\sigma }_{\mathrm{i}\mathrm{o}\mathrm{n},{\mathrm{O}}_2}+{\sigma }_{\mathrm{d}\mathrm{i}\mathrm{o}\mathrm{n},{\mathrm{O}}_2}+{\sigma }_{\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{t},{\mathrm{O}}_2}.(11)$

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

$P_i\ge R_{01}(12)$

(henceforth, R₀₁ 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

$\frac{{\sum }_{l=1}^{g-1}{\sigma }_{\mathrm{t}\mathrm{o}\mathrm{t},l}{n}_l}{{\sum }_{l=1}^{N_g}{\sigma }_{\mathrm{t}\mathrm{o}\mathrm{t},l}{n}_l}<R_{01}\le \frac{{\sum }_{l=1}^g{\sigma }_{\mathrm{t}\mathrm{o}\mathrm{t},l}{n}_l}{{\sum }_{l=1}^{N_g}{\sigma }_{\mathrm{t}\mathrm{o}\mathrm{t},l}{n}_l}(13)$

3) which particular collisional event c occurs among the total $N_c$electron-neutral $\mathrm{e}-g$ processes (elastic, electronic excitation, and ionization if $g$ is an atom and elastic, vibrational and electronic excitation, dissociation, ionization and dissociative ionization if $g$ is a molecule)

$\frac{{\sum }_{l=1}^{c-1}{\sigma }_{l,g}}{{\sum }_{l=1}^{N_c}{\sigma }_{l,g}}<R_{01}\le \frac{{\sum }_{l=1}^c{\sigma }_{l,g}}{{\sum }_{l=1}^{N_c}{\sigma }_{l,g}}.(14)$

The Monte Carlo module is called every 60 PIC cycles (${\Delta t}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}={60\Delta t}_{\mathrm{P}\mathrm{I}\mathrm{C}}$) 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\ll 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:

${\mathit{v}}^{\prime }=\stackrel{\sim }{\mathit{v}}\mathit{cos}\chi -\mathit{h}\mathit{sin}\chi (15)$

where $\stackrel{\sim }{\mathit{v}}=\mathit{v}{\left(1-\frac{E_{\mathrm{t}\mathrm{h}}}{E}\right)}^{1/2}$ and $\mathit{h}=\left[{\stackrel{\sim }{v}}_{\perp }\mathit{cos}\eta ,-\frac{{\stackrel{\sim }{v}}_r{\stackrel{\sim }{v}}_{\theta }\mathit{cos}\eta +\stackrel{\sim }{v}{\stackrel{\sim }{v}}_z\mathit{sin}\eta }{{\stackrel{\sim }{v}}_{\perp }},-\frac{{\stackrel{\sim }{v}}_r{\stackrel{\sim }{v}}_z\mathit{cos}\eta +\stackrel{\sim }{v}{\stackrel{\sim }{v}}_{\theta }\mathit{sin}\eta }{{\stackrel{\sim }{v}}_{\perp }}\right]$. Here, $\mathit{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 $\chi$ and $\eta$ 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]:

$\cos \chi =1-2R_{01}(\mathrm{16a})$

$\eta =2\pi R_{01}.(\mathrm{16b})$

In the case of an ionization event (single ionization for an atom and molecule or dissociative ionization for a molecule) the total post-collisional energy is equally distributed between the scattered (parent) and ejected (progeny) electrons:

$E_{\mathrm{s}\mathrm{c}}=\frac{R_{01}}{2}\left(E-E_{\mathrm{i}\mathrm{o}\mathrm{n}}\right)(\mathrm{17a})$

$E_{\mathrm{e}\mathrm{j}}=\left(1-\frac{R_{01}}{2}\right)\left(E-E_{\mathrm{i}\mathrm{o}\mathrm{n}}\right).(\mathrm{17b})$

where $E_{\mathrm{i}\mathrm{o}\mathrm{n}}$ is the ionization potential (Tables 2–4). The scattered electron is deflected at an angle $\chi$ assumed to be isotropic [Eq. 16a], while the ejected electron is deflected at an angle ${\chi }^{\prime }$, assumed to be orthogonal to the $\chi$ direction (${\chi }^{\prime }=\chi +\pi /2$).

The velocity components of the ion byproduct are sampled from a full Maxwellian distribution with the local temperature $T_g\left(r,z\right)$ 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}\left(r,z\right)=\langle {\sigma }_{l,g}\left(E\right)v\rangle (18)$

is calculated [here the brackets refer to the average over all electrons belonging to the home cell centered around the location $\left(r,z\right)$] and passed together with the electron density $n_{\mathrm{e}}\left(r,z\right)$ 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.

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 ${\Gamma }_{g,iw}\left(z\right)$ and ${\Gamma }_{g,ow}\left(z\right)$ 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 $\alpha =$1 (neutral emitted from the surface with wall temperature $T_{gw}=$500 K). The same recombination coefficients ${\gamma }_N={\gamma }_{N^{+}}=0.07$ and ${\gamma }_O={\gamma }_{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 ${\gamma }_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₂ with probability p = 0.07/2 = 3.5% (with the factor of two because two N atoms are needed to produce an N₂ 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.

2.2.1 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\cong n_{\mathrm{e}}\left(r,z\right){\Delta t}_{\mathrm{T}\mathrm{P}\mathrm{M}\mathrm{C}}{k}_{\mathrm{t}\mathrm{o}\mathrm{t},g}\left(r,z\right)=n_{\mathrm{e}}\left(r,z\right){\Delta t}_{\mathrm{T}\mathrm{P}\mathrm{M}\mathrm{C}}{\sum }_{l=1}^{N_c}{k}_{l,g}\left(r,z\right),(19)$

where the electron density $n_{\mathrm{e}}\left(r,z\right)$ and the total reaction rate coefficient $k_{\mathrm{t}\mathrm{o}\mathrm{t},g}\left(r,z\right)={\sum }_{l=1}^{N_c}{k}_{l,g}\left(r,z\right)$ 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.

Due to the low-pressure regime, the possible gas-gas processes are not considered.

3 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_{\mathrm{c}\mathrm{h}}=$ 1 cm, while the inner and outer radii are $r_{\mathrm{i}\mathrm{n}}=$ 0.5 cm and $r_{\mathrm{o}\mathrm{u}\mathrm{t}}=$ 1 cm, respectively. The power level is $P_{\mathrm{d}}=$80 W, corresponding to a discharge voltage $V_{\mathrm{d}}=$ 200 V and current $I_{\mathrm{d}}$∼ 0.4 A. For the different propellants analyzed, the mass flow rate is fixed at $\dot{m}=$ 1.0 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 bell-shaped 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 6

Figure 6. Magnetic field topology (streamlines in black and B field magnitude in Tesla) used for all analyzed cases.

Table 5

Table 5. Engineering parameters used for the SPT20 simulations.

3.1 Case A: Xenon propellant test case

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 $\varphi$(V), 2) electron temperature $T_{\mathrm{e}}$(eV), 3) electron density $n_{\mathrm{e}}$ (m⁻³), 4) atomic density $n_{\mathrm{X}\mathrm{e}}$ (m⁻³), 5)–6) electronic excitation and ionization rate coefficients, $k_{\mathrm{e}\mathrm{E},\mathrm{X}\mathrm{e}}$ and $k_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{X}\mathrm{e}}$ (m³s⁻¹) (Eq. 18), 7) ionization source term $S_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{X}\mathrm{e}}=k_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{X}\mathrm{e}}{n}_{\mathrm{X}\mathrm{e}}{n}_{\mathrm{e}}$ (m⁻³s⁻¹), and 8) ion flux ${\Gamma }_{{\mathrm{X}\mathrm{e}}^{+}}+$ (m⁻³s⁻¹) 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).

Figure 7

Figure 7. Two-dimensional maps of the most relevant HT plasma and gas parameters for the SPT20 nominal case using Xe (case A) as the propellant: (A) electric potential $\varphi$ (V), (B) electron temperature $T_{\mathrm{e}}$ (eV), (C) electron density $n_{\mathrm{e}}$ (m⁻³), (D) atomic density $n_{\mathrm{X}\mathrm{e}}$ (m⁻³), (E) and (F) electronic excitation and ionization rate coefficients, $k_{\mathrm{e}\mathrm{E},\mathrm{X}\mathrm{e}}$ and $k_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{X}\mathrm{e}}$ (m³s⁻¹) (Eq. 18), (G) ionization source term $S_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{X}\mathrm{e}}$ (m⁻³s⁻¹) and (H) ion flux ${\Gamma }_{X{e}^{+}}$ (m⁻²s⁻¹) 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).

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_{\mathrm{X}\mathrm{e}}<$ 5∙10¹⁹ m⁻³) 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_{\mathrm{e}}\sim$ 1.5∙10¹⁸ m⁻³) 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_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{X}\mathrm{e}}$ (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].

3.2 Comparisons of different air species propellants

This section compared cases using N₂ (case B), O₂ (case C) and N₂-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₂, and NO species entering the thruster channel discharge from the anode [49].

Table 6

Table 6. Global performance and main plasma parameters corresponding to the different propellants.

Figures 8–10 report the same plasma and gas parameters of 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₂) gas density averaged over the entire computational domain

$\alpha =\frac{n_G}{n_G+n_{G_2}},(20)$

while the mass utilization efficiency ${\eta }_{\mathrm{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

${\eta }_{\mathrm{m}}=\frac{{\displaystyle \sum }_i{m_i\Gamma }_i}{\dot{m},},(21)$

with $m_i$ and ${\Gamma }_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).

Figure 8

Figure 8. Two-dimensional maps of the most relevant HT plasma and gas parameters for the SPT20 nominal case using N₂ (case B) as the propellant: (A) electric potential $\varphi$ (V), (B) electron temperature $T_{\mathrm{e}}$ (eV), (C) electron density $n_{\mathrm{e}}$ (m⁻³), (D) atomic density $n_{{\mathrm{N}}_2}$ (m⁻³), (E) and (F) electronic excitation and ionization rate coefficients, $k_{\mathrm{e}\mathrm{E},{\mathrm{N}}_2}$ and $k_{\mathrm{i}\mathrm{o}\mathrm{n},{\mathrm{N}}_2}$ (m³s⁻¹) (Eq. 18), (G) ionization source term $S_{\mathrm{i}\mathrm{o}\mathrm{n},{\mathrm{N}}_2}$ (m⁻³s⁻¹) and (H) ion flux ${\Gamma }_{N_2^{+}}$ (m⁻²s⁻¹) 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).

Figure 9

Figure 9. Two-dimensional maps of the most relevant HT plasma and gas parameters for the SPT20 nominal case using O₂ (case C) as the propellant: (A) electric potential $\varphi$ (V), (B) electron temperature $T_{\mathrm{e}}$ (eV), (C) electron density $n_{\mathrm{e}}$ (m⁻³), (D) atomic density $n_{\mathrm{O}}$ (m⁻³), (E) and (F) electronic excitation and ionization rate coefficients, $k_{\mathrm{e}\mathrm{E},\mathrm{O}}$ and $k_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{O}}$ (m³s⁻¹) (Eq. 18), (G) ionization source term $S_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{O}}$ (m⁻³s⁻¹) and (H) ion flux ${\Gamma }_{O^{+}}$ (m⁻²s⁻¹) 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).

Figure 10

Figure 10. Two-dimensional maps of the most relevant HT plasma and gas parameters for the SPT20 nominal case using N₂ and O (case D) as the propellant: (A) electric potential $\varphi$ (V), (B) electron temperature $T_{\mathrm{e}}$ (eV), (C) electron density $n_{\mathrm{e}}$ (m⁻³), (D) atomic density $n_{\mathrm{O}}$ (m⁻³), (E) and (F) electronic excitation and ionization rate coefficients, $k_{\mathrm{e}\mathrm{E},\mathrm{O}}$ and $k_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{O}}$ (m³s⁻¹) (Eq. 18), (G) ionization source term $S_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{O}}$ (m⁻³s⁻¹) and (H) ion flux ${\Gamma }_{O^{+}}$ (m⁻²s⁻¹) 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).

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 ${\eta }_{\mathrm{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_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{g}}$ (m³s⁻¹) [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 (${\sigma }_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{X}\mathrm{e}}^{\max }\cong$ 6 × 10⁻²⁰ m²) than the corresponding cross-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₂ 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₂ shows the largest ionization cross-section peak ${\sigma }_{\mathrm{i}\mathrm{o}\mathrm{n},{\mathrm{N}}_2}^{\max }\cong$ 2 × 10⁻²⁰ m² (Figure 5B).

Figure 11

Figure 11. Axial profiles along the channel thruster centerline of the ionization reaction rate coefficient $k_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{G}}$ (m³s⁻¹) (Eq. 18) for the dominant gas species for the different cases (molecular nitrogen N₂ for case B and atomic oxygen O for cases C and D). The black vertical line corresponds to the exit plane location.

2) Nitrogen N₂ and oxygen O₂ 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 their chance of being ionized by the electron cloud. The corresponding ionization mean free path

${\lambda }_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{g}}=\frac{v_g}{n_{\mathrm{e}}{k}_{\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{g}}}(22)$

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_{\mathrm{e}}$ (Figure 12A) and ionization rate $k_{\mathrm{i}\mathrm{o}\mathrm{n},g}$ (Figure 11). Additionally, for case C (molecular oxygen O₂ propellant), the high dissociation rate ($\alpha =$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₂ dissociation degree.

Figure 12

Figure 12. Axial profiles along the channel thruster centerline with (A) electron density $n_{\mathrm{e}}$ (m⁻³) and (B) propellant ionization mean free path ${\lambda }_{\mathrm{i}\mathrm{o}\mathrm{n}}$ (cm) (Eq. 23). The black vertical line corresponds to the exit plane location.

To better understand this point, Figure 12 shows the axial profiles along the channel thruster centerline of the 1) electron density $n_{\mathrm{e}}$ and 2) the generalized propellant ionization mean free path for the different cases

${\lambda }_{\mathrm{i}\mathrm{o}\mathrm{n}}={\sum }_g{{\chi }_g\lambda }_{\mathrm{i}\mathrm{o}\mathrm{n},g}(23)$

with ${\chi }_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_{\mathrm{c}\mathrm{h}}$ twice as large as the ionization length ${\lambda }_{\mathrm{i}\mathrm{o}\mathrm{n}}$), 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₂ propellant (case C) and 10 for the cases of molecular nitrogen N₂ and air mixture N₂-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_{\mathrm{M}\mathrm{e}\mathrm{V}}\sim$ 10⁻² W) due to the relatively high electron temperature range in HT ($T_{\mathrm{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 $\frac{P_{\mathrm{A}\mathrm{i}\mathrm{o}\mathrm{n}}}{P_{\mathrm{A}\mathrm{e}\mathrm{E}}}$ = 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₂ 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₂ molecular ionization (almost 2 eV smaller).

Table 7

Table 7. Electron power loss (W) repartition in the main collisional molecular (MeE-electronic excitation, Mdis-dissociation, and Mion-ionization) and atomic (AeE-electronic excitation and Aion-ionization) channels for the different propellants.

Figure 13

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.

4 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₂).

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₂ and N₂/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₂ and molecular electronic excitation and dissociation (equally distributed) for N₂.

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 ion-induced 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.

Author contributions

FT developed the code and analyzed the results of the model. PM implemented the algorithm to solve the Poisson equation. All authors contributed to the writing of this article.

Funding

This work was funded by the Italian Ministry of University and Research (MUR) (project PON “CLOSE to the Earth”, No. ARS ARS01–00141).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. O’Reilly D, Herdrich G, Kavanagh DF. Electric propulsion methods for small satellites: A review. Aerospace (2021) 8(22). doi:10.3390/aerospace8010022

CrossRef Full Text | Google Scholar

2. Szabo J, Pote B, Paintal S, Robin M, Hillier A, Branam RD, et al. Performance evaluation of an iodine-vapor Hall thruster. J Propuls Power (2012) 28:848–57. doi:10.2514/1.B34291

CrossRef Full Text | Google Scholar

3. Schwertheim A, Knoll A. Low power thrust measurements of the water electrolysis Hall effect thruster. CEAS Space J (2021) 14:3–17. doi:10.1007/s12567-021-00350-y

CrossRef Full Text | Google Scholar

4. Schönherr T, Komurasaki K, Romano F, Massuti-Ballester B, Herdrich G. Analysis of atmosphere-breathing electric propulsion. IEEE Trans Plasma Sci (2015) 43(1):287–94. doi:10.1109/TPS.2014.2364053

CrossRef Full Text | Google Scholar

5. Zheng P, Wu J, Zhang Y, Biqi W. A comprehensive review of atmosphere-breathing electric propulsion systems. Int J Aerospace Eng (2020) 2020:1–21. doi:10.1155/2020/8811847

CrossRef Full Text | Google Scholar

6. Romano F, Espinosa-Orozco J, Pfeiffer M, Herdrich G, Crisp N, Roberts P, et al. Intake design for an atmosphere-breathing electric propulsion system (ABEP). Acta Astronaut (2021) 187:225–35. doi:10.1016/j.actaastro.2021.06.033

CrossRef Full Text | Google Scholar

7. Ferrato E, Giannetti V, Tisaev M, Lucca Fabris A, Califano F, Andreussi T. Rarefied flow simulation of conical intake and plasma thruster for very low earth orbit spaceflight. Front Phys (2022) 2022:823098. doi:10.3389/fphy.2022.823098

CrossRef Full Text | Google Scholar

8. Singh LA, Walker MLR. A review of research in low Earth orbit propellant collection. Prog Aerospace Sci (2015) 75:15–25. doi:10.1016/j.paerosci.2015.03.001

CrossRef Full Text | Google Scholar

9. Gurciullo A, Fabris AL, Cappelli MA. Ion plume investigation of a Hall effect thruster operating with Xe/N2 and Xe/air mixtures. J Phys D: Appl Phys (2019) 52:464003. doi:10.1088/1361-6463/ab36c5

CrossRef Full Text | Google Scholar

10. Marchioni F, Cappelli MA. Extended channel Hall thruster for air-breathing electric propulsion. J Appl Phys (2021) 130:053306. doi:10.1063/5.0048283

CrossRef Full Text | Google Scholar

11. Taploo A, Lin L, Keidar M. Analysis of ionization in air-breathing plasma thruster. Phys Plasmas (2021) 28:093505. doi:10.1063/5.0059896

CrossRef Full Text | Google Scholar

12. Garrigues L. Computational study of Hall-effect thruster with ambient atmospheric gas as propellant. J Propuls Power (2012) 28(2):344–54. doi:10.2514/1.B34307

CrossRef Full Text | Google Scholar

13. Pekker L, Keidar M. Analysis of airbreathing Hall-effect thrusters. J Propuls Power (2012) 28(6):1399–405. doi:10.2514/1.B34441

CrossRef Full Text | Google Scholar

14. Loyan AV, Maksymenko TA. Performance investigation of SPT-20M low power Hall effect Thruster. In: IEPC-2007-100 paper, 30th International Electric Propulsion Conference; September 17-20, 2007; Florence, Italy (2007).

Google Scholar

15. Tskhakaya D, Matyash K, Schneider R, Taccogna F. The particle-in-cell method. Contrib Plasma Phys (2007) 47(8-9):563–94. doi:10.1002/ctpp.200710072

CrossRef Full Text | Google Scholar

16. Taccogna F. Monte Carlo Collision method for low temperature plasma simulation. J Plasma Phys (2015) 81(1):305810102. doi:10.1017/S0022377814000567

CrossRef Full Text | Google Scholar

17. Taccogna F, Longo S, Capitelli M, Schneider R. Stationary plasma thruster simulation. Comput Phys Commun (2004) 164(1-3):160–70. doi:10.1016/j.cpc.2004.06.025

CrossRef Full Text | Google Scholar

18. Taccogna F, Longo S, Capitelli M, Schneider R. Surface-driven asymmetry and instability in the acceleration region of Hall thruster. Contrib Plasma Phys (2008) 48(4):375–86. doi:10.1002/ctpp.200810061

CrossRef Full Text | Google Scholar

19. Taccogna F, Minelli P. Three-dimensional particle-in-cell model of Hall thruster: The discharge channel. Phys Plasmas (2018) 25(6):061208. doi:10.1063/1.5023482

CrossRef Full Text | Google Scholar

20. Taccogna F, Minelli P, Asadi Z, Bogopolsky G. Numerical studies of the ExB electron drift instability in Hall thrusters. Plasma Sourc Sci Technol (2019) 28(6):064002. doi:10.1088/1361-6595/ab08af

CrossRef Full Text | Google Scholar

21. Taccogna F, Minelli P, Bruno D, Longo S, Schneider R. Kinetic divertor modeling. Chem Phys (2012) 398:27–32. doi:10.1016/j.chemphys.2011.04.004

CrossRef Full Text | Google Scholar

22. Jambunathan R, Levin DA. A self-consistent open boundary condition for fully kinetic plasma thruster plume simulations. IEEE Trans Plasma Sci (2020) 48(3):610–30. doi:10.1109/tps.2020.2968887

CrossRef Full Text | Google Scholar

23. Li M, Merino M, Ahedo E, Haibin T. On electron boundary conditions in PIC plasma thruster plume simulations. Plasma Sourc Sci Technol (2019) 28(3):034004. doi:10.1088/1361-6595/ab0949

CrossRef Full Text | Google Scholar

24. Taccogna F, Garrigues L. Latest progress in Hall thrusters plasma modelling. Rev Mod Plasma Phys (2019) 3(1):12–63. doi:10.1007/s41614-019-0033-1

CrossRef Full Text | Google Scholar

25. Marín-Cebrián A, Domínguez-Vázquez A, Fajardo P, Ahedo E. Macroscopic plasma analysis from 1D-radial kinetic results of a Hall thruster discharge. Plasma Sourc Sci Technol (2021) 30(11):115011. doi:10.1088/1361-6595/ac325e

CrossRef Full Text | Google Scholar

26. Szabo J, Warner N, Martinez-Sanchez M, Batishchev O. Full particle-in-cell simulation methodology for axisymmetric Hall effect thrusters. J Propuls Power (2014) 30(1):197–208. doi:10.2514/1.B34774

CrossRef Full Text | Google Scholar

27. Taccogna F, Longo S, Capitelli M, Schneider R. Self-similarity in Hall plasma discharges: Applications to particle models. Phys Plasmas (2005) 12(5):053502. doi:10.1063/1.1877517

CrossRef Full Text | Google Scholar

28. Balay S, Abhyankar S, Adams M, Brown J, Brune P, Buschelman K, et al. PETSc users manual (2019).

Google Scholar

29. Villemant M, Belhaj M, Sarrailh P, Dadouch S, Garrigues L, Boniface C. Measurements of electron emission under electron impact on BN sample for incident electron energy between 10 eV and 1000 eV. Europhys Lett (2019) 127:23001. doi:10.1209/0295-5075/127/23001

CrossRef Full Text | Google Scholar

30. Thorsteinsson EG, Gudmundsson GT. A global (volume averaged) model of a nitrogen discharge: I. Steady state. Plasma Sourc Sci Technol (2009) 18:045001. doi:10.1088/0963-0252/18/4/045001

CrossRef Full Text | Google Scholar

31. Gudmundsson JT, Kawamura E, Lieberman MA. A benchmark study of a capacitively coupled oxygen discharge of the oopd1 particle-in-cell Monte Carlo code. Plasma Sourc Sci Technol (2013) 22:035011. doi:10.1088/0963-0252/22/3/035011

CrossRef Full Text | Google Scholar

32. Liebermann LA, Lichtenberg AJ. Principles of plasma discharges and materials processing. 2nd ed. Hoboken, NJ, USA: Wiley (2005).

Google Scholar

33. Brunger MJ, Buckman SJ. Electron–molecule scattering cross-sections. I. Experimental techniques and data for diatomic molecules. Phys Rep (2002) 357:215–458. doi:10.1016/S0370-1573(01)00032-1

CrossRef Full Text | Google Scholar

34. Itikawa Y. Cross sections for electron collisions with nitrogen molecules. J Phys Chem Ref Data (2006) 35:31–53. doi:10.1063/1.1937426

CrossRef Full Text | Google Scholar

35. Tabata T, Shirai T, Sataka M, Kubo H. Analytic cross sections for electron impact collisions with nitrogen molecules. Data Nucl Data Tables (2006) 92:375–406. doi:10.1016/j.adt.2006.02.002

CrossRef Full Text | Google Scholar

36. Tabata T, Shirai T, Sataka M, Kubo H. Erratum to “Analytic cross sections for electron impact collisions with nitrogen molecules” [At. Data Nucl. Data Tables 92 (2006) 375–406]. Data Nucl Data Tables (2011) 98:74. doi:10.1016/j.adt.2011.06.002

CrossRef Full Text | Google Scholar

37. Kawaguchi S, Takahashi K, Satoh K. Electron collision cross section set for N₂ and electron transport in N₂, N₂/He, and N₂/Ar. Plasma Sourc Sci Technol (2021) 30:035010. doi:10.1088/1361-6595/abe1d4

CrossRef Full Text | Google Scholar

38. Laporta V, Little DA, Celiberto R, Tennyson J. Electron-impact resonant vibrational excitation and dissociation processes involving vibrationally excited N₂ molecules. Plasma Sourc Sci Technol (2014) 23:065002. doi:10.1088/0963-0252/23/6/065002

CrossRef Full Text | Google Scholar

39. BSR database. BSR database (2022). Available from: www.lxcat.net (Accessed June 9, 2022).

Google Scholar

40. Lewis BR, Gibson ST, Zhang W, Lefebvre-Brion H, Robbe J-M. Predissociation mechanism for the lowest ¹Π_u states of N₂. J Chem Phys (2005) 122:144302. doi:10.1063/1.1869986

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Itikawa Y. Cross sections for electron collisions with oxygen molecules. J Phys Chem Ref Data (2009) 38:1–20. doi:10.1063/1.3025886

CrossRef Full Text | Google Scholar

42. Vahedi V, Surendra M. A Monte Carlo collision model for the particle-in-cell method: Applications to argon and oxygen discharges. Comput Phys Commun (1995) 87:179–98. doi:10.1016/0010-4655(94)00171-w

CrossRef Full Text | Google Scholar

43. Vass M, Wilczek S, Lafleur T, Brinkmann RP, Donkó Z, Schulze J. Electron power absorption in low pressure capacitively coupled electronegative oxygen radio frequency plasmas. Plasma Sourc Sci Technol (2020) 29:025019. doi:10.1088/1361-6595/ab5f27

CrossRef Full Text | Google Scholar

44. Biagi database. Biagi database (2022). Available from: www.lxcat.net (Accessed June 9, 2022).

Google Scholar

45. Itikawa Y, Ichimura A. Cross sections for collisions of electrons and photons with atomic oxygen. J Phys Chem Ref Data (1990) 19:637–51. doi:10.1063/1.555857

CrossRef Full Text | Google Scholar

46. Laher RR, Gilmore FR. Updated excitation and ionization cross sections for electron impact on atomic oxygen. J Phys Chem Ref Data (1990) 19:277–305. doi:10.1063/1.555872

CrossRef Full Text | Google Scholar

47. Bell KL, Gilbody HB, Hughes JG, Kingston AE, Smith FJ. Recommended data on the electron impact ionization of light atoms and ions. J Phys Chem Ref Data (1983) 12:891–916. doi:10.1063/1.555700

CrossRef Full Text | Google Scholar

48. Nanbu K. Probability theory of electron-molecule, ion-molecule, molecule-molecule, and Coulomb collisions for particle modeling of materials processing plasmas and cases. IEEE Trans Plasma Sci (2000) 28(3):971–90. doi:10.1109/27.887765

CrossRef Full Text | Google Scholar

49. Parodi P. Analysis and simulation of an intake for air-breathing electric propulsion systems. Graduation thesis. Pisa, Italy: Aerospace engineering department, University of Pisa (2019).

Google Scholar

50. Morozov AI, Melikov IV. Similitude in Hall-current plasma accelerators. Soviet Physics-Technical Phys (1974) 19:340–2.

Google Scholar

51. Yamashita Y, Tsukizaki R, Nishiyama K. Importance of stepwise ionization from the metastable state in electron cyclotron resonance ion thrusters. J Electr Propuls (2022) 1:2. doi:10.1007/s44205-022-00002-1

CrossRef Full Text | Google Scholar

52. Lawton SA, Phelps AV. Excitation of the b1S+ g state of O₂ by low energy electrons. J Chem Phys (1978) 69:1055. doi:10.1063/1.436700

CrossRef Full Text | Google Scholar