Plasma Dynamics and Electron Transport in a Hall-Thruster-Representative Configuration with Various Propellants: I—Variations with Discharge Voltage and Current Density

Abstract

The results from a wide-ranging parametric investigation into the behavior of the collisionless partially magnetized plasma discharge of three propellants—xenon, krypton, and argon—are reported in this two-part article. These studies are performed using high-fidelity reduced-order particle-in-cell (PIC) simulations in a 2D configuration that represents an axial–azimuthal cross-section of a Hall thruster. In this part I paper, we discuss the effects of discharge voltage and current density (mass flow rate). Our parametric studies assess the spectra of the resolved instabilities under various plasma conditions. We evaluate the ability of the relevant theories from the literature to explain the variations in the instabilities’ characteristics across the studied plasma parameter space and for various propellants. Moreover, we investigate the changes in the electrons’ cross-magnetic-field transport, as well as the significance of the contribution of different momentum terms to this phenomenon across the analyzed cases. In terms of salient observations, the ion acoustic instability (IAI)-related modes are found to be dominant across the simulation cases, with the ion transit time instability also seen to develop at low current density values. Across the explored parameter space, the instabilities have the main contributions to the electrons’ transport within the plume region. The peak of the electric momentum force term, representing the effect of the instabilities, overall shifts toward the plume as either the current density or the discharge voltage increases. The numerical findings are compared against relevant experimental observations reported in the literature.

1. Introduction

In electric plasma propulsion technologies like Hall thrusters, the global performance and stability of operation are majorly influenced by the underlying plasma phenomena. This unique correlation between engineering and physics in these technologies has posed challenges for their further growth and transformation. This is particularly because some of the important physical processes, such as the instability-induced electrons’ transport or the plasma–wall interactions that dynamically evolve under the influence of instabilities, are still not fully understood [1,2]. In the absence of a thorough understanding of the multidimensional, multiscale plasma phenomena that play a determining role in the global behavior of the thrusters, the design and development of plasma propulsion technologies have been so far mostly reliant on extended, resource-demanding testing to identify and/or mitigate the impact of the plasma phenomena across the operational envelope. Fully comprehending the less-understood plasma processes, their dependencies, effects, and interactions can lead to the development of first-principles-based predictive reduced-order models that alleviate the reliance on expensive experiments.

Toward the above, the numerical parametric studies provide an effective pathway, in part by enabling us to untangle the complex, multifaceted plasma behaviors. Such studies are also essential to complement the experimental research. The simulations provide a more controllable environment of study, in which the effects of various phenomena that normally affect the plasma behavior simultaneously can be isolated and investigated more closely. In addition, the simulations enable access to information regarding the underlying plasma processes that might not be readily accessible in an experimental setting due to the spatiotemporal resolution limitations of the diagnostics and/or the physical accessibility constraints related to the device being tested. As a result, performing parametric sweeps of different plasma discharge parameters using high-fidelity simulations can aid explaining the experimental observations and/or enriching the results of the test campaigns. The parametric simulations also serve as a viable approach to stress-test the available theoretical descriptions of the plasma to determine the extent to which the theoretically derived correlations between the characteristics of the plasma phenomena and the discharge properties are consistent with the numerical observations.

To investigate the variations in the characteristics of the underlying processes, as well as in the global plasma response, over large parameter spaces a numerical tool is necessary that enables self-consistent and cost-efficient multidimensional simulations of the microscopic, kinetic phenomena that influence the larger-scale, macroscopic behaviors. In this regard, even though the particle-in-cell (PIC) method is a suitable choice from the fidelity perspective, it is notably computationally expensive for multidimensional simulations without any scaling of the physical/geometrical parameters. This makes using traditional PIC codes for extensive simulations over large parameter spaces impractical. Moreover, the more computationally affordable fluid and hybrid fluid–PIC simulations today are unable to capture the kinetic effects, and they do not feature predictive capability. Indeed, despite the multitude of efforts over the past decade [3,4,5,6,7,8,9], generalizable closure models for the system of equations in the fluid and hybrid codes have yet to be established.

The reduced-order PIC scheme [10,11], developed by the authors at Imperial Plasma Propulsion Laboratory (IPPL) and extensively verified in different configurations [10,11,12,13,14,15], allows performing high-fidelity kinetic plasma simulations at a drastically reduced computational cost compared with conventional PIC methods. We emphasize that there have been other recent efforts as well in the community aimed at tackling the computational cost issue of PIC simulations, among which a few examples are reported in refs. [16,17,18,19]. Nonetheless, the reduced-order PIC scheme offers unique flexibility and extensibility characteristics that set it distinctly apart from the approaches pursued by other authors toward lowering the cost of PIC models.

Therefore, in this two-part article, we employed the reduced-order PIC approach to carry out a series of parametric investigations in a 2D configuration that represents an axial–azimuthal cross-section of a Hall thruster. The parameters whose influence we assessed in this work are the discharge voltage, total ion current density (equivalent to the mass flow rate), and the magnetic field topology in terms of its peak intensity along the axial direction as well as the axial gradients. The influence of all these parameters is investigated for three propellants—xenon, krypton, and argon.

The simulations are performed using the IPPL-Q2D code, which is based on the reduced-order PIC scheme. This code has already been used for several physics studies, including investigations of the parametric variation in plasma phenomena along the radial–azimuthal coordinates of a Hall-thruster-like geometry [20,21,22,23]. For the axial–azimuthal simulations in this article, the domain decomposition associated with the reduced-order PIC [10] is carried out using 40 vertical “regions” ($M$) along the axial direction and 20 horizontal regions ($N$) along the azimuthal direction. Prior verifications of this specific approximation of the 2D axial–azimuthal problem, reported in refs. [10,15], showed that it has a computational cost of only about 8% of a full 2D simulation, while featuring an L2 error of about 5–10% across plasma parameters. Thus, considering the physics-oriented scope of this effort and the extended parameter space to be investigated, the approximation order of the Q2D PIC simulations reported here provides a reasonable trade-off between the results’ accuracy and the computational cost. We elaborate further on these remarks within the rest of the article (Section 2 and Appendix A).

The objective of the presented parametric studies is to understand (1) variations in the macroscopic plasma properties, (2) changes in the characteristics of the plasma fluctuations, particularly the ion-acoustic and the ion transit time instabilities, and the extent to which the available theories in the literature can describe the observed changes, (3) variations in the electrons’ cross-field transport and the significance of the contribution of various momentum terms to this phenomenon, and (4) changes in the distribution function of the ion species.

In this part I of the article, the focus is on the effects of the discharge voltage and the ion current density (mass flow rate). The studies presented here build upon and expand the insights provided by a number of prior works from the literature.

The physics of the instabilities and the instability-induced electrons’ transport has been numerically investigated in 2D Hall-thruster-representative axial–azimuthal configurations in refs. [5,24,25,26,27,28,29,30,31,32,33,34,35]. Among these works, Boeuf and Garrigues [24] and Charoy [26,27] used a simplified simulation set-up that featured an imposed ionization source for parametric studies into the effects of the peak ion current density on the plasma discharge for the xenon propellant. Boeuf and Garrigues [24] presented a detailed numerical–theoretical analysis of the Electron Drift Instability (EDI) and its transition to ion-acoustic instability (IAI). They also observed the development of ion transit time instability at low current densities and argued that this is correlated with reduced electron transport at low current density conditions [24]. In the same problem set-up, Charoy [26] analyzed the contributions of the force terms in the electrons’ azimuthal momentum equation to the electrons’ transport for different current densities and compared the contribution of the instabilities against the theoretical relation of the enhanced transport via electron–ion friction force proposed by Lafleur et al. in ref. [4]. Petronio carried out a number of parametric studies in ref. [34] using an axial–azimuthal simulation set-up with a self-consistent description of the ionization and a 1D fluid treatment of the neutral species dynamics [36]. The author particularly investigated the effects of changes in the discharge voltage, the mass flow rate, and the propellants (xenon, krypton, and iodine) on the global performance parameters of the simulated Hall thruster and the characteristics of the breathing mode oscillations [34]. With a similar scope to the studies carried out by Petronio in ref. [34], Reza et al. performed self-consistent axial–azimuthal simulations of a high-power magnetically shielded Hall thruster using the single-region implementation [12,13] of the reduced-order PIC code over a few operating conditions of the thruster [35]. The studied operating conditions represented the variations in the mass flow rate and the magnetic field peak intensity at a fixed discharge voltage [35]. The impacts on the thruster’s performance, the characteristics of the breathing mode oscillations, the azimuthal instabilities, and the electrons’ mobility were investigated [35].

