Hybrid plasma simulations of a magnetically shielded Hall thruster

Magnetic shielding (MS) of a Hall effect thruster (HET) chamber has been proven an effective technique to limit both wall erosion, due to high-energy ion bombardment, and heat loads, thus enabling the design of the next generation of HETs, featuring enhanced performances and operational lifetimes.¹ Due to their recent development, just a few prototypes of MS-HETs have been experimentally tested to date. In the mid to high power range, there are the 4.5 kW BPT-4000,² the 12.5 kW HERMeS,³ the 6 kW H6MS,⁴ the 9 kW H9,⁵ the 20 kW NASA-300MS,⁶ and SITAEL’s 5 kW HT5k⁷ and 20 kW HT20k;⁸ then, in the low-power range (<1 kW) there are the MaSMi-60⁹ and the ISCT-200.¹⁰ Yet, a lower number of studies comparing experimental and simulation results have been realized. Relevant ones were performed with the multi-fluid simulation code Hall2De¹¹ for BPT-4000,¹ H6,^12,13 HERMeS,^14,15 H9,¹⁶ and MaSMi-60.¹⁷ 

The relatively low numbers of MS-HET prototypes and studies, together with the lack of predictive models of HET discharges (due to the open problems on plasma–wall interaction and electron turbulence, particularly), make uncertain the direct extrapolation of results and trends from one thruster system to another. This situation affects, for instance, the design and development of new optimized electric propulsion architectures such as the direct-drive power concept,^18–20 which requires a precise characterization of the thruster performances and the cathode–anode electrical coupling over the nominal operation range for its appropriate integration with other subsystems.²⁰ Therefore, advances on the validation of simulation tools against experimental data, capable not only of providing a full characterization of a HET plasma discharge and performances, but also of addressing thermal, electrical, and material issues related to its operation, are of central interest in HET research.

This work presents 2D numerical simulations of the 5 kW HT5k,⁷ and a partial validation of the numerical model with the limited experimental data existing for this prototype. The simulations are carried out with the code HYPHEN^21,22 and constitute its first test with an MS topology. HYPHEN is a multi-thruster, hybrid-formulation code, which uses a particle-in-cell (PIC) model for heavy species and a magnetized drift-diffusion fluid model for electrons. Contrary to other hybrid codes, relying on a quasi-1D electron model,^23–25 whose application to MS topologies is rather complicated, HYPHEN, as Hall2De, adopts a full-2D electron fluid model on a magnetic field aligned mesh (MFAM), thus allowing a complete characterization of the electron currents, which are crucial to study, for instance, plasma–wall interaction effects. The version of HYPHEN used here is also the first one incorporating a “wall” cathode instead of a “volumetric” cathode. This last one worked fine for HETs with laterally located cathodes (except for the code being axisymmetric, although electron emission becomes quickly homogeneous azimuthally²⁶); however, it was not very suitable for the centrally mounted cathode of the HT5k.

HYPHEN simulates the slow-dynamics, axisymmetric transport of the plasma discharge. This implies that two processes, involving kinetic, non-symmetric, and high-frequency aspects, need to be modeled phenomenologically in the electron fluid equations. The first one is the estimation of the particle and energy fluxes to the thruster dielectric walls. These depend on the electron velocity distribution function (VDF), which is non-Maxwellian because of the partial depletion of the collected VDF’s high-energy tail and the secondary electron emission (SEE) by the walls. Kinetic studies^27–29 and experimental evidence^30,31 are, therefore, used to tune the fluxes to the walls. The second phenomenological model intends to reproduce the slow, turbulent transport of the electron fluid, resultant from averaging (on the azimuth and the high-frequency time scale) the oscillating, azimuthal forces emanating from nonlinear instabilities.^32–36 Since there is not yet an established theory of this averaged turbulence force, authors have treated it as an anomalous collisional force, calibrated, when possible, with experimental data. Still, simulations differ much on the selected strength and shape of the anomalous collision frequency.^{15,23,25,37–40}

The central goal of this work is the analysis of the 2D profiles of the discharge, the identification of the main aspects related to the MS topology and the central cathode, and the effects on thruster performances. The discussion aims at improving the understanding of the plasma physics peculiar to an MS-HET and revealing its potential advantages over a HET with a conventional magnetic topology. The document is organized as follows. Section II describes the HT5k thruster unit and the experimental data supporting this work. Section III describes the simulation model and the main settings. Subsection IV A discusses the tuning of the turbulent parameters to match the experimental data. Then, Subsections IV B and IV C analyze the 2D plasma discharge and the thruster performances. Finally, conclusions are drawn in Sec. V.

SITAEL began to develop the HT5k [Fig. 1(a)] in 2013. This 5 kW-class thruster consists of two main elements: the thruster itself and the HC20 hollow cathode. Along the last few years, SITAEL designed, manufactured, and tested three different development models (DMs) of the HT5k, as part of several national and international programmes. The different prototypes underwent several technical investigations,^7,24 which permitted to demonstrate low erosion, high performance, direct-drive operations, as well as performance stability in high-vacuum conditions (pressure $<10−5$ mbar). The design and manufacturing process of the HT5k thruster unit engineering qualification model started in 2019, when the thruster was chosen as the main propulsive unit for the orbit raising and station keeping of the Ital-GovSatCom geostationary platform.⁴¹ The HT5k DM3 prototype considered here implements a centrally mounted cathode and a non-conventional magnetic field topology, which is beneficial from the point of view of discharge channel erosion. Design and previous testing efforts were dedicated to enhance critical components and to optimize the thruster thermal behavior. The thruster tests were run in SITAEL’s IV10 facility, reaching pressures of the order of $7×10−6$ mbar (Xe) while firing at 4.4 kW of discharge power. The DM3 demonstrated competitive performance and showed stable and efficient operation in the 3–7 kW discharge power range, featuring anodic thrust efficiencies up to $60%$⁠.

Figure 1(b) displays a sketch of the geometry of the thruster chamber and the near plume region. The geometrical parameters $Lc$ and $Hc$ correspond to the thruster chamber length and width, respectively. The plasma domain to be simulated with HYPHEN corresponds to the cylindrical axisymmetric half meridian plane, including the annular thruster chamber and the near plume region. The latter extends axially from the chamber exit plane up to $6Lc$ and radially from the symmetry axis up to $6Hc$⁠. The position of the cathode exit plane, at the thruster symmetry axis, is also indicated in Fig. 1(b). Figure 1(c) shows a scheme of the resistor-inductor-capacitor (RLC) filter unit connecting anode and cathode including the power source voltage $Vs$⁠. The discharge voltage $Vd$ is set between the anode and the cathode. The discharge current $Id$ flowing between anode and cathode is indicated in the sketch according to the flow direction of electrons. The values of the RLC filter elements are R = 4.7 $Ω$⁠, L = 360 $μ$H, and C = 94 $μ$F.

The experimental setup is detailed in Refs. 24 and 42. The experimental data used in the simulations, with xenon as propellant, are listed in Table I. Time-averaged values of the discharge current $Id$ and the thrust $F$ are available for five operation points, hereafter referred to as cases 1–5, defined by a pair (⁠$Vs$⁠, $m˙A$⁠), where $Vs$ is the power source voltage and $m˙A$ is the propellant mass flow injected through the anode to the thruster chamber. The former ranges from 300 to 400 V and the latter from 10 to 14 mg/s. For all cases, a 7.5% of $m˙A$ is injected through the cathode. Data repeatability (i.e., standard deviation) of the measurements is 5%. The value of the time-averaged $Isp$ ranges from 1900s to 2100s from case to case. In addition, the operation point (⁠$Vs$⁠, $m˙A$⁠) = (300 V, 14 mg/s) features an anodic thrust efficiency of 58.2%.

Figure 2(a) shows a schematic representation of the HYPHEN structure and simulation loop, which is briefly outlined next. HYPHEN is an axisymmetric, hybrid, OpenMP-parallelized code built modularly. The code version for HET simulations consists of three main modules: the Ion module (I-module), which follows a Lagrangian approach for simulating the dynamics of the PIC macroparticles of heavy species; the Electron module (E-module), which solves a fluid model for the magnetized electron population; and the Sheath module (S-module), which provides the proper coupling between the quasineutral plasma bulk, and the thruster walls. The E-module assures automatically plasma quasineutrality in the simulation domain. Thus, the Debye sheaths managed by the S-module are, in fact, discontinuity surfaces adjacent to the thruster walls. The three modules are coupled within a time-marching sequential loop.

The I-module operates on a structured mesh of the simulation domain, shown in Fig. 2(b). On the contrary, and in order to limit the numerical diffusion arising from the strong anisotropic transport on magnetized electrons, the E-module uses an unstructured MFAM,⁴³ defined by the externally applied magnetic field $B$ and shown in Fig. 2(c). The magnetic configuration of the MS-HT5k, which features a null magnetic point inside the channel, attempts to screen well all its internal walls.

Two reference frames are considered: one is the cylindrical frame ${1z,1r,1θ}$⁠, with coordinates $(z,r,θ)$⁠; and the second one is the magnetically aligned frame ${1⊥,1∥,1θ}$⁠, with $1∥=B/B$ and $1⊥=1∥×1θ$⁠, and coordinates $(λ,σ,θ)$⁠. The orthogonal magnetic coordinates $λ(z,r)$ and $σ(z,r)$ of the MFAM, Fig. 2(c), are obtained from solving $∇⋅B=0$ and $∇×B=0$⁠.