With respect to the above efforts, which are of most relevance to the present study, our work has a number of advances/novelties: First, we explored a broader parameter space in terms of the discharge voltage and the ion current density values, as well as the ion mass. Second, we extensively tested the relevant theories of plasma fluctuations with numerical observations concerning the characteristics of the resolved instabilities across the explored parameter space. Third, compared in particular against the works by Boeuf and Garrigues [24] and Charoy [26,27], which performed their parametric investigations with a cathode boundary condition based on the current-continuity approach [27], our parametric studies are carried out with a quasi-neutrality cathode condition [11]. This is an important distinction due to the documented effect that the cathode boundary condition has on the axial–azimuthal PIC simulation results [27]. Fourth and last, the spectral analyses of the plasma instabilities are complemented by the novel Dynamic Mode Decomposition approach [37,38], which enables a comprehensive simultaneous spatial-temporal characterization of the dominant instability modes.

2. Overview of the Reduced-Order PIC Scheme

As was pointed out in Section 1, the simulations discussed in this article are carried out using the IPPL-Q2D (quasi-2D) code, which is developed based on the “reduced-order” or “quasi-multidimensional” PIC scheme [10,11]. In this section, we present a brief description of the concept behind the method.

The reduced-order PIC relies on a novel dimensionality-reduction treatment that reduces the computational cost of conventional PIC simulations by up to three orders of magnitude, depending on the desired order of approximation of the problem and the number of involved dimensions.

The dimensionality reduction underlying the reduced-order PIC is enabled through the definition of a new computational grid system, which we term the “reduced-order grid”. The adoption of this grid system in PIC allows increasing the cell size beyond the Debye length scale along one or two simulation dimension(s) at a time (in 2D and 3D simulation, respectively) by an arbitrary amount while keeping the cell size along the other dimension(s) compatible with the ordinary PIC’s requirements. Consequently, the number of computational cells and, hence, the total number of macroparticles are reduced, resulting in lower computations. The schematics of the reduced-order grid in 2D, also known as the “reduced-dimension” or quasi-2D grid [14], is illustrated in Figure 1.

In the reduced-order grid system, the computational domain is decomposed into multiple rectangular (in 2D) or cubic (in 3D) “regions”, which can be thought of as the discretization of the domain using a “coarse” grid (Figure 1(left)). Within any specific region, denoted as $\mathsf{\Omega }$ on the left-hand side of Figure 1, we discretize each simulation dimension separately using 1D “elongated” cells (orange-shaded elements in Figure 1(right)) [14,22]. The criterion for the fine grids’ cell sizes along the $x$ and $y$ direction, $\mathsf{\Delta }x$ and $\mathsf{\Delta }y$, is the same as that in a conventional 1D PIC, i.e., smaller than the Debye length [14,22]. These fine-discretization grids yield a decoupling between the different coordinates in each region, allowing us to remarkably reduce the number of cells and, hence, the required total number of macroparticles [14,22]. The vertical and horizontal fine grids in each region (Figure 1(right)) capture separately the variations in the plasma properties along the $x$ and $y$ directions, respectively [14,22].

The described grid system establishes a new framework enabling quasi-multidimensional numerical simulations, which provide an approximation of the full multidimensional problem with arbitrary closeness to the corresponding full multidimensional solution. The order-reduction technique, by construction, provides flexibility for the user to determine the desired trade-off between the level of predictions’ accuracy and the computational cost. Indeed, the degree of “order reduction” (or approximation order) with respect to full 2D/3D and, consequently, the computational cost of a reduced-order simulation depends on the fineness of the coarse grid, i.e., the number of regions. In other words, the aspect ratios of the horizontal and the vertical cells, $l_x/\mathsf{\Delta }x$ and $l_y/\mathsf{\Delta }y$ are directly proportional to the speedup factor achieved through this reduced-order scheme.

The incorporation of the reduced-order grid in PIC leads to the development of the reduced-order PIC (RO-PIC) scheme. The resulting PIC code based on the RO-PIC scheme is a quasi-multidimensional code, which we refer to as quasi-2D (Q2D) when accounting for two dimensions. The reduced-order PIC overall follows the standard algorithm of the particle-in-cell method [39,40,41]. The distinction between the reduced-order PIC code (e.g., the Q2D code) and a conventional one lies in the computational grid used for solving the electric potential and the electric field as well as for scattering the particle-based data (such as the electric charge) onto and gathering the grid-based data (such as the electric field) from it onto the particles’ position. Therefore, the formulation and implementation of the functions and modules using and/or interacting with the grid must be adjusted, which includes the Poisson solver. The RO-PIC’s Poisson solve function resolves a “reduced-dimension” Poisson’s equation. The formulation of the reduced-dimension Poisson’s equation is reported in refs. [10,11,12,13], where the results of extensive verifications of the RO-PIC’s solution approach are also provided.

To conclude, in Appendix A, we summarized the results from prior verifications of IPPL-Q2D in the same axial–azimuthal case as that used for the studies here. The results include a comparison of the IPPL-Q2D’s predictions with various numbers of regions against the conventional 2D PIC results, as well as a plot of error convergence vs. the number of regions. Based on the results in Appendix A, the simulation with 20 regions along the azimuthal dimension and 40 regions along the axial dimension ($N=20$ and $M=40$) demonstrates a great agreement with the full 2D results, while yielding about 13-fold speedup compared with the full 2D PIC simulation.

Accordingly, as was also underlined in Section 1, the use of Q2D simulations with $N=20$ and $M=40$ for the presented studies in this work is warranted.

3. Description of the Simulations’ Set-Up and Conditions

The simulations are performed in a set-up that largely resembles that of the benchmarking effort in ref. [25]. This set-up represents a 2D Cartesian cross-sectional plane of a Hall thruster along the axial–azimuthal coordinates and is a well-established case in the community.

A schematic visualization of the simulation domain and the coordinate system is shown in Figure 2a, where $x$ and $z$ denote the axial and the azimuthal directions, respectively.

Referring to Figure 2a, the applied magnetic field is purely radial. The left-hand side boundary represents the anode side, whereas the cathode boundary is located on the right-hand side of the domain. The domain’s axial extent ($L_x$) is 2.5 cm, corresponding to 0.75 cm of a Hall thruster’s discharge channel and 1.75 cm of the near plume. The azimuthal length of the domain ($L_z$) is 1.28 cm. The number of nodes along $x$ and $z$ are by default 500 and 256, respectively. This corresponds to an equal cell size of $\mathsf{\Delta }x=\mathsf{\Delta }z$ = 50 $\mathsf{\mu }\mathrm{m}$. The default timestep of the simulations is $5\times {10}^{-12}$ s. These values of the cell size and the timestep are the same as those of the benchmark case [25]. For the simulations in which the peak ion number density and/or the peak electron temperature was observed to be larger than the corresponding values in the benchmark, the cell size and the timestep were adjusted accordingly to ensure that these simulation parameters meet the stability criteria of explicit PIC schemes [30,41]. The total simulated time in all cases is 30 $\mathsf{\mu }\mathrm{s}$, and the collisional processes are not considered.