Let $Zs$⁠, $ns$⁠, $us$⁠, and $js=eZsnsus$ be the charge number, particle density, macroscopic velocity, and current density of the plasma species $s$ (i.e. electrons $e$⁠, neutrals $n$⁠, singly charged ions $i1$⁠, and doubly charged ions $i2$⁠); $E=−∇ϕ$ be the electric field, with $ϕ$ being the electric potential; and $Te$ being the electron temperature. Every simulation step, the I-module takes as inputs $B$⁠, $E$⁠, and $Te$ and performs the following tasks: (i) the propagation of macroparticles one time step $Δt$ forward, according to the electromagnetic fields acting on them; (ii) the injection of new macroparticles into the domain and the removal of exiting ones; (iii) the interaction of macroparticles with the thruster walls, such as neutral reflection and ion recombination; (iv) the generation of new ion macroparticles due to the ionization of neutrals; and (v) the computation, through a particle-to-mesh weighting process, of the macroscopic properties characterizing each heavy species. Further details can be found in Refs. 21, 44, and 45.

The E-module, taking these heavy-species magnitudes as inputs, solves a quasineutral, drift-diffusion fluid model for the magnetized electron population, obtaining $ϕ$⁠, $Te$⁠, $je$⁠, and the electron heat flux vector $qe$⁠. The electron fluid model equations are^37,46,47

Equations (1) and (2), with $ji=∑s≠e,neZsnsus$⁠, correspond to plasma quasineutrality and the plasma current conservation equation, and the right hand sides are inputs from the I-module. The momentum equation (3) neglects electron inertia, assumes the pressure tensor to be isotropic, and includes the resistive force $Fres$ and the turbulent (or anomalous) force $Ft$⁠. The resistive force satisfies

with $me$ being the electron elementary mass, $νe=∑s≠eνes$ being the total momentum transfer frequency due to collisions with all heavy species, $νes$ being the individual contributions for each heavy species $s$⁠, and

is an equivalent heavy species collisional current density.²¹ The turbulent force $Ft$⁠, accounting for azimuth-averaged, wave-based anomalous transport, is modeled phenomenologically as^23,37,38,48

with $νt$ being the turbulent collision frequency, $ωce=eB/me$ being the electron gyrofrequency, and $αt(z,r)$ being the phenomenological function representing the local turbulence level.^23,37 Turbulent-based contributions to the axial and radial momentum equations are negligible compared to the rest of forces there.

The main assumption of the electron drift-diffusion model is to neglect inertia in the momentum equation. This is justified as long as the electron kinetic energy is much less than their thermal energy, i.e., $meue2≪Te$⁠, a condition well satisfied in HET discharges, except in localized regions and certain operation points. The great advantage of the drift-diffusion model is that the momentum equation (3) reduces to a generalized Ohm’s law for $je$ (with an electric conductivity tensor $σe¯¯$⁠), which is much easier to treat numerically than the whole differential equation on $ue$⁠.

The component of Ohm’s law along the local cross-field direction $1⊥$ is⁴⁹ 

where $σe=e2ne/(meνe)$ is the parallel electric conductivity, $χ=ωce/νe$ is the classical Hall parameter, and $χt=χ/(1+αtχ)$ is the reduced Hall parameter when including $νt$⁠. Therefore, the effective Hall parameter is $χχt$ and scales as $∝αt−1/2$ if turbulent transport dominates. The last term on the right side of Eq. (9) represents collisions with heavy species and is negligible generally.

Equation (4) is the electron energy equation for an isotropic pressure tensor in the inertialess limit, where the second term in the left side gathers the enthalpy and heat fluxes, and the right side includes the work of the electric field and the power losses from inelastic collisions (e.g., excitation and ionization). Fourier’s law for the heat flux, Eq. (5), includes the thermal conductivity tensor $Ke¯¯=5Teσe¯¯/(2e2)$⁠, and corresponds to the drift-diffusion limit of the evolution equation for $qe$⁠.⁴⁶ 

The numerical treatment of Eqs. (2)–(5) in the unstructured, irregular MFAM was developed in Refs. 43, 49, and 50. As shown in Fig. 2(c), the MFAM is composed of inner and boundary cells. Inner cells are those enclosed by $B$-parallel (blue) and $B$-perpendicular (red) lines. Boundary cells, however, contain at least one boundary face aligned with a domain boundary, which is not a magnetic line generally. Centroids (or computational points) of both cell and faces correspond to the magnetic center or the geometric center, when the former is not available (e.g., at boundary cells). A finite volume scheme on the MFAM cells is applied to conservation Eqs. (2) and (4); and gradient reconstruction schemes, ad hoc for the unstructured MFAM, are applied to Ohm’s and Fourier’s vector laws (3) and (5). For each variable $ϕ$ and $Te$⁠, this generates one algebraic equation per MFAM cell. Additionally, boundary conditions specifying $jne=1n⋅je$ and $qne=1n⋅qe$ (and discussed below) are imposed at each MFAM boundary face, with $1n$ being the outward unit normal vector. This yields matrix equations for both $ϕ$ and $Te$ at the centroids of the cells and of the boundary faces of the MFAM. A direct solver for sparse linear systems is used for the parallelized computation of the solution.^51,52

The reference $ϕ=0$ for the potential is set at the cathode boundary faces, so the anode potential is $ϕ=Vd$⁠. The RLC filter unit, Fig. 1(c), relates $Vd$ to the imposed source voltage, $Vs$⁠, and the varying discharge current, $Id$⁠, through