From Figure 2b, the radial magnetic field ($B$) is noticed to have a Gaussian profile along the axial direction, similar to ref. [25], with its peak intensity of 10 mT occurring at the channel’s exit plane ($x$ = 0.75 cm). The set-up also features a temporally invariant ionization source ($S_i$) [25,26,27] with a cosine-shaped distribution along $x$ and a uniform distribution along $z$. As per the detailed descriptions available in refs. [25,26,27], the peak value of this ionization source is calculated such that the total ion current density from the system will be equal to a preset value of $J$. Varying the value of $J$ amounts to changing the operating mass flow rate of a Hall thruster in a real-world scenario. Note that a fixed current density condition would translate into different operating anode mass flow rates for various propellants. Figure 2b shows the profiles of the imposed ionization source for the simulations that we carried out to investigate the effects of the ion current density on plasma discharges of xenon, krypton, and argon. For these simulations, the value of $J$ is changed from 50 to 800 ${\mathrm{Am}}^{-2}$.

The curve with $J$ value of 400 ${\mathrm{Am}}^{-2}$ (yellow curve in Figure 2b) corresponds to the ionization source of the benchmark simulation case [25]. This is the profile of the source used in our other simulations dedicated to the study of the effects of the discharge voltage ($V$) in this part, and the effects of the magnetic field topology in part II [42].

To study the effects of varying $V$, we have swept this parameter over the range of 100 to 600 V. We refer to the $V$ value of 200 V as the baseline value since it is the same discharge voltage as that adopted in the benchmark [25]. This baseline $V$ value is applied between the anode and the cathode for various current density simulations in the present paper, as well as for the simulations in part II [42].

Table 1. List of the Q2D simulations performed with various propellants, summarizing the values of the total ion current density ($J$) and the discharge voltage ($V$) used in each case. $J_0$ and $V_0$ denote, respectively, the $J$ and the $V$ values of the benchmark case (baseline) from ref. [25].

Simulation Case	J Value [Am−2]	V Value [V]	Propellant
1	1/8 $J_0$ = 50	$V_0$	Xenon, Krypton, Argon
2	1/4 $J_0$ = 100	$V_0$
3	1/2 $J_0$ = 200	$V_0$
4 (Baseline)	$J_0$ = 400	$V_0$ = 200
5	2 $J_0$ = 800	$V_0$
6	$J_0$	1/2$V_0$ = 100
7	$J_0$	3/2$V_0$ = 300
8	$J_0$	2$V_0$ = 400
9	$J_0$	3$V_0$ = 600

summarizes the values of the ion current density and the discharge voltage used for various simulation cases reported in this paper.

To initialize the simulations, the electrons and the ions are loaded uniformly throughout the simulation domain as electron–ion pairs with a density of $5\times {10}^{16} {\mathrm{m}}^{-3}$. The plasma species are sampled from Maxwellian distribution functions at 10 eV for the electrons and 0.5 eV for the ions. The simulations are initiated with a total macroparticle count of 1,024,000, which corresponds to a macroparticle-per-cell count of 100.

Regarding the particles’ boundary conditions, both the ions and the electrons that reach either the anode or the cathode side are removed from the simulation. However, to maintain the discharge, electrons at each timestep are injected into the domain from the cathode side. The number of the injected electrons is obtained using the quasi-neutrality approach, as described in refs. [11,27]. We chose to follow the quasi-neutrality approach for the electrons’ reinjection because ref. [27] demonstrated that this approach is less sensitive to the temperature of the reinjected electrons. The reinjected electrons are sampled from a Maxwellian distribution at the electrons’ initial temperature of 10 eV. The azimuthal direction is periodic. Hence, the particles leaving the domain at one end are injected back from the other end.

As for the electric potential boundary conditions, along the azimuthal direction, a periodic boundary condition is considered. Along the axial direction, Dirichlet boundary conditions of $V$ volts at the anode side (see Table 1 for the exact value across the different simulations) and 0 V at the cathode side of the domain are applied. These boundary conditions are implemented in the reduced-order PIC, as per the approach described in refs. [10,11].

Regarding the set-up of the problem, we would emphasize that it represents a simplification to a real-world Hall thruster in several aspects.

First, the ionization is imposed through a fixed injection source rather than being captured self-consistently. As a result, the low-frequency discharge oscillations associated with breathing mode [2,43] are neither represented nor are the focus of this work. Moreover, it is noted that, in a real setting, the strength of the ionization and the relative positioning of the ionization zone with respect to the acceleration region change over a breathing cycle [43]. These variations influence the instantaneous plasma properties (conditions), as well as the high-frequency instabilities [2]. Although these effects are not captured within individual simulations here, variations across the simulations with different plasma conditions can cast light on some of these effects. In other words, the conditions established within each simulation can represent the plasma state at a different stage of the discharge’s breathing cycle. The quasi-static nature of the present set-up further enables the effective isolation of various effects.

Second, the collisions are not accounted for in our simulations. Thus, the electrons’ classic collisional scattering is absent. Neglecting collisions additionally implies that any difference observed in the plasma behavior with various propellants is due to the difference in the ion mass. This is in fact desired in this case since our intention has been to distinguish the direct potential effects of the ion mass from the effects of collisional processes (which are due to differences in the collision cross-sections of the propellants). Moreover, the absence of collisions facilitates comparing the numerical results of the instabilities’ characterizations against the published theories.

Third and last, because of the absence of some of the energy loss mechanisms, such as ionization and excitation as well as wall losses, the predicted temperatures are expected to be overestimated. Also, given that the simulations do not account for the radial direction and wall interactions, the term “exit plane”, which is referenced throughout the text and in the captions of the figures, merely represents the location of the peak of the magnetic field profile and does not have any physical significance in terms of altering boundary conditions. Subsequently, the terms inside the “channel” and in the “plume” refer to the axial locations before and after the B-field’s peak, respectively. Finally, the near-anode region (near the gas distributor element in a typical Hall thruster) is not simulated, and the magnetic field values are relatively high near the anode-side (left) boundary of the domain.

4. Brief Notes from the Theory of the Electron Drift and the Ion-Acoustic Instabilities

Along the axial–azimuthal plane of Hall thrusters and similar plasma technologies, which feature a perpendicular orientation of the applied electric (E) and magnetic (B) fields—the so-called “E × B” configuration, the Electron Cyclotron Drift Instability (ECDI) [28,44] and the ion-acoustic instability [45,46] are shown to develop as dominant azimuthal fluctuation modes. This section provides a brief overview of the theoretical background on these instabilities.

The solution of the kinetic dispersion relation of ECDI [47,48] in a plane perpendicular to the magnetic field (axial–azimuthal plane) indicates the presence of unstable modes near the electron cyclotron resonances ($k_{z, res}=n\frac{{\omega }_{ce}}{V_{de}} ; n=1, 2,\dots$). The resonance broadening caused by the wave–particles interactions will eventually smooth out the cyclotron resonance, leading to the transition of the ECDI to the ion-acoustic instability [45,46]. The frequency ($\omega$) and the growth rate ($\gamma$) of the ion-acoustic instability are given by Equations (1) and (2) [3,24,45,46]

$$\omega \approx k_x{V}_{di}\pm {\displaystyle \frac{k{\lambda }_{De}{\omega }_{pi}}{\sqrt{1+k^2{\lambda }_{De}^2}}} ,$$

(1)

$$\gamma \approx \pm \sqrt{{\displaystyle \frac{\pi m_e}{8m_i}}}{\displaystyle \frac{k_z{V}_{de}}{{\left(1+k^2{\lambda }_{De}^2\right)}^{\frac{3}{2}}}} ,$$

(2)