This equation is integrated within the E-module over time with a first order numerical scheme.

The simulation domain extends up to the sheath edge of the quasineutral plasma, represented by the MFAM boundary faces. The solution for the (infinitely) thin Debye sheaths is detailed in the  Appendix. There, Eqs. (A2) and (A7) provide the appropriate conditions at each MFAM boundary face in the form of nonlinear relations for $jne$ and $qne$ vs the potential jump across the sheath, $Δϕsh$⁠.

At a dielectric wall, expressions for $jne$ and $qne$ are quite straightforward: a zero-collected electric current yields directly $jne=−jni$⁠; Eq. (A2) is then solved for $Δϕsh$⁠, and Eq. (A7) yields $qne$⁠. At the current-driving anode, with known potential $Vd$⁠, the determination of $jne$ and $qne$⁠, at each anode face of the MFAM, requires to compute previously $Δϕsh$ in the following way: Eq. (A2) is combined with Ohm’s law (3) for $jne$ yielding a non-linear implicit equation for $ϕ=Vd+Δϕsh$ at each anode face. These equations are linearized and introduced in the matrix system yielding all $ϕ$⁠; if needed, iterations on the implicit equations are run.

At the cathode boundary faces, the discharge current $Id$ divided by the cathode area defines the electron current density $jne(>0)$⁠. The electron energy flux at each cathode face, expressed as $qne−(5/2)Tejne/e$⁠, is set equal to $−2Tcjne/e$⁠, with $2Tc$ being the average emission energy per electron. At the (quasineutral) axis, symmetry conditions imply $jne=0$ and $qne=0$⁠; since the I-module yields $∂ne/∂r=0$⁠, one has $∂pe/∂r=0$ and $∂ϕ/∂r=0$ too. Then, at the plume downstream (quasineutral) boundary, a current-free condition $jne=−jni$ and a Maxwellian electron heat flux $qne=−2Tejne/e$ are imposed.^49,50

The time discretization of the electron equations follows a semi-implicit scheme,⁵³ with a sub-time step $Δte=Δt/Ne$ and $Ne=O(1)$⁠. This scheme allows to keep a linear system for $Te$ while reducing the value of $Ne$ required for convergence. Finally, note that mesh interpolation of plasma variables between ion and electron modules is required.

Figures 2(b) and 2(c) show the PIC mesh and the MFAM used to simulate the HT5k. The main characteristics of the meshes and the relevant simulation parameters are listed in Table II. A singular magnetic point (⁠$B=0$⁠) is located inside the thruster chamber. The green line at the left boundary in Fig. 2(b) represents the annular anode wall, while the small black box in Fig. 2(c) indicates the position of the central cathode boundary.

The prescribed xenon mass flow $m˙A$ is injected from a Maxwellian reservoir through the whole annular anode featuring a flat profile with a sonic axial velocity based on its own temperature (Table II). Considering the same injection properties, a neutral mass flow $m˙C=0.075m˙A$ is injected through the cathode boundary together with an electron current equal to $Id$⁠. The emission energy of electrons at the cathode is set to $2Tc=4.5$ eV;⁵⁴ just as a sensitivity check, simulations run with $2Tc=2.25$ eV have shown no observable differences in the discharge, except naturally at the cathode neighborhood.

Wall recombination of ions contributes to the neutral density. Singly and doubly charged ions are generated volumetrically by electron-neutral collisions. Single ionization rates are obtained from the BIAGI database,⁵⁵ while double ionization rates follow the Drawin model,⁵⁶ including the reactions $A+e→A+++3e$ and $A++e→A+++2e$⁠. Neutrals from ion recombination at walls are re-emitted diffusely considering complete ion energy accommodation at the wall, as suggested by several authors.^57,58 Thus, the neutral emission energy is only given by the wall temperature, which is set to 850 K.⁵⁹ Neutrals are reflected diffusely at the wall with zero energy accommodation; Refs. 21 and 45 provide further details on the interaction of heavy-species macroparticles with walls.

The simulations monitor independently the populations of neutrals, singly charged, and doubly charged ions. Each species population is controlled, setting a target number of 500 macroparticles per cell with a $±$10% of tolerance.²¹ 

The (ion) time step in Table II is set so that a typical doubly charged ion takes at least two time steps to cross the smallest PIC cell. The simulations are started by injecting neutrals through the anode and the cathode and considering a minimum background plasma density to trigger the discharge.⁴⁵ Every simulation features a total of 60 000 time steps (equivalent to 900 $μ$s of simulation time) so that $Id$ undergoes a sufficiently large number of low-frequency (i.e., breathing mode) oscillation cycles. Five sub-time steps per ion time step (⁠$Ne=5$⁠) are used to integrate electron equations.²¹ All the results shown in Secs. IV A–IV C are time-averaged over several $Id$ cycles.

In order to complete the electron model, the turbulence function $αt(z,r)$ in Eq. (8) must be chosen. Since for each case in Table I only two experimental parameters, $Id$ and $F$⁠, are known, the function $αt$ will be of the axial “step-out” type shown in Fig. 2(d), with two fitting parameters only, $αt1$ and $αt2(>αt1)$⁠, applying, respectively, and approximately, inside and outside the channel. These step-out profiles have provided good fittings in previous studies.^39,58,61,62 Hereafter, a particular step-out profile is referred to as (⁠$αt1$⁠, $αt2$⁠).

For each of the cases of Table I, the pair (⁠$αt1$⁠, $αt2$⁠) has been tuned to reproduce $(Id,F)$ with a relative error smaller than 5%, thus consistent with the repeatability of the experimental data. The results are in Table III and Fig. 3. Turbulent transport is larger in the plume by nearly one order of magnitude, in line with the existing literature. Also aligned with previous studies for MS-HETs¹³ and traditional HETs,⁶² there is a moderate change of the turbulence parameters with the operation point. The parameter $αt2$ features the largest variation, increasing as $Vs$ decreases and $m˙A$ increases. The obtained turbulence fitting for $(Id,F)$ turns out to be also good reproducing the main oscillations of $Id$⁠, those known as the breathing mode^34,48,63,64 and listed in Table III. For instance, for case 1, simulations yield an oscillation frequency $fd=20.1$ kHz, close to the 22 kHz reported in experiments. For all cases, $ΔId$ ranges between 4% and 13% of $Id$⁠, also consistent with experiments.

It is known that the thrust and, mainly, the discharge current are rather sensitive to the turbulence parameters. This is illustrated here by case 2* in Table III, which corresponds to the operation point of case 2 but using the turbulent fitting of case 1. Simulation errors with respect to experimental data are about 18% for $Id$ and 9% for $F$⁠.

Nonetheless, these differences in performances do not change the main trends of the 2D discharge. Figures 4(a)–4(d) show the 1D axial profiles of $ϕ$⁠, $Ez≡−∂ϕ/∂z$⁠, $Te$⁠, and $ne$⁠, along the thruster channel midline and for cases 1, 2, and 2*. For case 1, most of the ion acceleration occurs in the near plume along 2–3 channel lengths. This agrees with existing experimental measurements for a previous HT5k prototype.⁶⁵ In case 2, with higher $Vs$⁠, the total potential fall spans over a broader axial region, which would explain the slightly higher plume divergence we observe and is in line with previous numerical¹³ and experimental^9,10 studies. The peaks of $Te$ and $Ez$ are outside the channel and move slightly downstream for higher $Vs$⁠, in line with previous numerical studies too.¹³ 

This behavior of $Ez$ follows some experimental results reported in Ref. 10 but differs from those observed in Ref. 66, suggesting that the observed mild trends with $Vs$ depend, at least partially, on the prototype and the operational conditions.

The change of turbulent collisionality from case 2 to 2* in Fig. 4, results mainly in steeper gradients of $ϕ$ and $Te$ in the region where $αt$ is transitioning from $αt1$ to $αt2$⁠, but plasma variables behave very similarly for cases 2 and 2*, and typical MS effects are reproduced in both cases, in line with the numerical studies of Refs. 17 and 67.

Figures 4(e)–4(h) show the effect on the discharge of moving the transition point $(αt1+αt2)/2$ upstream or downstream. For case 1 and at the mid-channel radius, that point was placed at $z/Lc=1.46$ (for reference, the peak of the magnetic field is at $z/Lc=1.29$⁠). Cases $2−$ and $2+$ displace the transition point to $z/Lc=1.28$ and 1.59, respectively. The variations of $Id$ are small, within the accepted range of $5%$⁠, while the value of $F$ remains practically constant. Changes in the 2D discharge inside the chamber are barely perceptible, but mild changes on the location and values of the maximum $Ez$ are observed. This suggests that the location of the transition point can be used to fit the location of the peak of $Ez$⁠.

The strength of the maximum $Ez$ can be adjusted using a third fitting parameter $αt3$ around the peak of the magnetic field, as it is done in “quenched turbulence models,” in general with $αt3<αt1,αt2$⁠.^68–70 One stage ahead in fitting techniques is due to Mikellides and Lopez Ortega,¹⁵ who adjust a multi-piecewise function $αt$ with a larger set of experimental plasma data. Interestingly, the resultant function $αt$ resembles a “step-out plus quenched model,” with $αt3<αt1<αt2$⁠.

In any case, all these are just numerical attempts to adjust time-averaged experimental data, without a firm theoretical basis (not even for extrapolating the function $αt$ from one operation point to another). The fitting problem is more severe for fully 2D MS topologies than for the traditional ones, which feature quasi-radial magnetic lines. Fortunately, this weakness in modeling accurately the slow turbulent transport of electrons does not preclude the study of other central features of the plasma discharge.