where ${\lambda }_{De}$ is the Debye length, ${\omega }_{pi}$ is the ion plasma frequency, $k_x$ and $k_z$ are the axial and azimuthal wavenumber components, respectively. $V_{di}$ is the ions’ drift velocity, $V_{de}$ electrons’ drift velocity, and $m_e$ and $m_i$ are the electron and ion particle mass, respectively. It has been observed that the axial wavenumber of the instability is significantly smaller than the azimuthal wave number ($k_x\ll k_z$) [12,26,32]. Thus, in Equations (1) and (2), we can assume $k\approx k_z$.

The dominant mode of the IAI corresponding to the mode with the maximum growth rate satisfies ${\displaystyle \frac{\partial \gamma }{\partial k_z}}=0$ condition. Accordingly, the characteristics of the instability’s fastest growing mode, namely the azimuthal wavenumber ($k_{z, max}$), the frequency (${\omega }_{max}$), and the phase velocity ($v_{ph,max}$), are obtained as

$$k_{z, max}\approx {\displaystyle \frac{1}{\sqrt{2} {\lambda }_{De}}} , {\omega }_{max}\approx {\displaystyle \frac{{\omega }_{pi}}{\sqrt{3}}} , v_{ph,max}\approx \sqrt{{\displaystyle \frac{2}{3}}}{C}_s ,$$

(3)

where C_s is the ion sound speed. Ion-trapping is proposed as the saturation mechanism for the IAI [45,46,49]. The instability grows quasi-linearly until it reaches a sufficiently large amplitude that can trap significant population of the ions. At the saturation state, the approximate RMS amplitudes of the azimuthal electric field ($\delta E_{RMS}$) and electron number density ($\delta n_{e, RMS}$) fluctuations corresponding to the fastest growing mode of IAI are given by the expressions in Equation (4) [24,30,50]

$$\delta E_{RMS}\approx {\displaystyle \frac{T_e}{12{\lambda }_{De}}} , \delta n_{e, RMS}\approx {\displaystyle \frac{n_e}{6\sqrt{2}}} ,$$

(4)

where $T_e$ and $n_e$ represent the steady-state electron temperature and number density, respectively.

5. Results

In this section, we present the results of the Q2D simulations to show how the discharge voltage and the total ion current density affect the time-averaged macroscopic plasma properties, as well as the electron and the ion currents [12,13]. The results are shown for xenon, krypton, and argon gases to additionally demonstrate the influences of the ion mass.

All results discussed in this section and the subsequent one correspond to the quasi-steady state of the system. The simulations reach a quasi-steady state (characterized by almost periodic oscillations) within less than $10 \mathsf{\mu }\mathrm{s}$. Hence, the time interval of 20–30 $\mathsf{\mu }\mathrm{s}$ is taken as the window for the analysis of data.

Figure 3 and Figure 4 show, respectively, the time-averaged (over 20–30 $\mathsf{\mu }\mathrm{s}$) axial profiles of plasma properties—ion number density ($n_i$), axial electric field ($E_x$), and electron temperature ($T_e$)—for the simulation cases with different discharge voltages ($V$) and various total ion current densities ($J$).

Overall, we observe that the typical features of an E×B discharge along the axial–azimuthal plane are represented across all the simulated conditions. These include the occurrence of the maximum number density around the center of the ionization zone upstream of the location of the magnetic field’s peak. The ionization zone is defined by the axial extent of the ionization source and is highlighted by the grey areas in the figures. Moreover, for all plasma conditions, the peak of the axial electric field occurs slightly downstream of the magnetic field’s peak, and the acceleration region (which correspond to the region of high $E_x$ values) has an overlap with the ionization zone. A significant portion of the acceleration occurs after the B-field’s peak, i.e., outside the channel (beyond the vertical dashed black line). Also, it is noticed that the peak of the electron temperature nearly aligns with the peak of the applied magnetic field in all the cases.

We now examine more closely the trends in Figure 3 and Figure 4. We start with the variations in the plasma properties with the discharge voltage (V). Increasing $V$ from 100 V to 600 V expectedly leads to a consistent increase in the $E_x$ profile. This accelerates the ions to higher velocities, thereby causing a monotonic decrease in the overall $n_i$ profile (See also Figure 5, left column). The rate of increase in peak $E_x$ with voltage seems to decrease at high voltages (from $V=400$ to $600 \mathrm{V}$). This implies the reduced effectiveness of the magnetic barrier in confining the axial motion of the electrons, possibly due to higher electron temperatures at high voltages. In particular, at $V=600 \mathrm{V}$, the reduced electron confinement is also evidenced by the sudden increase in the electron axial current reaching the anode boundary. This will be shown shortly in this section in a figure presenting the currents’ variation trends. In practice, the magnetic field intensity is typically increased for operation at higher voltages to ensure effective electron confinement [51,52].

Moreover, as $V$ increases, the width of the acceleration zone (the region of high $E_x$ magnitude) progressively expands, and the location of the peak of the $E_x$ profile shifts downstream. The elevated axial electric fields associated with higher voltages generate a more pronounced Joule heating effect [27], which leads to higher electron temperatures. These trends are similar among the three studied propellants.

Experimental measurements have shown that increasing voltage leads to the broadening of and the upstream shift in the acceleration zone [53,54,55]. Gawron et al. [56], however, observed that the spatial extent of the high electric field region decreases when $V$ is increased. In ref. [55], the authors emphasized that the obtained $E_x$ profile from LIF measurements of the ions’ axial velocity is quite sensitive to the data postprocessing procedure. Additionally, ref. [57] argued that the extent of the acceleration layer is subject to misinterpretation, as it depends on the chosen definition of such extent, leading to previous data being erroneously interpreted as a broadening of the acceleration region. Accordingly, they observed that in the PPSX000 thruster [57], the width of the acceleration layer remains unchanged with increasing $V$, whereas, in the PPS100-LM thruster [57], the acceleration region narrows down when $V$ increases. Ref. [57], however, reported that the peak of the $E_x$ profile shifts upstream toward the anode, consistent with the other measurements [53,54,55,56]. It was speculated that the upstream movement of the acceleration zone (layer) is linked to the increased electron temperatures at higher voltages [57]. At higher temperatures, electrons have a larger Larmor radius, leading to less confinement by the magnetic field. As a result, they penetrate deeper into the magnetic barrier, causing the electric field to shift towards the anode. Ref. [53], nonetheless, hypothesized that the upstream shift of the acceleration zone is because of the narrowing down of the ionization zone due to the higher temperatures at high voltage values.

Noting the above, we observe from Figure 3 that, except for the voltage range of $100$ to $200 \mathrm{V}$, the peak of the $E_x$ profile shifts slightly downstream as the voltage increases. It is underlined that, here, we consider the variation in the axial location of the peak $E_x$ as a measure of the upstream/downstream movement of the acceleration zone. The opposite trend to the experiments is most likely due to the fixed ionization source in our simulations. Given that, in the simulations, the electrons’ transport is self-consistently resolved, observing an opposite trend compared with the experiments undermines the increased electron temperature (and, thus, less efficient confinement) as a plausible explanation for the upstream shift of the peak of $E_x$. Instead, this upstream shift must be attributed to the ionization-related processes.