Case 1 in Table III is chosen for this analysis. Figure 5 shows the 2D (⁠$z$⁠,$r$⁠) maps of relevant time-averaged plasma properties inside the thruster chamber (left column) and in the whole simulation domain (right column). The neutral density map, Figs. 5(a) and 5(b), illustrates the effects of gas ionization plus the wall-born neutrals from ion recombination. The neutral density exhibits a decrease of about two orders of magnitude inside the thruster chamber, corresponding to a large propellant utilization. In addition, plot (b) shows well the injection and partial ionization of the neutral gas injected through the cathode. Then, the enhanced electron-neutral collisionality favors the coupling of cathode electrons with the main ion beam. The plasma density [Figs. 5(c) and 5(d)] inside the chamber is higher and more uniform than in conventional HETs, due likely to the acceleration region being moved outwards from the chamber. As usually, the plasma density presents its maximum at a central location inside the chamber, here close to the magnetic null point. As shown in Fig. 4, the increase in the $Te$ peak at higher $Vs$ enhances gas ionization, augmenting slightly the peak of $ne$ inside the chamber. In contrast, the lower $ne$ in the near plume at higher $Vs$ is due to the larger ion acceleration. The right plot shows the secondary plasma plume created by ionization of cathode neutrals, which merges downstream with the main plume.

The electric potential [Figs. 5(e) and 5(f)] also presents a maximum close to the singular magnetic point. Potential variations are rather small inside the chamber, as expected in an MS-HET configuration. As commented already, the MS moves the acceleration region outside the chamber and equipotential lines follow approximately the magnetic lines there. Inside, the magnetic null point and the pressure gradients uncouple equipotential lines from magnetic ones. The electron temperature [Figs. 5(g) and 5(h)] peaks in the near plume, rather close to the maximum $Ez$⁠. Temperature isolines follow closely magnetic lines and temperature gradients are very pronounced when the electron flow enters into the chamber. The temperature is below 5 eV near all chamber walls, which is indeed one of the main achievements of MS topologies^12,13,25 and leads to small energy losses to the walls.

Figure 6 depicts 2D vector maps for the longitudinal components of the ion, electron, and electric currents, defined as $ȷ~i=ji−jθi1θ$ and so on. The ion current is obtained from a particle-to-mesh weighting algorithm, while the electron current comes out directly from the fluid model. The ion streamlines, Figs. 6(a) and 6(b), reflect the existence of backward, forward, and lateral ion flows. Although there is a point with $ȷ~i=0$⁠, notice that the ionization source is distributed in the whole channel volume and the streamlines represent the ion macroscopic behavior. Ions, practically unmagnetized, follow the electric field and are not prevented from impacting the channel walls. Interestingly, plot (a) shows nearly wall-parallel ion current streamlines close to the channel exit, a fact also observed experimentally.⁵ Ion streamlines moving into the far plume present the expected divergence outside the thruster chamber. Plot (b) shows again how the secondary ion beam created around the cathode combines with the main beam downstream.

The longitudinal electron streamlines, Figs. 6(c) and 6(d), are all born from the cathode but they are split into one beam directed to the far downstream region, which neutralizes the main ion beam, and one beam moving into the interior of the thruster. This second beam follows first the magnetic lines, but then enters the chamber almost perpendicular to the magnetic field, and finally, it is again channeled by the convergent–divergent magnetic lines. Interestingly, this last configuration does not constitute at all a magnetic nozzle⁷¹ (the magnetic strength is null at the “throat,” instead of maximum) but it channels the electrons anyway. Still, some electron streamlines are driven into lateral walls to cancel the ion flows there. It has also been realized that the neutral injection through the cathode enhances much the diffusive cross-field transport of the inward beam and the electric coupling of the cathode with the main beam. Finally, Figs. 6(e) and 6(f) show the 2D streamlines of the longitudinal electric current. The zero net collected current condition imposed at the dielectric walls and along the downstream boundary (which is assumed current-free locally) yields the expected current loops connecting anode and cathode.

To complete the description of the discharge, Figs. 7(a)–7(c) plot the main plasma magnitudes along the thruster internal walls. The abscissa length $s$ runs from the inner corner at the exit of the chamber toward the anode corner, continues along the anode, and finishes along the outer, interior wall of the chamber. Figure 7(a) shows the electron temperature at the sheath edge and the sheath potential fall, $Δϕsh$⁠, this last one computed from the S-module. The temperature near all internal walls is rather low (⁠$∼$3–4 eV) and the ratio $eΔϕsh/Te$ is between 1 and 3. This ratio is affected by the SEE and the partially depleted VDF of primary electrons, at the dielectric walls, and by the need to control the local flux of electrons, at the anode. The low $Te$ makes the SEE yield small: it is $δs≈$ 0.2 for $Te=5$ eV.

Figure 7(b) plots the electron and ion currents toward the walls (which are constant across the Debye sheaths). The plotted electron current density is the net one, i.e., the difference between the currents of primary and secondary electrons at the wall. Ion and electron current densities are identical at dielectric walls. At the anode, the backward ion current density is around a 20% of the electron current, which is a percentage larger than desirable. Figure 7(c) plots the average ion and electron energy at the wall, per net particle impacting the wall, computed from Eqs. (A5) and (A8). The low values of $Δϕsh$ yield ion-impact energies, $Ei,wall$⁠, well below typical threshold values for wall erosion,^13,72 which is the main advantage expected from MS topologies.

The ion current balance at steady state can be expressed as

where $Iprod$ is the current of ions generated by ionization in the simulation domain; $IiD$⁠, $IiA$⁠, and $IiC$ are the ion currents impacting the dielectric, anode, and cathode walls, respectively [and defined in Fig. 2(b)]; and $Ii∞$ is the ion beam current leaving the domain at plume boundaries, the only one contributing to thrust. All currents are defined as positive; $Iprod$ comes out from a volumetric integration, and the other ones are computed from surface integrals at the domain boundaries.

Table IV details $Iprod$ for cases 1–5 and how it is distributed among the different boundaries; $IiC$ is about one order of magnitude lower than $IiA$ and has not been included. The table also includes the propellant utilization, the current efficiency, and the charge efficiency, defined as

respectively. Here, $m˙=m˙A+m˙C$ and $m˙i∞$ is the total ion mass flow across the plume boundaries. The trends of these partial efficiencies with source voltage and mass flow are the usual ones. Compared to a typical discharge in a conventional HET, $ηu$ and $ηcur$ are rather good, and the relative amount of doubly charged ions, measured by $ηch$⁠, is similar. The relative current losses to the lateral walls, $IiD/Iprod$⁠, are similar to conventional HETs, because of the compensation of two trends: a higher plasma density and a lower electron temperature. For this prototype, the relative current losses to the anode, $IiA/Iprod$⁠, are high, likely due to the existence of the null point not far away from the anode. Although further analyses would be needed, the conclusion here on the current balance is that the MS topology of the HT5k does not present clear advantages over a more conventional one. However, this is not going to be the conclusion with the power balance.

The plasma power balance for the steady state discharge is

where $P=IdVd+PC$ is the total power deposited into the plasma discharge, sum of the discharge power and the net power delivered through cathode electron emission (⁠$PC$ amounting to 1%–2% of $P$⁠); $P∞$ is the plasma energy flow through the plume boundaries; $PD$ and $PA$ are the power losses at the dielectric walls and anode, respectively; and $Pinel$ corresponds to the power losses due to inelastic (ionization and excitation) collisions. All powers are defined as positive. $Pinel$ is obtained from a volumetric integral, $P∞$ comes from a surface integral at the plume boundary, and $PD$ and $PA$ are computed from surface integrals at the respective walls (not at the Debye sheath edges). Equations (11) and (13) have also served to check that the numerical errors given by HYPHEN simulations are acceptable (below 2%).

Inelastic losses correspond entirely to electrons and they are roughly proportional to $Iprod$ (notice that single and double ion creation has different ionization energies). Energy losses to walls D and A present contributions of ions and electrons, which were discussed in Sec. IV B, and they are proportional to their respective currents $IiD$⁠, $IiA$⁠, and $IeA=Id+IiA$⁠. The downstream energy flow $P∞$ corresponds mainly to ions but it also includes the residual electron energy flow in the far plume.

Table V lists the main contributions to the power balance and some related magnitudes for cases 1–5. Observe that while the current losses to the lateral walls amount to about a 40% of $Iprod$⁠, the energy losses to these walls amount to a mere 7%. Adding the energy losses to the anode, the total energy losses to the walls are just 9%–12%, which can be considered an important achievement of MS topologies. Inelastic losses are consistent with the ion production: the ratio $Pinel/Iprod$ yields 23 eV as effective single-ionization cost (including the contribution from excitation collisions).

The thrust efficiency is defined and then factorized as

where $F$ is the thrust measured from plasma properties at the plume boundary. The energy, divergence, and dispersion efficiencies are defined, respectively, as

with $Pz∞$ being the flow of axial plasma (ion and electron) energy across the plume boundaries. Here, $ηene$ quantifies the relative power in the downstream plume, $ηdiv$ assesses the plume divergence based on axial energy and total energy flows, and $ηdisp$ quantifies the level of velocity dispersion of all plasma species (which would be one for a mono-velocity gas). Plume energy flows include the residual energy of electrons, coming from their incomplete expansion in the finite simulation domain, but this is quite low: about 4% of $P∞$ for case 1.

Setting $cos2⁡αdiv=ηdiv$⁠, the half-divergence angles in these simulations are $αdiv∼$ $22°$–$25°$⁠. As commented already, plume divergence increases slightly for higher $Vs$⁠. The corresponding decrease in $ηdiv$ is compensated by $ηene$ and $ηdisp$⁠, thus yielding mild variations of $η$⁠. Previous studies for a different HET prototype also reported slight changes of $η$ in the range $Vs=300−−700$ V.¹³ Here, we find $η∼56%$ (or $∼$61% if the “anodic” thrust efficiency is used), the overestimate relative to measured values being just attributed to the slight overestimate of the simulated $F$⁠.