We now turn our focus to the variation trends with respect to the imposed ion current density ($J$). It is noticed from Figure 4 that the profiles of the $n_i$ and the $T_e$, as well as the peak of the $E_x$, rise monotonically with increasing $J$, albeit with much lower rates at low $J$ values ($J=50 \mathrm{to} 200 {\mathrm{Am}}^{-2}$). The increase in the $E_x$ peak is the consequence of Generalized Ohm’s law [52], which suggests a direct relationship between the current density and the established electric field. Later in this section, we will show that the axial electron current density consistently increases with $J$. Nonetheless, as will be demonstrated in Section 6.1 and Section 6.2, the higher plasma density at larger Js leads to the development of stronger azimuthal plasma fluctuations, which induce larger axial electron transport. Although the enhanced transport (equivalently, a larger plasma conductivity or a weaker electron confinement by the magnetic field) contributes to the reduction in the established electric field, this reduction does not compensate for the increase in the electrons’ axial current. Accordingly, the net result of these two opposing factors is the rise in the $E_x$ peak values. Furthermore, it is noted that the acceleration zone becomes narrower as the $J$ increases, reducing the overlap between the acceleration and ionization zones. The variations in the peak of $E_x$ across the simulated Js progressively decrease from heavier to lighter propellants (from xenon to argon). The elevated electron temperature at higher $Js$ is attributed to the increased Joule heating [27] as well as the enhanced particle heating caused by stronger instability waves.

Most of the above observations are consistent with the LIF measurements of Hall thruster plasma. The experiments indeed showed that when the anode mass flow rate (equivalent to $J$ here) increases, the peak intensity of the $E_x$ increases, and the width of the zone of high electric field magnitudes tends to become narrower [54,57]. However, ref. [55] observed no modification in the shape of the $E_x$ with the mass flow rate. Regarding the position of the region of high $E_x$ values, Hargus et al. and Mazouffre et al. [54,57] noticed an upstream movement whereas, in ref. [55], a downstream shift was observed for the $E_x$ profile with increasing mass flow rate. Such a nonunique trend might be due to the different operating conditions tested. In fact, our simulations also show a nonmonotonic behavior. At lower Js ($J=50 \mathrm{to} 200 {\mathrm{Am}}^{-2}$), the peak $E_x$ moves upstream, whereas at higher $J$s ($J=400 \mathrm{to} 800 {\mathrm{Am}}^{-2}$), the peak $E_x$ shifts downstream (this can be better seen from the last figure of Section 6.1).

Moving on, Figure 5 shows the variations in the spatiotemporally averaged ion number density and ions’ axial drift velocity ($V_{ix}$) in the plume ($1.5 \mathrm{cm}<x<2.5 \mathrm{cm}$). This extent corresponds to where the ions are fully accelerated, and the axial profiles of their properties are almost flat. The ideal exhaust velocities ($v_{exhaust}$) of the ions, calculated using Equation (5), are superimposed on the $V_{ix}$ plots

$$v_{exhaust}=\sqrt{{\displaystyle \frac{2e{V}_d}{m_i}}} .$$

(5)

In Equation (5), $e$ is the elementary charge. As is seen on the right column of Figure 5, the $V_{ix}$ values from the simulations deviate from the ideal theoretical trends. This discrepancy is due to the fact that Equation (5) does not account for energy losses. In reality, the voltage in Equation (5) shall be the effective voltage available for the acceleration of the ions, i.e., the beam voltage $V_b={\eta }_e{V}_d$, with ${\eta }_e$ being the energy efficiency [52]. The main energy loss mechanisms, which can subtract from the beam voltage include the ionization (and excitation), wall losses, and voltage losses at the cathode (non-zero cathode coupling voltage). However, considering the absence of these loss mechanisms in this problem, the main source of acceleration inefficiency is the overlap between the acceleration and the ionization zones. This overlap causes the ions that are injected within the acceleration region to accelerate only through a partial potential drop, occurring downstream of their injection location, which leads to velocities below the nominal velocities estimated from Equation (5). However, there is a counter effect related to the sheath formed at the anode side of the domain, which can increase the potential drop through which the ions are accelerated. This is manifested as a bump followed by a dip in the $E_x$ profiles in Figure 3 and Figure 4. The interplay between these two opposing factors more completely explains the discrepancy between the ions’ axial velocities from the simulation and those estimated from Equation (5).

Particularly, the deviation of the ion velocities from the ideal values ($v_{exhaust}$) consistently increases toward higher discharge voltages. This is because, at lower voltages, the potential rise within the plasma bulk due to the anode sheath can almost compensate for (or, e.g., exceed at $V=100 \mathrm{V}$ with xenon) the acceleration inefficiency due to the overlap of ionization and acceleration zones. However, as $V$ increases, the effect of the anode sheath becomes less significant compared with the discharge potential drop, hence, failing to balance the losses due to the acceleration inefficiencies. For the same reason, the deviation of the $V_{ix}$ from the $v_{exhaust}$ becomes more substantial from xenon to argon due to the consistent decrease in the anode sheath potential drop from xenon to argon, as is seen in Figure 3 (middle row). Despite these observations, experimental measurements [54,57] indicated a constant voltage loss (about 60–75 V) across all tested voltages, suggesting the independence of the voltage losses from the operating voltage. This discrepancy between our numerical results and the experiments in terms of the voltage loss trend is most likely due to not accounting for all the loss mechanisms in the simulations, including the ionization, the wall losses, and the cathode coupling voltage.

Concerning the variations in $V_{ix}$ with $J$, we observe from Figure 5 that the discrepancy from ideal values diminishes from $J=50 {\mathrm{Am}}^{-2}$ to $J=400 {\mathrm{Am}}^{-2}$ because of the narrowing of the acceleration region, thereby a reduced overlap with the ionization zone, as well as the increase in the anode sheath potential drop, which compensates for increasingly larger portions of the acceleration nonidealities. Similar to the cases with various voltages, the discrepancy between the simulated and the ideal ions’ velocities becomes larger from xenon to krypton, which is attributed to the progressive decrease in the sheath potential from xenon to krypton (middle row of Figure 4). At $J=800 {\mathrm{Am}}^{-2}$ with xenon and krypton; however, the anode sheath potential is so large that it overcompensates the acceleration inefficiencies, leading to ions’ velocities beyond the ideal values. In the same condition with argon, the two effects just described are almost balanced, and the $V_{ix}$ is thus close to the $v_{exhaust}$.

In Figure 6, we illustrated the variation trends of the electron current ($I_e$), the ion beam current ($I_i$), and the ratio of these two current terms ($I_e/I_i$) against the simulated plasma conditions. Note that the $I_e/I_i$ ratio represents the electron cross-field transport and, hence, is a measure of the inefficacy of the magnetic field barrier in confining the electrons, which in turn impacts the efficiency of the Hall thruster.

The ion current has been obtained from the simulations as the flux of the ions crossing the cathode boundary of the domain (${\Gamma }_{\mathrm{ai}}$). The electron current is calculated by subtracting the ion current from the discharge current, $I_d$ ($I_e=I_d-I_i$), where $I_d$ is determined as the difference between the fluxes of the electrons and the ions reaching the anode boundary (${\Gamma }_{ae}-{\Gamma }_{\mathrm{ai}}$).

The top row in Figure 6 shows that, while $I_e$ remains almost constant for xenon across the voltages of $200$ to $400 \mathrm{V}$, krypton and argon exhibit greater variability within this voltage range. From $400$ to $600 \mathrm{V}$, $I_e$ jumps by a factor of about 1.5 which, as mentioned before, is due to the weaker electron confinement by the fixed magnetic field. $I_i$ appears to reach a plateau beyond $V=300 \mathrm{V}$, indicating the maximum ion beam current that can be extracted with the imposed injection source intensity.

From the bottom row plots in Figure 6, it is observed that, as $J$ increases, the $I_e/I_0$ and the $I_e/I_i$ ratios both become larger. This is mainly due to the enhanced electrons’ axial transport induced by stronger azimuthal instabilities at higher $J$ values. It is highlighted though that the unrealistically high electron currents at elevated $J$ conditions are due to the fact that the baseline magnetic field of the set-up does not necessarily correspond to an optimal Hall thruster design.