The numerical simulations of the HT5k prototype with HYPHEN have allowed the testing of this code with MS topologies and a centrally located cathode and its partial validation with the limited experimental data available.

Turbulence transport has been modeled with a step-out profile of a phenomenological anomalous frequency, which has been fitted with experimental data of $Id$ and $F$⁠. The observed trends and the moderate changes of the fitting with the operation point tend to agree with previous works in the literature. Oscillations of the discharge current agree well with experimental ones. A larger set of experimental data would allow a finer tuning of the turbulent transport, but it would probably not alter the present comprehension of the physics of the MS-HET.

A detailed picture of the plasma discharge has been shown for one operation point. The discussion has been focused on those aspects highlighting the effects introduced by the MS topology and the central cathode, over conventional thrusters with lateral cathodes. The MS topology shifts the acceleration region outwards and keeps a high-density, low-temperature plasma inside the chamber. This combination of plasma properties explains that plasma flows to walls are similar to those in conventional HETs. The low plasma temperature in the chamber implies low electric fields, small Debye sheaths around lateral walls, and low SEE, leading to a small energy deposition at the walls and small ion impact energies, well below the usual threshold for material sputtering.

In the near plume, a central cathode with both electron and neutral emission has been simulated. A secondary plasma plume is generated, which merges with the main one downstream. A fraction of the electron emission drifts toward the interior of the thruster perpendicular to the magnetic lines but, once inside, this current is again channeled by the magnetic lines until crossing the null point, except for a small part leaking with ions to the lateral walls.

Finally, detailed experimental measurements of the plasma discharge, in progress, would allow to validate more firmly the present simulation tool, especially those aspects where phenomenological approaches are still important, such as turbulent transport, plasma–wall interaction, and downstream boundary conditions.

This work has been supported by the EDDA project, funded by the European Union’s Horizon 2020 Research and Innovation Program, under Grant Agreement No. 870470.

The authors have no conflicts of interest to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The S-module of HYPHEN solves a planar, unmagnetized, collisionless, kinetic model of the thin Debye sheaths develop around the walls. The model includes SEE and retains other non-Maxwellian features of the electron VDF.⁷³ The S-module provides the appropriate boundary conditions for the quasineutral electron fluid equations, specifying, at each quasineutral MFAM boundary face (i.e., Debye sheath edge), the perpendicular electron current and heat flux.

Taking as reference a ceramic material with large SEE, the electric current density of primary (p) electrons from the quasineutral plasma into the wall is assumed to follow

which corresponds to a partially depleted, partially reflected Maxwellian VDF.⁷³ Here, $Δϕsh$ is the potential fall in the sheath, $c¯e=8Te/(πme)$⁠, $δr$ is the fraction of primary electrons reaching the wall but being reflected back elastically by it, and $σrp$ estimates the replenishment fraction of the VDF tail corresponding to impacting electrons. If $δs$ is the SEE yield, the net, local electron current density to the wall is

Kinetic studies of plasma–wall interaction^27,28 show that $σrp<1$ because of the weak electron collisionality. Simulations results with $σrp$ between 0.05 and 0.30 have revealed small effects on thruster performances and relevant plasma profiles. This is in agreement with experimental results evidencing a reduced influence of plasma–wall interaction on the plasma discharge in magnetically shielded thrusters.^30,31 Based on this, $σrp=0.1$ is set for all simulations here, leaving further sensitivity analyses or more accurate models to further work.

For a conducting wall, such as the anode, we just take $δs,δr≈0$⁠. For the lateral ceramic walls, $δr$ and $δs$ are modeled according to Refs. 73–77,

with $δr0$⁠, $Er$⁠, and $E1$ being material dependent parameters and $δs∗$ being the effective upper-bounded SEE yield, corresponding to a space-charge limited (SCL) sheath. For the boron nitride walls of our thruster, the following values were taken: $δr0=0.4$⁠, $Er=40$ eV, $E1=50$ eV, and $δs∗=0.986$⁠.

The local net power density deposited by the whole electron population at the wall is

where $Ee,wall$ is the average energy per net collected electron and $2Ts$ is the average energy per wall-emitted electron (equal to 4 eV in the simulations here). Then, the net power density of electrons at the sheath edge is

At that edge, the kinetic electron model inside the sheath is matched with the fluid electron model of the outer quasineutral domain. The matching of $jne$ is obvious, but for the energy fluxes, the “kinetic” one, $Pne,edge″$ of Eq. (A6), is matched to the “fluid” one, $−5Tejne/(2e)+qne$⁠. This yields the heat flux at the sheath edge as

which is used as boundary condition on the quasineutral domain for the electron fluid.

For ions, the PIC solution in the quasineutral domain is also matched to the kinetic sheath model at the sheath edge. Dedicated particle-to-surface weighting schemes for the PIC formulation of ions^21,78,79 yield directly the net ion power density at (the PIC side of) the sheath edge, $Pni,edge″$⁠. Then, the power density deposited at the wall is $Pni,wall″=Pni,edge″+jniΔϕsh$ and the average energy per wall-impacting ion is