We additionally note that the slopes of the variations in $I_e/I_0$ and $I_e/I_i$ are different among various propellants. Another noteworthy observation is that, at low $J$ values of 50 and 100 ${\mathrm{Am}}^{-2}$, the $I_i$/$I_0$ is near unity, meaning that most of the ions introduced into the domain by the injection source are accelerated downstream toward the cathode and are extracted at the cathode side. However, toward larger $J$s, the ratio of $I_i/I_0$ becomes increasingly less than 1 for all the studied propellants, implying that an increasing fraction of the ions drift toward the anode. A similar observation has also been reported in ref. [24] from the xenon-propellant axial–azimuthal simulations.

6. Discussion

More detailed analyses of the simulations’ results are presented and discussed in this section.

In Section 6.1, we analyze how the variations in the discharge voltage and the ion current density affect the axial profiles of the electrons’ transport and the contributions of the different force terms in the azimuthal momentum equation of the electron species.

In Section 6.2, we examine the influences of $V$ and $J$ on the characteristics of the resolved instabilities from the simulations. The instabilities’ spatiotemporal behavior is primarily investigated by analyzing the dispersion map of the instabilities in addition to using the Dynamic Mode Decomposition (DMD) approach. We have recently demonstrated in refs. [37,38] that the DMD (Optimized DMD variant to be precise) is a powerful spectral analysis technique for plasma systems in lieu of or in complementarity to the fast Fourier transform (FFT) method. Where relevant, the variations with the plasma conditions ($V$ or $J$) of the wavenumber, the frequency, and the RMS amplitude of the fluctuations are compared against the theories overviewed in Section 4.

Finally, in Section 6.3, we discuss how the changes in the $V$ and $J$ affect the velocity distribution function of the ion species.

6.1. Variations in the Electrons’ Transport and the Contributing Momentum Force Terms

In Figure 7, we compared the time-averaged axial profiles of the force terms in the azimuthal momentum equation of the electrons (Equation (6)) across the different simulations carried out.

$$-q{n}_e{v}_{e,x}B={\partial }_t\left(m{n}_e{v}_{e,z}\right)+{\partial }_x\left(m{n}_e{v}_{e,x}{v}_{e,z}\right)+{\partial }_x\left({\Pi }_{e,xz}\right)+q{\tilde{n}}_e{\tilde{E}}_z$$

(6)

In Equation (6), $q$ is the unit charge, $m$ is the electron mass, $v_{e,x}$ the axial drift velocity, $v_{e,z}$ the azimuthal drift velocity, ${\tilde{n}}_e$ the number density fluctuations, and ${\tilde{E}}_z$ the azimuthal electric field fluctuations. Each term on the right-hand side of Equation (6) represents a mechanism affecting the electrons’ mobility toward the anode, which, following the convention from refs. [12,27,30], is defined as follows: $F_t={\partial }_t\left(m{n}_e{v}_{e,z}\right)$ is the temporal inertia, $F_I={\partial }_x\left(m{n}_e{v}_{e,x}{v}_{e,z}\right)$ is the convective inertia, $F_{\Pi }={\partial }_x\left({\Pi }_{e,xz}\right)$ is the viscous force term, and $F_E=q{\tilde{n}}_e{\tilde{E}}_z$ is the electric force term. The aggregate effects of these mechanisms on the electrons’ azimuthal momentum are represented by the magnetic force term, $F_B=-q{n}_e{v}_{e,x}B$. The definition of all the force terms in Equation (6) as the moments of the electrons’ distribution function can be found in ref. [12]. The $F_t$ term was found to be always vanishingly small; therefore, it is not discussed here.

Referring to the plots in Figure 7, a number of common trends and behaviors can be noticed across most of the simulation cases:

(i)

The $F_{\Pi }$ term (orange curve) becomes a dominant contributor to the electrons’ transport near the anode side of the domain, whereas, near the location of the magnetic field peak almost immediately after the exit plane, the $F_{\Pi }$ values become negative.

(ii)

The $F_I$ term exhibits an almost symmetrical behavior with respect to the axial distribution of the $F_{\Pi }$ term. The $F_I$ profile peaks near the exit plane and becomes negative toward the anode side. Away from the near anode, the $F_I$ values are nearly the negative counterparts of $F_{\Pi }$ values such that they almost cancel each other out.

(iii)

Apart from the near-anode region, in most cases, the $F_E$ term has the dominant contribution to the electrons’ transport, especially from near the magnetic field peak region toward the plume.

Another interesting aspect concerns the variation in the axial locations of the peak values of the $F_B$ and the $F_E$ terms denoted as ${x_F}_{B,max}$ and ${x_F}_{E,max}$, respectively. These trends are visualized in Figure 8. As a reference, the axial locations of the peak intensity of the axial electric field (${x_E}_{x,max}$) are also shown in this figure.

Referring to Figure 7 and Figure 8, the following specific observations can be made regarding the variation trends of the momentum force terms with the discharge voltage and the total ion current density.

Across various voltages ($\mathit{V}$): the values of the $F_B$ term remain almost constant in the $V$ range of $100$ to $400 \mathrm{V}$. However, for krypton, some variability in the peak value of the $F_B$ is noticed. At $V=600 \mathrm{V}$, the $F_B$ values notably increase for all the propellants. The profile of the $F_{\Pi }$ term develops an additional peak near the exit plane with increasing $V$. This increases the extent over which the $F_{\Pi }$ term is dominant. Furthermore, as $V$ increases, the contributions of the $F_{\Pi }$ and the $F_{\mathrm{I}}$ terms become increasingly comparable with that of $F_E$ term. This observation can be traced to the generally larger electron drift and thermal velocities at higher voltages. Nevertheless, due to the symmetric behavior of the $F_{\Pi }$ and the $F_{\mathrm{I}}$ terms, their contributions almost cancel out upstream the channel exit plane, leaving $F_E$ as a primary determining factor in electrons’ transport. Finally, whereas the $x_{F_{B,max}}$ remains relatively unchanged across various voltages, the $x_{F_{E,max}}$ consistently shifts toward downstream for all the propellants as $V$ increases.

Across various total ion current densities ($\mathit{J}$): the behaviors of the different force terms are overall similar across $J$ values. Increasing $J$ mostly results in the magnitude of the terms to become almost proportionally larger. Moreover, the relative contribution of each term to the overall electrons’ transport at each axial location remains almost the same across various $J$s. Lastly, both the $x_{F_{B,max}}$ and the $x_{F_{E,max}}$ overall exhibit consistent downstream shifts with increasing $J$.

In addition to the above observations, a “wavy” axial distribution is noted for the $F_E$ and the $F_I$ terms in the plume at $V=100 \mathrm{V}$ and $J=800 {\mathrm{Am}}^{-2}$ for the xenon propellant, the exact root cause of which is not yet known.

As a final remark here, the relative “symmetry” between the viscous term ($F_{\Pi }$) and the convective inertia term ($F_I$) is an interesting observation that merits further investigation. This can be important for the plasma fluid simulations because if only one of the effects (and not both) is considered within the simulation, its contribution to electrons’ transport may remain unbalanced, potentially leading to inaccuracies in the results. However, it is important to highlight that the assumption of a prescribed ionization source may alter the physics of the problem and, consequently, affect the transport terms. Hence, a deeper assessment needs to be performed in the presence of self-consistent ionization and collisional processes regarding the axial distributions of the momentum force terms $F_{\Pi }$ and $F_I$.

6.2. Numerical and Theoretical Analyses of the Instabilities’ Characteristics

To set the context for the following discussions, we recall that one of the important instabilities developing along the azimuthal direction of a Hall thruster (and similar E×B plasma configurations) is the ECDI, which is a short-wavelength, high-frequency instability that features high growth rates [24,28,58,59,60]. Previous studies in the literature [25,26,27], as well as in our preceding investigations [12,13], have shown that, in the adopted axial–azimuthal set-up of the simulations in this work, the rapid nonlinear growth of the ECDI cannot be captured. Instead, the simulations resolve the pseudo-saturated state of the ECDI, which has been reported in the literature [45,46] to be associated with a transition to relatively longer-wavelength IAI modes of lower growth rates [61]. Hence, the dominant azimuthal wave modes from our simulations are expected to represent ion-acoustic-like characteristics.

To assess the validity of this statement across the simulated conditions, we compared the numerical dispersion maps of the azimuthal instabilities (obtained from the simulations) with the theoretical dispersion relation of the IAI (given by Equation (1)). These comparisons are provided for the xenon propellant as an example in Figure 9. We observe from this figure that, for most cases, the numerical dispersion of the instabilities within the channel (at $x=0.375 \mathrm{cm}$) aligns with the IAI’s nonlinear dispersion relation. However, in the plume (at $x=1.625 \mathrm{cm}$), the instabilities’ dispersion matches rather well the linear dispersion relation associated with the ion sound waves ($\omega =k_z{C}_s$). In addition, the instabilities exhibit dispersion across a broader range of frequencies and azimuthal wavenumbers in the plume compared with inside the channel.

Nonetheless, we note that there are some exceptions to the observed behaviors described above. In particular, at very low current density ($J=50 {\mathrm{Am}}^{-2}$), the instabilities inside the channel display significant deviation from the IAI’s dispersion relation.

Having observed an overall agreement between the instabilities’ numerical dispersion plots and the IAI’s theoretical dispersion relation (Equation (1)), we now aim to test the efficacy of the existing theories on the IAI in estimating the characteristics of the numerically observed dominant wave modes across the studied set-up conditions. To this end, we have compared the predicted characteristics of the azimuthal instabilities by the simulations with their theoretical estimations at two axial locations: (1) at $x=0.375 \mathrm{cm}$, which corresponds to the middle point of the discharge channel along the axis, and (2) at $x=1.625 \mathrm{cm}$, which is situated midway within the plume region.

In Figure 10, the azimuthal wavenumber ($k_z$) and the real frequency ($f$) of the dominant IAI wave mode from the simulations at the two probed axial locations are plotted. Note that, in Figure 10, the frequencies are in Hertz (Hz, ${\mathrm{s}}^{-1}$) and not in radians per second (rad/s). The simulated $k_z$ and $f$ of the dominant IAI mode are determined by identifying the highest peaks in the spatial and temporal FFTs of the azimuthal electric field signal [22,23]. The theoretical values, which are associated with the fastest-growing IAI mode, are calculated from Equation (3), with the ${\lambda }_{De}$ and the ${\omega }_{pi}$ in the relations of Equation (3) estimated using the time-averaged plasma properties at the respective axial locations from the simulations.

From Figure 10, it is evident that, in general, the characteristics of the instabilities vary from inside the channel toward the plume. A possible explanation for this variation is the downstream convection (with the ion beam) of the waves excited in the region near the anode boundary. Given the significant variation in the plasma conditions from the anode to the plume regions, the instabilities’ dispersion relation changes in accordance with the local plasma conditions as the waves are convected out with the ion beam [26].

Concerning the variations vs. $V$, we see from Figure 10 that the theoretical relations predict a monotonic decrease in the $k_z$ and the $f$ of the fastest-growing IAI mode with $V$, which is a consequence of the trends of the time-averaged plasma properties vs. $V$ as were seen in Figure 3. However, it is noticed that the simulated wavenumber and frequency of the dominant wave mode do not exhibit a similar trend with $V$. At low voltages, the simulated $k_z$ values are quite disparate from the theoretical ones, both within the channel and in the near-plume whereas, at higher voltages, the agreement with the theory is seen to relatively improve. In the frequency plots, we notice that, within the channel, the agreement with the theory is rather acceptable. However, in the plume, the theoretically predicted frequency of the dominant IAI mode is consistently lower than the simulated value across all the voltages and the propellants, except for the xenon propellant at $200 \mathrm{V}$ (the baseline case), where the theory and the simulation are in close agreement.

The theory and the simulations agree on the overall increasing trend of the $k_z$ and the $f$ with $J$, both within the channel and in the plume. Nevertheless, their values do not always exactly match. The most agreement is observed between the frequency values inside the channel up to $J=400 {\mathrm{Am}}^{-2}$. Also, for the xenon propellant in the plume, the theoretical frequencies are quite close to the simulated ones for $50{ \mathrm{Am}}^{-2}\le J\le 400 {\mathrm{Am}}^{-2}$. It is interesting that, in the plume and within the ion current density range of 50–400 ${\mathrm{Am}}^{-2}$, the wavenumber of the dominant mode for almost all the three propellants is rather constant with $J$. Also, for $J\le 400{ \mathrm{Am}}^{-2}$, $k_z$ either increases or remains almost unchanged from channel to the plume. However, at $J=800{ \mathrm{Am}}^{-2}$ for xenon and krypton, the $k_z$ is notably smaller in the plume compared with inside the channel.

Regarding the saturation of the azimuthal waves, to verify that the trapping of the ions by the waves occurred in our simulations, we look at the ions’ azimuthal phase-space plots in Figure 11, shown for the xenon propellant as an example. We observe that except for the conditions of low current density ($J\le 100 {\mathrm{Am}}^2$), the ions’ distribution is strongly affected by the azimuthal waves, featuring characteristic loops in $z-v_z$ space that are indicative of the ion-wave trapping. At $J=50$ and $100 {\mathrm{Am}}^{-2}$, the amplitude of the waves is not strong enough to substantially impact the ions’ distribution function.

In Figure 12, we compared the RMS amplitudes of the azimuthal electric field ($\delta E_{RMS}$) and the electron number density ($\delta n_{e,RMS}$) fluctuations in the plume (spatiotemporally averaged across $x=1.25$ to $2.25 \mathrm{cm}$ and over 20–30 μs) from the simulations against the theoretical predictions from Equation (4). This comparison casts light on how closely the rather simple theoretical expressions of the oscillation’s amplitude, estimated based on the ion-wave trapping saturation mechanism, can explain the numerical observations concerning the waves’ strength across a broad range of set-up conditions.

Looking at Figure 12(right column), the theory predicts a constant value for the normalized amplitudes of the number density fluctuations. This constant value overall falls in the middle of numerically observed values of ${\displaystyle \frac{\delta n_{e,RMS}}{n_e}}$ across various set-up conditions. Furthermore, the simulated amplitudes of the electric field fluctuations ($\delta E_{RMS}$) show a substantial discrepancy from the theoretical predictions in some conditions. It is recalled that the $\delta E_{RMS}$ relation in Equation (4) assumes the theoretical characteristics of the IAI’s fastest-growing mode [30]. Therefore, the observed discrepancies between the simulated and the theoretical saturation amplitudes can be partly attributed to the differences in the characteristics of the modes between the simulations and the theory.

As for the specific trends and behaviors with $V$ and $J$ values, the simulated $\delta E_{RMS}$ progressively deviates from the theoretical predictions toward higher voltages. In particular, at $V=600 \mathrm{V}$, there is a substantial difference between the $\delta E_{RMS}$ from the simulations and the theory, with the deviation being more pronounced for xenon and krypton. Across the different values of $J$, there is generally a good agreement between the numerical and the theoretical values of $\delta E_{RMS}$.

In addition to the above results, the DMD analyses performed to isolate the spatiotemporal structures of the dominant instability modes (see Appendix B) revealed a long-wavelength instability with both axial and azimuthal wavenumbers at the low current density limit, i.e., $J=50$ and $100 {\mathrm{Am}}^{-2}$, with the frequencies of about 1 MHz and 0.6 MHz, respectively, for each $J$ value. The spatial structure and the frequency of this instability mode resemble those of the ion transit time instability (ITTI) [62,63,64]. Our observation is also consistent with the findings in ref. [24], where the authors had noted the presence of ITTI at $J$ = 50 ${\mathrm{Am}}^{-2}$ with a frequency on the order of 1 MHz. We see in Section 6.3 that the ITTI has indeed influenced the ions’ distribution function in the low current density conditions.

Another noteworthy observation from the DMD results in Appendix B is that the resolved instability modes feature nonzero axial wavenumbers inside the channel whereas, downstream of the channel exit plane, the instabilities primarily constitute purely azimuthal oscillations. Additionally, as was noticed from Figure 10, the azimuthal wavenumber of the fluctuations ($k_z$) undergoes a change upstream and downstream of the channel exit plane in line with the variation in the local dispersion relation of the instabilities (determined by the local plasma properties).

6.3. Variations in the Ion Species’ Distribution Function

In this section, we look at the velocity distribution functions of the ion species to understand the impacts that the variations in the discharge voltage and the total ion current density have on the ions’ distributions.

In this regard, Figure 13 presents the ion velocity distribution functions (IVDFs) for various set-up conditions. The IVDFs are integrated over the plume region ($1.5 \mathrm{cm}<x<2.5 \mathrm{cm}$) downstream of the acceleration zone. Additionally, the 1D1V ions’ distribution functions along the axial direction under different $J$ conditions are illustrated in Figure 14. In these plots, we also show the ions’ ideal exhaust velocities given by Equation (5) for reference.

From Figure 13, the broadening of the IVDFs toward lower energies can be, in general, attributed to the overlap of the ionization and acceleration zones, and the degree of the broadening is correlated with the extent of this overlap. The spread of the velocities is also evident in the plume from the plots in Figure 14. In particular, we notice that, with increasing $J$, the spread of the velocities decreases and the ion population in the plume becomes increasingly closer to a beam-like distribution. At low $J$ values ($J=50\text{–}100 {\mathrm{Am}}^{-2}$), however, the IVDFs show a notable spread over an extended range of energies and feature an oscillatory distribution as well. These oscillations are visible in the ions’ distribution functions in the plume (Figure 14) and are also noticeable (although not as evidently) in Figure 13 for $J=50\text{–}100 {\mathrm{Am}}^{-2}$. These observations are the signature evidence of the presence of ITTI instability, consistent with what was inferred from the DMD-mode plots for the low $J$ values (Appendix B).

The perturbations in the IVDFs and the 1D1V distribution functions root in the so-called “wave-riding” of the ions [24] because of their interactions with the ITTI waves. In more detail, the electric field oscillations associated with the ITTI create a back-and-forth movement of the acceleration region. When the acceleration region moves upstream, the ions already generated further downstream of the channel experience a reduced portion of the potential drop. Contrarily, when the acceleration zone moves downstream, these ions encounter a larger portion of the potential drop along their path. Accordingly, these oscillations result in the ions being accelerated to velocities below and above the ideal speed.

7. Conclusions

In this article, we reported the results of a parametric study within an axial–azimuthal Hall-thruster-representative configuration and with three propellants of high applied interest—xenon, krypton, and argon. We used our Q2D reduced-order PIC code, IPPL-Q2D, at a high-fidelity approximation order of the 2D problem to investigate how the changes in the values of the discharge voltage and the total ion current density (equivalently, the mass flow rate) affect the macroscopic properties of the plasma and the characteristics of the electrons’ transport, as well as the dominant contributing force terms, the spectra of the instabilities, and the kinetic properties of the ion population. We employed a novel spectral analysis technique—the Optimized Dynamic Mode Decomposition—to enable the characterization of the spatial structure of the resolved fluctuations. We also extensively compared our numerical results across the explored parameter space with the available theories from the literature in terms of the wavenumber, the frequency, and the amplitudes of the azimuthal ion-acoustic instability modes.

In terms of the most salient findings, we first observed that, generally speaking, the numerical dispersion of the instabilities inside the channel and in the plume matches, respectively, the nonlinear and linear theoretical dispersion relations of the IAI. This was overall true over the various studied set-up conditions except for the very low current density values of $J=50 \mathrm{and} 100 {\mathrm{Am}}^{-2}$, where indications of the presence of ITTI were noted. The development of strong ITTI modes at these $J$ conditions can be correlated with the insufficient levels of electrons’ transport which, according to ref. [29], can trigger this instability mode. Evidence of the presence of ITTI in these conditions included the detection of modes with nonzero axial wavenumber components from the DMD analyses, as well as the observation of the “wave-riding” ions in the 1D1V ions’ distribution functions and notable broadening of the axial IVDFs.

Second, ion-wave trapping, which acts as the saturation mechanism for the instabilities, was evidently observed across all set-up conditions except for the lowest current densities ($J\le 100 {\mathrm{Am}}^2$).

Third, regarding the contributions of the individual momentum force terms to the electrons’ axial transport, in most conditions, the viscous effects ($F_{\Pi }$) were the dominant term near the anode boundary. However, in the rest of the domain, the $F_{\Pi }$’s contribution was nearly offset by the convective inertia ($F_I$) term, leaving the force term representing instabilities’ effects ($F_E$) as the primary contributor to the transport.

Fourth, the theoretical predictions of the azimuthal wavenumber and the frequency of the dominant IAI mode had up to 200% and up to 150% error, respectively, compared with the corresponding simulated values across the investigated plasma conditions and propellants. Moreover, the theoretical predictions of the saturated amplitude of the IAI-induced plasma density fluctuations in the plume based on the ion-wave trapping saturation mechanism were about 12% of the steady-state plume density, whereas the simulated values varied between 8 to 15% of the plume density over the different set-up conditions.

Besides the above findings, the main observed variation trends with the discharge voltage ($V$) were the following:

(1)

The width of the acceleration zone progressively expanded with increasing $V$, and the location of the peak of the $E_x$ profile shifted downstream.

(2)

The deviation of the ion velocities in the plume ($V_{ix}$) from the ideal values ($v_{exhaust}$) consistently increased toward higher discharge voltages, with the deviation of the $V_{ix}$ from $v_{exhaust}$ becoming more substantial from xenon to argon.

(3)

The electron current ($I_e$) exhibited small variability across the voltages of 200 to 400 V but showed a jump in values by a factor of about 1.5 from 400 to 600 V, which was traced to the increased role of the viscous force term at these voltages. The ion current density ratio ($I_i/I_0$) appeared to reach a plateau beyond $V=300 \mathrm{V}$, indicating the maximum ion beam current fraction that can be extracted with the imposed injection source intensity.

(4)

The instabilities’ characteristics showed minor variability across the simulated voltage range.

(5)

The $F_{\Pi }$ term became increasingly more significant within a larger segment of the channel region, overtaking the instabilities’ contribution to transport particularly at 600 V.

Furthermore, the most important variations with the total ion current density ($J$) were as follows:

(1)

The overlap between the acceleration and the ionization zones increased at small values of the current density ($J$), which reduced the ions’ axial drift velocity in the plume.

(2)

The fraction of the ion current collected at the cathode ($I_i/I_0$) became increasingly smaller than 1 as the J increased.

(3)

The azimuthal wavenumber, frequency, and amplitude of the instabilities increased with $J$ and so did their induced electrons’ transport.

(4)

The ions’ velocity distribution functions at the low $J$ values of 50 and $100{ \mathrm{Am}}^{-2}$ showed the distinct features of the ions’ interactions with the ITTI.

As the final point, it is acknowledged concerning the observations above that the imposed ionization source in the simulations may have to some extent affected the results. The consideration of an ionization source has been of course consistent with our aims and objectives in this article. It is, nonetheless, important to further verify the outcomes from this work in a controlled experimental setting and/or using 3D PIC simulations that self-consistently resolve the collisional processes (including the ionization) as well as the interactions with the radial walls. The insights derived from the work presented here serve as a crucial foundation and can guide these future efforts.