Coupled plasma transport and electromagnetic wave simulation of an ECR thruster

For the heavy species, a PIC formulation with Monte Carlo collisions is used [41]. Macroparticles evolve in a structured mesh subject to (i) a particle mover that propagates the trajectory of the macroparticles with a leap-flog scheme subject to the quasi-steady fields −∇ϕ and B ₀; (ii) various types of collisions, including ionization (e.g. single, double, and single to double) by electron bombardment and charge-exchange collisions among neutrals and ions [49], (iii) the interaction with the different channel walls and their Debye sheaths, which includes neutral accommodation at the walls and re-emission, ion recombination and re-emission as neutrals, and the fulfilment of the Bohm criterion at the sheath edge [50]. At the position of the injector, a prescribed mass flow rate $\dot {m}$ of neutrals is injected into the domain. Heavy particles reaching the downstream open boundary are removed. Neutrals, singly-charged, and doubly-charged ions are treated as different heavy species and the corresponding macroparticles are kept in different computational lists, which are used for weighting their macroscopic properties onto the PIC mesh. Figure 3(a) plots the structured cylindrical-type mesh used in the I-module. The mesh is finer close to the walls to improve the characterization of gradients in the vicinity of walls. An statistical population control mechanism is implemented to each species looking for an optimal number and size of macroparticles in each cell [51]. A typical number of macroparticles per cell for this type of simulations is around 300. In order to improve the population control in the plume, the mesh cell-size is increased downstream.

Zoom In Zoom Out Reset image size

In their low-frequency response, electrons are treated as a magnetized diffusive fluid [47]. The set of fluid equations is the following:

$$\begin{equation}\nabla \cdot {\boldsymbol{j}}_{\mathbf{e}}=-\nabla \cdot {\boldsymbol{j}}_{\mathbf{i}},\end{equation} \tag{ 1 }$$

$$\begin{equation}0=-\nabla {p}_{\mathrm{e}}+e{n}_{\mathrm{e}}\nabla \phi +{\boldsymbol{j}}_{\mathrm{e}}{\times}{\boldsymbol{B}}_{0}+{\boldsymbol{F}}_{\text{coll}}+{\boldsymbol{F}}_{\text{turb}},\end{equation} \tag{ 2 }$$

$$\begin{align}& \frac{\partial }{\partial t}\left(\frac{3}{2}{p}_{\mathrm{e}}\right)+\nabla \cdot \left(\frac{5}{2}{p}_{\mathrm{e}}{\boldsymbol{u}}_{\mathrm{e}}+{\boldsymbol{q}}_{\mathrm{e}}\right)={\boldsymbol{u}}_{\mathrm{e}}\cdot \nabla {p}_{\mathrm{e}}\\ & \qquad -{\boldsymbol{u}}_{\mathrm{e}}\cdot \;\left({\boldsymbol{F}}_{\text{coll}}+{\boldsymbol{F}}_{\text{turb}}\right)+{Q}_{\mathrm{a}}+{Q}_{\text{coll}},\end{align} \tag{ 3 }$$

$$\begin{align}& 0=-\frac{5}{2e}\nabla \left({p}_{\mathrm{e}}{T}_{\mathrm{e}}\right)-\left({\boldsymbol{q}}_{\mathrm{e}}+\frac{5}{2}{p}_{\mathrm{e}}{\boldsymbol{u}}_{\mathrm{e}}\right){\boldsymbol{B}}_{0}\\ & \qquad +\frac{5{p}_{\mathrm{e}}}{2}\nabla \phi -\frac{{m}_{\mathrm{e}}{\nu }_{\mathrm{e}}}{e}{\boldsymbol{q}}_{\mathrm{e}}+{\boldsymbol{Y}}_{\text{turb}}.\end{align} \tag{ 4 }$$

Equation (1) states the conservation of electric current, with j _i = ∑_s≠e Z_s en_s u _s the total positive current provided by the I-module. Equation (2) is the momentum equation, where electron inertia has been neglected; p_e = n_e T_e is the (isotropic) electron pressure; and F _coll = −m_e n_e∑_s≠e ν_es ( u _e − u _s ) is the resistive force due to collisions, with ν_es = ν_es (T_e) the effective collision frequency of electrons at temperature T_e with species s, which includes ionization and excitation collisions and elastic electron collisions with neutrals and ions.

Finally, the term F _turb = 1_θ α_t j_θe B₀ in equation (2) is a phenomenological turbulence model, with α_t an empirical parameter adjusted to match experimental results. Turbulence is believed to be the main source of anomalous transport in Hall thrusters, and this model (which results in a net enhancement of the electron transport perpendicular to the magnetic field) has been used successfully in the past to explain experimental measurements [52, 53]. Evidence of turbulence on electrodeless plasma thrusters like the ECRT is still very limited, but existing results [54] indicate that it is a pervasive element of these devices (and HPTs) as well. As the ECRT also features an E × B plasma like Hall thrusters, the same model is used here to introduce the effect of turbulent transport in the simulation. After a sensitivity analysis on α_t, which showed it affects substantially the peak electron temperature, but little the propulsive performance figures (such as thrust and specific impulse), a value α_t = 0.02 was found to match well the experimentally available T_e measurements at the axis downstream [55].

Equation (3) is the internal energy equation, where: q _e is the electron heat flux; Q_coll = −n_e ∑_s≠e ν_es ɛ_es , (with ɛ_es the energy yield of the collision) is the power sink due to collisions with heavy species, mainly ionization and excitation of neutrals. Equation (4) is a diffusive model for the electron heat flux, structurally similar to the momentum equation, with ν_e = ∑_s≠e ν_es , and Y _turb = −1_θ α_t q_θe B₀ a turbulence-based contribution, here enhancing the perpendicular heat flux. Additionally, the quasineutrality condition

$$\begin{equation}{n}_{\mathrm{e}}=\sum\limits _{s\ne \mathrm{e}}{Z}_{s}{n}_{s}\end{equation} \tag{ 5 }$$

is imposed in the bulk of the simulation domain, i.e., except in the sheaths. In this expression, the right-hand side is the positive charge density, provided by the I-module, with Z_s the charge number of each species.

Boundary conditions for the electron fluid equations depend on the type of surfaces. On the axis, symmetry conditions are imposed, implying that the radial electric current and electron heat flux are zero. On the current-free downstream boundaries, the net normal electric current, j_n = j ⋅ 1_n, is set to zero and the normal electron heat flux satisfies q_ne = q _e ⋅ 1_n = 2T_e n_e u_ne. At the thruster walls the normal electric current density j_n and the normal heat flux q_ne are computed by the sheath solver described next.

The Sheath solver relates the electron magnitudes at the sheath edge and the thruster wall. In the thruster simulated here, walls are covered by boron nitride, which has a high SEE yield, defined as the ratio of secondary-to-primary electron fluxes. The macroscopic SEE yield, δ_w(T_e) is modelled according to reference [56], which is based on the following parameters:

$$\begin{equation}{\delta }_{\mathrm{w}}={\delta }_{\mathrm{w}\mathrm{s}}+{\delta }_{\mathrm{w}\mathrm{r}},\qquad {\delta }_{\mathrm{w}\mathrm{s}}=\frac{2{T}_{\mathrm{e}}}{{E}_{\mathrm{s}}},\qquad {\delta }_{\mathrm{w}\mathrm{r}}=\frac{{\delta }_{r0}{E}_{\mathrm{r}}^{2}}{{\left({T}_{\mathrm{e}}+{E}_{\mathrm{r}}\right)}^{2}},\end{equation} \tag{ 6 }$$

where δ_ws(T_e) and δ_wr(T_e) are yields for, respectively, (true) secondary electrons (emitted with a small temperature T_se), and elastically reflected primary electrons. For high electron temperatures δ_ws is limited to a maximum of 0.983, when charge saturation is expected to happen [57]. The energies E_s and E_r are δ_r0 are material dependent. Additionally, the sheath model takes into account the partial depletion of the high-energy tail of primary electrons lost into the walls, with the replenishment parameter σ (equal to 100% for full replenishment). The net normal electric current to the dielectric wall, j_n, with contributions of ions and primary and secondary electrons, is zero. As a result the sheath local potential fall between sheath edge Q and (dielectric) wall W, fulfills the relation

$$\begin{equation}e\frac{{\phi }_{Q}-{\phi }_{W}}{{T}_{\mathrm{e}}}=\mathrm{ln}\left[{\left.\frac{e{n}_{\mathrm{e}}{\bar{c}}_{\mathrm{e}}}{4{j}_{\mathrm{n}\mathrm{i}}}\right\vert }_{Q}\left(1-{\delta }_{\mathrm{w}\mathrm{s}}\right)\left(1-{\delta }_{\mathrm{w}\mathrm{r}}\right)\sigma \right].\end{equation} \tag{ 7 }$$

The normal heat flux to the wall, q_ne, is the difference of the normal flux of total electron energy, adding the contributions of primary and secondary electrons, minus the normal flux of convective energy of the electron fluid [56]. Settings here are for boron nitride: E_s = 50 eV, E_r = 40 eV, δ_r0 = 0.4, T_se = 2 eV and σ = 0.3 [52, 56, 58].

In order to avoid numerical diffusion caused by the large magnetic anisotropy, the electron fluid equations are solved in the magnetically-aligned mesh plotted in figure 3(b). Specific numerical algorithms have been developed for this [47]. For the conservation laws (1) and (3), these are based on finite volume schemes. For equations (2) and (4), which are actually state equations relating j _e and q _e with ϕ and T_e gradients respectively, the algorithms are based on finite difference schemes. Irregular cells, next to boundaries require special attention for their accurate solution.

The W-wave module is an axisymmetric, frequency-domain full-wave model. This model takes as inputs the maps of electron density n_e and collisionality ν_e from the E-module. As outputs, it returns the high-frequency EM fields $\boldsymbol{E}\left(z,r,\theta ,t\right)=\mathfrak{R}\left[\hat{\boldsymbol{E}}\left(z,r,\theta \right)\mathrm{exp}\left(-\mathrm{i}\omega t\right)\right]$, $\boldsymbol{B}\left(z,r,\theta ,t\right)=\mathfrak{R}\left[\hat{\boldsymbol{B}}\left(z,r,\theta \right)\mathrm{exp}\left(-\mathrm{i}\omega t\right)\right]$, where variables with hat are the complex field amplitudes. From these fields, the map of power absorption density, Q_a, is then computed.

The frequency-domain of Faraday and Ampère–Maxwell's laws determine the complex amplitudes of the fields:

$$\begin{equation}\nabla {\times}\hat{\boldsymbol{E}}=\mathrm{i}\omega \hat{\boldsymbol{B}},\end{equation} \tag{ 8 }$$

$$\begin{equation}\nabla {\times}\hat{\boldsymbol{B}}=-\mathrm{i}\omega {\mu }_{0}\hat{\boldsymbol{D}}+{\mu }_{0}{\hat{\boldsymbol{j}}}_{\mathrm{a}},\end{equation} \tag{ 9 }$$

where μ₀ is the magnetic permeability of free space, and ${\hat{\boldsymbol{{\jmath}}}}_{\mathrm{a}}$ and $\hat{\boldsymbol{D}}$ represent the complex amplitudes of the excitation current and the electric displacement field, respectively. Note that $\hat{\boldsymbol{B}}$ and $\hat{\boldsymbol{D}}$ automatically satisfy $\nabla \cdot \hat{\boldsymbol{B}}=0$ and $\nabla \cdot \hat{\boldsymbol{D}}=-\mathrm{i}\nabla \cdot {\hat{\boldsymbol{{\jmath}}}}_{\mathrm{a}}/\omega$.

The field $\hat{\boldsymbol{D}}$ includes the contribution of vacuum permittivity and the susceptibility of the plasma. A cold plasma model is used to relate $\hat{\boldsymbol{D}}$ linearly to $\hat{\boldsymbol{E}}$,

$$\begin{equation}\hat{\boldsymbol{D}}={\varepsilon }_{0}\bar{\bar{\kappa }}\cdot \hat{\boldsymbol{E}},\end{equation} \tag{ 10 }$$

where ɛ₀ is the permittivity of free space and $\bar{\bar{\kappa }}$ is the dielectric tensor, given in the {1_z , 1_r , 1_θ } vector basis by

$$\begin{equation}\bar{\bar{\kappa }}=\left[\begin{matrix}\hfill \mathrm{c}\beta \hfill & \hfill -\mathrm{s}\beta \hfill & \hfill 0\hfill \\ \hfill \mathrm{s}\beta \hfill & \hfill \mathrm{c}\beta \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right]\cdot \left[\begin{matrix}\hfill \mathcal{P}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \mathcal{S}\hfill & \hfill -\mathrm{i}\mathcal{D}\hfill \\ \hfill 0\hfill & \hfill \mathrm{i}\mathcal{D}\hfill & \hfill \mathcal{S}\hfill \end{matrix}\right]\cdot \left[\begin{matrix}\hfill \mathrm{c}\beta \hfill & \hfill \mathrm{s}\beta \hfill & \hfill 0\hfill \\ \hfill -\mathrm{s}\beta \hfill & \hfill \mathrm{c}\beta \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right],\end{equation} \tag{ 11 }$$

where cβ and sβ represent the cosine and sine of the angle β formed by the applied magnetic field B ₀ with the z axis, and [59]

$$\begin{equation}\mathcal{S}=1-\sum\limits _{s}\frac{{\omega }_{ps}^{2}\left(\omega +\mathrm{i}{\nu }_{s}\right)}{\omega \left[{\left(\omega +\mathrm{i}{\nu }_{s}\right)}^{2}-{\omega }_{cs}^{2}\right]},\end{equation} \tag{ 12 }$$

$$\begin{equation}\mathcal{D}=\sum\limits _{s}\frac{{\omega }_{cs}{\omega }_{ps}^{2}}{\omega \left[{\left(\omega +\mathrm{i}{\nu }_{s}\right)}^{2}-{\omega }_{cs}^{2}\right]},\end{equation} \tag{ 13 }$$

$$\begin{equation}\mathcal{P}=1-\sum\limits _{s}\frac{{\omega }_{ps}^{2}}{\omega \left(\omega +\mathrm{i}{\nu }_{s}\right)},\end{equation} \tag{ 14 }$$

where:

$$\begin{equation}{\omega }_{ps}^{2}=\frac{{e}^{2}{Z}_{s}^{2}{n}_{s}}{{{\epsilon}}_{0}{m}_{s}},\qquad {\omega }_{cs}=\frac{{Z}_{s}e{B}_{0}}{{m}_{s}},\end{equation} \tag{ 15 }$$

are, respectively, the plasma frequency and (signed) cyclotron frequency of species s, and ν_s its damping frequency. At the frequency range of the device (a few GHz), the electron response dominates, and the other species (singly- and doubly-charged ions) can be safely neglected. The influence of ν_e on EM field propagation and absorption in the context of ECRTs was analyzed in [60]. There, it was shown that the actual value of the damping constant does not modify substantially the integrated power absorption at the resonance, and its only major effect was to thicken the absorption region. In the present model, ν_e is the effective electron collision frequency (i.e. the sum of the collision frequencies with the heavy species) from the E-module.

Combining equations (8) and (9) to eliminate $\hat{\boldsymbol{B}}$, one finds the inhomogeneous Maxwell's wave equation also referred to as the double-curl wave equation for the electric field:

$$\begin{equation}\nabla {\times}\left(\nabla {\times}\hat{\boldsymbol{E}}\right)-{k}_{0}^{2}\bar{\bar{\kappa }}\cdot \hat{\boldsymbol{E}}=\mathrm{i}\omega {\mu }_{0}{\hat{\boldsymbol{{\jmath}}}}_{\mathrm{a}},\end{equation} \tag{ 16 }$$

where ${k}_{0}=\omega \sqrt{{\varepsilon }_{0}{\mu }_{0}}$ is the wave number in vacuum.

The weak form of equation (16) is [61]:

$$\begin{equation}\begin{aligned}\hfill & {\int }_{{\Omega}}\left[\left(\nabla {\times}{\hat{\boldsymbol{W}}}^{{\ast}}\right)\cdot \left(\nabla {\times}\hat{\boldsymbol{E}}\right)-{k}_{0}^{2}{\hat{\boldsymbol{W}}}^{{\ast}}\cdot \bar{\bar{\kappa }}\cdot \hat{\boldsymbol{E}}\right]\mathrm{d}V\hfill \\ \hfill & \qquad +{\int }_{\partial {\Omega}}{\hat{\boldsymbol{W}}}^{{\ast}}\cdot \left[{\boldsymbol{1}}_{\mathrm{n}}{\times}\left(\nabla {\times}\hat{\boldsymbol{E}}\right)\right]\mathrm{d}S=\mathrm{i}{\mu }_{0}\omega {\int }_{{\Omega}}{\hat{\boldsymbol{W}}}^{{\ast}}\cdot \hat{{\boldsymbol{j}}_{\mathrm{a}}}\enspace \mathrm{d}V,\hfill \end{aligned}\end{equation} \tag{ 17 }$$

where Ω, ∂Ω are the simulation domain and its enclosing boundary, $\hat{\boldsymbol{W}}$ is a trial (i.e. weight) function, and ∗ denotes the complex conjugate transpose. In this work, we decompose the fields in azimuthal modes of mode number m, such that $\hat{\boldsymbol{E}}\left(z,r,\theta \right)=\tilde {\boldsymbol{E}}\left(z,r\right)\mathrm{exp}\left(\mathrm{i}m\theta \right)$. As the ECRT exhibits axisymmetry, the only azimuthal mode of interest for this study is the purely axisymmetric mode m = 0. The field $\tilde {\boldsymbol{E}}$ is expanded using a mixed basis of shape functions:

$$\begin{equation}\tilde {\boldsymbol{E}}=\sum\limits _{i}{a}_{i}{\boldsymbol{N}}_{i}\left(z,r\right)+\sum\limits _{j}{b}_{j}{P}_{j}\left(z,r\right){\mathbf{1}}_{\theta },\end{equation} \tag{ 18 }$$

where N _i are Nédélec (edge) functions and P_j are Lagrange (nodal) functions, i, j are indices that run for all edges and nodes, respectively, and a_i , b_j are the complex expansion coefficients for the meridian and azimuthal electric field. This decomposition is known to avoid spurious solutions present in other FE formulations of Maxwell equations [60–62].

The second term of equation (17) must be evaluated at the boundary of the domain. On the thruster walls the effect of the thin layer of boron nitride on the high-frequency fields is ignored, and the metal behind is treated as a perfect electric conductor,

$$\begin{equation}{\boldsymbol{1}}_{\mathrm{n}}{\times}\tilde {\boldsymbol{E}}=\mathbf{0},\end{equation} \tag{ 19 }$$

resulting in homogeneous Dirichlet conditions on both N_j and P_j . Condition (19) is also applied to the core and shield of the simulated segment of coaxial cable. Natural boundary conditions are applied on the downstream and lateral boundaries of the plasma plume:

$$\begin{equation}{\boldsymbol{1}}_{\mathrm{n}}{\times}\left(\nabla {\times}\tilde {\boldsymbol{E}}\right)=\mathbf{0}.\end{equation} \tag{ 20 }$$

At the axis, axisymmetric boundary conditions [62, 63] are prescribed:

$$\begin{equation}{\tilde {E}}_{\mathrm{r}}={\tilde {E}}_{\theta }={\left(\nabla {\times}\tilde {\boldsymbol{E}}\right)}_{\mathrm{r}}={\left(\nabla {\times}\tilde {\boldsymbol{E}}\right)}_{\theta }=0.\end{equation} \tag{ 21 }$$

Lastly, to model the input of high-frequency EM power into the domain, a lumped port condition [64] is set upstream of the coaxial cable (at z = −2λ₀), between the core and the shield conductors, which are electrically connected to the inner rod element and the lateral thruster wall, respectively. As a result, the transverse electromagnetic mode (TEM) is excited in the coaxial cable that delivers power to the plasma inside the thruster, and allows for any reflected power to return unimpeded through the coaxial.

Substituting equation (18) into (17), and using Galerkin's approach so that the same functional space is chosen for $\tilde {\boldsymbol{E}}$ than for $\tilde {\boldsymbol{W}}$, we construct a sparse square linear problem for the coefficients a_i and b_j . The modular finite element methods (MFEM) library [65] is used to construct the matrix and obtain the solution.

The absorbed power density Q_a can be computed as the real part of the scalar product between the high-frequency plasma current and the high-frequency electric field, which after some manipulation yields:

$$\begin{align}\hfill {Q}_{\mathrm{a}}& =-\frac{\omega {{\epsilon}}_{0}}{2}\mathfrak{I}\left({\tilde {\boldsymbol{E}}}^{{\ast}}\cdot {\bar{\bar{\kappa }}}^{{\ast}}\cdot \tilde {\boldsymbol{E}}\right)=\frac{\omega {{\epsilon}}_{0}}{2}\frac{{\tilde {\boldsymbol{E}}}^{{\ast}}\cdot {\bar{\bar{\kappa }}}^{\mathrm{A}}\cdot \tilde {\boldsymbol{E}}}{i}\hfill \\ \hfill & =\frac{\omega {{\epsilon}}_{0}}{2}\left[\mathfrak{I}\left(\mathcal{S}\right)\left(\vert {\tilde {E}}_{\perp }{\vert }^{2}+\vert {\tilde {E}}_{\theta }{\vert }^{2}\right)+\mathfrak{I}\left(\mathcal{P}\right)\vert {\tilde {E}}_{{\Vert}}{\vert }^{2}\right]\hfill \\ \hfill & \quad +\omega {{\epsilon}}_{0}\enspace \mathfrak{I}\left(\mathcal{D}\right)\mathfrak{I}\left({\tilde {E}}_{\perp }^{{\ast}}{\tilde {E}}_{\theta }\right),\hfill \end{align} \tag{ 22 }$$

where ${\bar{\bar{\kappa }}}^{\mathrm{A}}=\left(\bar{\bar{\kappa }}-{\bar{\bar{\kappa }}}^{{\ast}}\right)/2$ is the antihermitian part of $\bar{\bar{\kappa }}$ associated to the presence of ν_e in equations (12)–(14) and is the only contribution of the dielectric tensor to include power absorption by the plasma.

The W-module uses an unstructured triangular mesh that allows reproducing complex geometries and also to locally refine the mesh as needed. A sketch of the W-module mesh is shown in figure 3(c). The open-source software Gmsh [66] is used for the mesh generation. As shown in sections below, the high-frequency field solution features substantial differences among regions of the domain due to the existence of qualitatively distinct propagation regimes, which depend on the electron density n_e and applied magnetic field B ₀. The transition from one region to another is delimited by a cut-off and/or resonance surface. In each propagation regime, zero, one or two EM modes may exist, and the characteristic wavelength may vary substantially for each mode and with the propagation direction. In order to keep the number of mesh elements per wavelength N_λ constant across the simulation domain, the mesh size is adapted for each region. A local minimum wavelength is estimated at each point of the domain from the dispersion relation of each principal propagating mode. The most restrictive regions in terms of minimum wavelength are seen to be the vicinity of the ECR and then the upper-hybrid resonance (UHR) region [59]. The mesh size is then computed according to this minimum wavelength, so that N_λ > 20 across the simulation domain. This scheme is seen to stabilize and converge the fields close to these resonance surfaces, which at low mesh resolution can display numerical noise. A finer mesh size is also used near domain boundary corners and the coaxial to thruster chamber transition. The refinement and numerical verification of the W-module, discussed in appendices A and B respectively, characterize the error committed with the chosen formulation and mesh approach. A rougher mesh is used during the first iterations of the code to speed up convergence. As can be noticed, the number of elements of the W-module mesh is much greater than the ones used for the transport modules, in at least one order of magnitude.

The ECR thruster used for this simulation is that of the MINOTOR prototype being matured at ONERA, as shown in figure 1. The applied magnetic field generated by an annular permanent magnet is shown in figure 3(b), whose local unit vector points downstream, defining magnetic streamlines. The operational point is the nominal one for the device [11, 16]. The geometric, operational, and numerical parameters are listed in table 1. The ECR surface is located approximately 3 mm downstream the backplate.

Table 1. Geometrical, operational and simulation parameters.

Parameter	Name	Units	Value
lr	Inner rod length	cm	2
rr	Inner rod radius	cm	0.115
L	Lateral wall length	cm	1.51
R	Lateral wall radius	cm	1.375
Lp	Plume domain length	cm	4
Rp	Plume domain radius	cm	2.75
zinj	Injection surface center z	cm	0.0
rinj	Injection surface center r	cm	0.5735
tinj	Injection surface width	cm	0.229
rc	Coaxial shield radius	cm	0.3
f	Wave frequency	GHz	2.45
$\dot {m}$	Neutral mass flow rate	mg s−1	0.2
Pf	Input power	W	30
ΔtI	I-module timestep	ns	10
ΔtE	E-module timestep	ns	0.25
ΔtW	W-module update timestep	μs	50
uinj	Propellant injection velocity	m s−1	300
Tinj	Propellant injection temperature	eV	0.02
αt	Anomalous transport coefficient	—	0.02
Nel,I	PIC mesh's cell number	—	3360
Nel,E	MFAM's cell number	—	2550
Nel,W	Wave mesh's cell number	—	120 782

The instantaneous plasma properties computed in the coarser PIC mesh and the MFAM are smoothed using a Gaussian filter when interpolated to the EM mesh, and conversely, spatial anti-aliasing is applied when transporting Q_a from the W-mesh to the MFAM.

The I-module has a timestep of Δt_I, which ensures that the fastest ion particles do not cross more than half a cell of the PIC mesh per timestep. The E-module is run 40 times per I-module step. Since the electron density n_e and effective collisionality ν_e vary only slowly as the steady-state is approached, the W-module is only run every 5000 I-module steps.

To initialize the simulation, the I-module is run alone for a number of steps to fill the domain with particles, using a simplified, Boltzmann-relation model to emulate the electron response. Then, the E-module is activated with a uniform Q_a map for an additional number of initialization steps. Finally, the coupled I-, E-, and W-modules are run together until steady state conditions are reached. Convergence of the plasma variables to steady state conditions is achieved already after 1.5 ms of simulation time. The simulation results of next section are shown at the end of the simulation, at 3.5 ms. Figure 4 illustrates the initialization and steady-state convergence of the simulation by plotting the computed thrust contributions of the different species. The updates of the power absorption profile are indicated by vertical dashed lines.

Zoom In Zoom Out Reset image size

In the following sections, all plasma variables shown have been averaged over 500 ion module time steps after reaching the converged state to reduce numerical noise.

In the steady-state, the neutral propellant density n_n peaks at 6.5 × 10¹⁹ m⁻³ near the injection port (see figure 5(a)). Then, it decreases away from this location featuring a nearly homogeneous profile inside the thruster chamber, adding the contribution of ion recombination at the walls. Neutral density decreases as neutrals expand downstream and ionization takes place, falling more than two orders of magnitude downstream. This decrease is most evident close to the injector, where the ionization is dominant. Neutral density is lowest in the peripheral plume, essentially out of reach for the ballistic neutrals. The fraction of propellant mass flow rate leaving the discharge chamber as ions is roughly 50%.

Zoom In Zoom Out Reset image size

The singly- and doubly-charged ion densities, n_i1 and n_i2, are shown in figures 5(b) and (d). As it can be noted, n_i1 is about one order of magnitude greater than n_i2, and both follow the same overall behavior. The maximum singly- and doubly-charged ion densities are of the order of 9.5 × 10¹⁷ m⁻³ and 9.6 × 10¹⁶ m⁻³, respectively. The peak densities are reached in the vicinity of the injector, where most of the ionization takes place. Ion densities decrease toward the thruster walls, and ions impacting there are recombined into neutrals. The ion densities drop noticeably around the inner rod element and in the magnetic tube that originates in this region, which coincides with the region of higher electron temperature as discussed further below. This feature of the plasma discharge was observed to be robust to a sensitivity analysis of the simulation. Consistently with n_i2 ≪ n_i1, the electron density n_e, shown in figure 5(c), is essentially equal to n_i1, peaking around 1.1 × 10¹⁸ m⁻³. The electron density decreases downstream along magnetic lines as the plasma expands, reaching values of 2.7 × 10¹⁶ m⁻³ at the symmetry axis at the end of the simulated MN domain. This value is consistent with estimations performed with other models as the quasi 1D model presented in reference [55].

The quasi-steady map of electron temperature T_e is shown in figure 5(e), with dashed lines indicating magnetic lines for reference. Due to the strong electron magnetization, effective energy transport rates of electrons along and across magnetic field lines are markedly different. While there exists large temperature gradients perpendicular to the applied magnetic field B ₀, a near-isothermal behavior is observed in the parallel direction, at least in the finite domain simulated here. Gradual parallel cooling of electrons is expected further downstream, an aspect of the electron expansion that requires a kinetic treatment to be modelled consistently [67].

The maximum electron temperature is 43 eV and is found near the dielectric window that connects the thruster chamber with the coaxial cable, and in the magnetic tube that emanates from this location. As discussed later, the neighborhood of this window is also the location of the maximum power absorption. The electron temperature close to the symmetry axis is around 28 eV which is consistent with the experimental values reported by Lafleur et al [68] and later by Correyero et al [55] for 0.2 mg s⁻¹ and, 20 and 30 W of absorbed power, respectively. T_e goes below 4.3 eV at the lateral wall of the thruster and close to 0.2 eV at the top of the external wall.

Further insight on the electron dynamics can be gained by inspecting the electron pressure p_e = n_e T_e, shown in figure 5(f). This variable is directly related to the thrust of the device, and with the quasi-steady electric field and azimuthal plasma currents, as discussed further below. While n_e displays a minimum in the neighborhood of the axis of symmetry and T_e presents a strong radial gradient, p_e shows a relatively smooth behavior in the whole discharge. Looking at z = const sections of the thruster, the peak pressure takes place at a mid radius inside the annular chamber. Further downstream, in the MN, the maximum pressure is found closer to the axis of symmetry.

The quasi-steady plasma potential ϕ in figure 5(d) presents an axial fall of roughly 50 V within the simulation domain. It must be noted that the expansion will continue further downstream of this domain to infinity with collisionless electron cooling and further potential drop. The maximum value of ϕ is found inside the discharge chamber, close to the inner rod element of the thruster, right before the quasineutral pre-sheath decreases its magnitude toward the thin non-neutral wall sheaths. The potential decreases axially, along the MN. Note that the position of the maximum does not coincide with either the maximum of n_e, T_e, nor p_e. In the direction perpendicular to the magnetic lines, and always in the quasineutral domain, the potential is essentially flat, except close to the inner rod element where it drops several tens of volts. A minor, secondary maximum of ϕ can be seen close to the lateral inner corner of the discharge chamber, which is attributed to this corner region of being essentially magnetically isolated from the rest of the plasma by the applied field lines, which means that electron transport into this region is severely limited; the small potential rise in this part of the device enhances electron transport from the rest of the plasma toward this region. Similarly, the small increase of ϕ at the top-left corner of the plume domain is attributed to the lack of electrons in this second magnetically-isolated region, although in this case it is considered a numerical side-effect of the finite simulation domain and the presence of the lateral plume boundary.

Another aspect of interest of the ϕ-map is that the potential fall along each magnetic line is not proportional to the corresponding value of T_e, as could perhaps be expected. If this were true, a higher potential fall would be observed in the magnetic tube close to the inner rod element of the thruster. On the contrary, the potential fall is smooth and roughly the same across magnetic lines, indicating that the expansion is a global phenomenon for the whole plasma, despite electrons being well-magnetized. This behavior is consistent with the profiles of n_e and T_e of figures 5(c) and (e). Indeed, to leading order, the projection of equation (2) along 1_∥ yields a Boltzmann-like relation for each magnetic line within the simulated domain,

$$\begin{equation}\frac{{n}_{\mathrm{e}}}{{n}_{\mathrm{e}0}}\simeq \mathrm{exp}\frac{e\left(\phi -{\phi }_{0}\right)}{{T}_{\mathrm{e}}},\end{equation} \tag{ 23 }$$

where the subindex 0 represents values at a reference location for each line (e.g. upstream). As Δϕ = ϕ − ϕ₀ is roughly identical for all magnetic lines, the relative drop of n_e is smaller along the lines where T_e is larger.

The singly-charged ion meridian velocity u _i1 and their streamlines are shown in figure 5(h). While electrons are magnetized, the heavy ions are essentially unmagnetized and thus their dynamics is dominated mainly by the quasi-steady potential ϕ, which accelerates them downstream. Observe that ion streamlines (which roughly correspond to individual ion trajectories, since ions are relatively cold) do indeed not coincide with magnetic lines, and that the ion flow detaches inwardly from the magnetic field [69], resulting in less divergence than the expected one for fully-magnetized ions [70]. The resulting mean velocities at the last section of the simulation domain range from 3 km s⁻¹ (at the axis) to 6 km s⁻¹ (at the periphery). This difference is attributed to the location where the ions are created and the path they follow inside the ionization chamber before entering the MN region. For comparison, laser induced fluorescence (LIF) measurements of the ECRT operating at nominal conditions [55] indicate ion velocities of about 7 km s⁻¹ and 11 km s⁻¹ at 4 and 12 cm downstream the thruster exit, respectively. In addition, laser induced fluorescence (LIF) measurements obtained at Michigan [71] with a similar prototype report velocities around 9 km s⁻¹ 4 cm downstream the thruster exit plane, at background pressure of 13 μTorr.

As a derived quantity, the meridian Mach number of singly-charged ions, based on the approximate local sound velocity, ${M}_{\mathrm{i}1}=\vert {\boldsymbol{u}}_{\mathrm{i}1}\vert /\sqrt{{T}_{\mathrm{e}}/{m}_{\mathrm{i}1}}$, features its maximal values (around M_i1 = 3) in the peripheral region of the plume, while at its core M_i1 is only slightly larger than 1. This radial variations are due to the differences in u _i1, compounded with the variations in electron temperature T_e. Arguably, the supersonic expansion in the MN of the device is uneven. It seems plausible that the presence of the rod element in the ECRT discharge chamber induces a dip in the ion Mach number in its wake. This behavior was seen to be robust against the sensitivity analysis carried out, and will be subject of future studies.

Lastly, it can be observed that inside the discharge chamber part of the ion streamlines precipitate onto the walls of the thruster guided by the electrostatic potential fall near those surfaces, giving rise to wall losses, which are discussed in section 3.3.

The quasi-steady state azimuthal electron and ion currents are an essential aspect of the operation of the device, as their interaction with the applied magnetic field B ₀ is responsible for the generation of magnetic thrust and the cross-field confinement of the plasma. Figure 6 displays the azimuthal current of electrons j_θe = −en_e u_θe, which is more than two orders of magnitude greater than the azimuthal current of any ion species.

Zoom In Zoom Out Reset image size

To leading order, in the perpendicular direction the electron pressure gradient −∇p_e is compensated by the electrostatic force density en_e∇ϕ and the magnetic force, j_θe B₀ 1_⊥ (see equation (2)),

$$\begin{equation}\frac{\partial {p}_{\mathrm{e}}}{\partial {\mathbf{1}}_{\perp }}\simeq e{n}_{\mathrm{e}}\frac{\partial \phi }{\partial {\mathbf{1}}_{\perp }}+{j}_{\theta \mathrm{e}}{B}_{0}.\end{equation} \tag{ 24 }$$

Given that the quasi-steady electric potential is essentially flat across the magnetic field lines, the current j_θe is responsible for most of the cross-field confinement of the electron pressure p_e. Consequently, the sign of j_θe switches mid-radius, close to where the maximum p_e is found in the perpendicular direction, following the sign of −∇p_e.

The reaction to the axial magnetic force density, j_θe B_r0, is felt on the coils of the thruster, and is termed magnetic thrust. Since in the present simulation B_z0, B_r0 > 0 in the domain, only negative j_θe results in positive thrust generation. It can be seen that the outermost part of the plasma is the main contributor to magnetic thrust, while the electrons close to the inner rod element of the device cause some negative magnetic drag in order to confine the electron pressure away from this element. This is a consequence of the maximum electron pressure being located at an intermediate radius rather than the axis of symmetry, and is likely an effect of the presence of the inner rod element. Downstream of this rod element, a small localized region with j_θe < 0 is found.

The solution to the high-frequency EM fields depends on the excitation frequency ω, the applied magnetic field B ₀, the electron density n_e, and the collisionality ν_e. In particular, the first three variables determine the field propagation regimes, the cut-offs, and the resonances. The last variable, ν_e, which ranges from a maximum value of 4.5 × 10⁸ Hz inside the source and decreases downstream, affects mainly to the homogeneity of the power absorption map. Notably, while B ₀ is an input to the simulation, n_e and ν_e depend dynamically on the plasma transport solution.

The different EM parametric regions in steady-state are shown in figure 7(a), following the numerical nomenclature (I to VIII) of Stix [59]. The right-hand circularly polarized wave (R) resonates when ω = eB₀/m_e which determines the critical value of B₀ (875 G) and the location of the ECR surface. This value separates regions I through V from regions VI through VIII. The plasma frequency cut-off occurs when ω = ω_pe, which determines the critical value of n_e (=7 × 10¹⁶ m⁻³). Regimes I, II, III, and VI take place below this value of n_e, and are considered underdense regimes. Regimes IV, V, VII and VIII are overdense regimes. Region I (low B₀, low n_e) is topologically equivalent to propagation in vacuum, with both R and L (i.e. left-hand circularly polarized) propagating modes. It is separated from region II by the R wave cut-off. In region II, only L waves propagate. As B₀ and n_e increase, eventually the UHR is crossed into region III, where two propagating solutions exist at directions other than B ₀. Region IV (low B₀, high n_e) exhibits propagating L waves only. In region V (separated from region IV by the L wave cut-off) no propagating solutions exist, and all EM fields are evanescent. Region VI (high B₀, low n_e) features both R and L polarizations. Regions VII and VIII (high B₀, high n_e) feature whistler R waves, which propagate in directions close to B ₀; additionally, region VII also has regular L-waves.

Zoom In Zoom Out Reset image size

As it can be noticed in figure 7(a), in this ECRT simulation all these regions (i.e. I through VIII) are present within the simulation domain and in close proximity to each other. The contour lines detail the location of the different cutoffs and resonances, which act as the boundary surfaces of the different propagation regimes. In particular, the presence of the ECR and UHR stands out. Whereas the location of the ECR is fixed only by the magnetic field, the UHR is determined by a combination of both the magnetic field and the electron density.

As the EM power enters the discharge chamber through the coaxial line, it enters regions VI–VIII (upstream of the critical B₀), and soon reaches the ECR surface. Most of the discharge chamber is overdense (ω_pe > ω), and therefore downstream of this surface there is a region V where only evanescent fields are allowed. However, a thin low-density channel exists near the inner wall corresponding to regions III and IV. Power can propagate along this channel, and also tunnel through a short evanescent region. Further downstream, as both B₀ and n_e decrease, all other regions are crossed, until eventually region I is reached.

The existence of all these regions implies wide changes of the local wavelengths in the domain. Figure 7(b) illustrates this by plotting the minimum wavelength in the perpendicular and parallel directions for all propagating modes. The smallest wavelengths occur near the ECR and the UHR surfaces, and the W-module mesh has been refined accordingly, as described in section 2.3 and appendix A.

The amplitudes and complex phases of the three components of the fast electric field E are shown in figure 8. The resulting electric field inside the coaxial cable is a strong radial field, characteristic of the TEM mode. As the power flows into the plasma, the radial component of the EM wave electric field ${\tilde {E}}_{\mathrm{r}}$ still dominates on both sides of the ECR surface, and continues to be large in the thin regions III and IV that exist around the rod element of the thruster. It can be observed that an axial electric field ${\tilde {E}}_{z}$ also develops before the ECR surface and downstream from the rod element decaying after the resonance, and that the azimuthal field ${\tilde {E}}_{\theta }$ plays only a minor role. Small magnitude fields also exist downstream, around the UHR surface.

Zoom In Zoom Out Reset image size

The phase plots in figure 8 indicate that no major wave structures exist in the domain, which would be evidenced by alternating phase structures, given its small size compared to the electric length in vacuum; the only exception is in the narrow propagation channel close to the inner rod element, where a short-wavelength mode is observed in ${\tilde {E}}_{z}$. No major standing-wave structures are observed within the simulation domain, indicating that reflection at the free surfaces downstream is small.

The power absorption density Q_a is shown in figure 9. Absorption occurs mainly at around the ECR surface, as expected. The power absorption has its maximum at the ECR, near the inner rod element of the thruster. Some power absorption also takes place around the inner rod, in the propagating channel. Overall, most of the power is deposited in the magnetic tubes close to the inner rod, which is also where electron temperature features its peak values. The absorption at the downstream UHR surface is non-zero, but orders of magnitude smaller.

Zoom In Zoom Out Reset image size

The absorption map can be explained by inspecting the magnitude of the different terms in equation (22). As the applied magnetic field is mostly axial in the device, ${\tilde {E}}_{{\Vert}}\simeq {\tilde {E}}_{z}$ and ${\tilde {E}}_{\perp }\simeq {\tilde {E}}_{\mathrm{r}}$. At the ECR surface, the imaginary part of $\mathcal{S}$ features its maximum and it is two orders of magnitude greater than the imaginary part of $\mathcal{P}$. As $\vert {\tilde {E}}_{\mathrm{r}}\vert \gg \vert {\tilde {E}}_{\theta }\vert$, the power absorption density Q_a is therefore dominated by the term

$$\begin{equation}{Q}_{\mathrm{a}}\approx \frac{\omega {{\epsilon}}_{0}}{2}\mathfrak{I}\left(\mathcal{S}\right)\vert {\tilde {\boldsymbol{E}}}_{\mathrm{r}}{\vert }^{2}.\end{equation} \tag{ 25 }$$

This term reaches its peak value close to the exit of the coaxial line, where the radial electric field is larger and the ECR surface is found, and consequently this is where the higher Q_a occurs. Lastly, it is noted that while absorption at the resonance requires ν_e > 0, the actual value of this parameter only influences the thickness of the resonance layer, while the integrated power deposited in it is essentially unaffected [60].

Figure 10(a) displays the particle flux to the thruster walls for each species, as a function of arc-length variable ζ, introduced in figure 1. This variable starts at the tip of the inner rod, goes around clockwise along the thruster wall meridian section, finishing at the end of the external wall.

Zoom In Zoom Out Reset image size

The internal thruster walls sustain a much larger particle flux than the external walls, differing in several orders of magnitude. Singly- and doubly-charged ion particle fluxes evolve similarly with ζ, being the former always greater than the latter. Singly-charge ion particle flux features values between 10¹⁹ and 10²¹ m⁻² s⁻¹ flux, peaking around the middle of the rear wall, where the injector is located and ions density is maximum. Neutral flux to the walls is shown in the figure for comparison. Amongst all the regions inside the thruster chamber, the corner region between the rear and lateral wall collects the least amount of ions. All corners exhibit a discontinuity motivated by the change in normal vector but this corner region features enhanced magnetic shielding. As a result insufficient electron temperature in the vicinity of this region, ionization is poor there, and combined by the strong magnetization results in a repelling electrostatic field that limits the ion flux to the walls significantly. The primary electron flux follows a similar behavior to the singly-charged ions both at the lateral wall and backplate. Nevertheless, approaching the inner rod wall, the flux increases by more than an order of magnitude reaching values of about 10²² m⁻² s⁻¹.

At each location along the (dielectric) walls, the net flux of negative and positive charges to the wall must be equal. The large difference between the fluxes of electrons and ions toward the inner rod wall is indicative of the large SEE that takes place there, where indeed the yield δ_ws, motivated by the considerable electron temperature, is close to unity. The difference of the fluxes decreases to near zero along the back, lateral, and external walls, where electrons are cooler and SEE becomes less important.

Figure 10(b) shows the total net energy flux to the walls of each species. This magnitude is relevant to understand the thermal loads endured by the exposed walls. The maximum energy load is exhibited at the inner rod. Electrons dominate the energy flux to this element, due to the high SEE, and are comparable to the ion energy flux in other walls. Singly-charged ions have a larger energy flux with respect to doubly-charged ions everywhere due to their significantly higher density. The maximum ion energy flux is 0.8 W/cm², and occurs at the injector location, near the region where most ionization takes place. As could be expected, the corner region of the thruster is clearly well shielded magnetically, not only from a particle flux viewpoint as reported above, but also from an energy viewpoint.

Finally, figure 10(c) portraits the average impact energy per particle of each species. This variable is relevant in erosion studies, as the impact energy must be greater than a material dependent threshold energy for sputtering to take place. In this case, we see that doubly-charged ions are responsible for the highest per-particle energy deposition to the walls. This is due mainly to their double-charge acceleration in the wall sheaths; indeed, the sheath accounts for most of the impact energy for both singly- and doubly-charged ions everywhere. Energies ranging from roughly 10–20 eV to 200–300 eV are computed. The maximum impact energy takes place at the end part of the rod element, where experiments report the largest and fastest erosion under operation. Specific electron energy is small except on the inner rod element. Neutrals impact energy is negligible.

The thrust F of a MN-based thruster is generated partially inside the source, but mainly in the MN region, which extends to infinity. In the present simulation, F is computed by integration of the plasma momentum on the free boundaries of the domain ∂Ω_f,

$$\begin{equation}F=\sum\limits _{s}{\int }_{\partial {{\Omega}}_{\mathrm{f}}}\left({m}_{s}{n}_{s}{u}_{zs}{\boldsymbol{u}}_{s}\cdot {\mathbf{1}}_{\mathrm{n}}+{n}_{s}{T}_{s}{\mathbf{1}}_{z}\cdot {\mathbf{1}}_{\mathrm{n}}\right)\mathrm{d}S,\end{equation} \tag{ 26 }$$

where the sum on s extends to all species. The thrust F can be divided in two contributions:

$$\begin{equation}{F}_{\mathrm{p}}=\sum\limits _{s}{\int }_{\partial {{\Omega}}_{\mathrm{w}}}\left({m}_{s}{n}_{s}{u}_{zs}{\boldsymbol{u}}_{s}\cdot {\mathbf{1}}_{\mathrm{n}}+{n}_{s}{T}_{s}{\mathbf{1}}_{z}\cdot {\mathbf{1}}_{\mathrm{n}}\right)\mathrm{d}S,\end{equation} \tag{ 27 }$$

and

$$\begin{equation}{F}_{\mathrm{m}}={\int }_{{\Omega}}-{j}_{\theta }{B}_{\mathrm{r}}\enspace \mathrm{d}V,\end{equation} \tag{ 28 }$$

which equal the sum of all the axial dynamic pressure force of each species s on the thruster walls ∂Ω_w and the axial magnetic force produced by the azimuthal plasma current on the thruster magnets, respectively [21]. It should be noted that the expansion continues further downstream, beyond the computational domain, and therefore a small part of the generated thrust is missed by the simulation. Indeed, the major part of ion acceleration in this thruster occurs in the first few cm after the thruster chamber exit [55], and therefore is captured here. The specific impulse is calculated from this magnitude as ${I}_{\mathrm{s}\mathrm{p}}=F/\dot {m}$.

The microwave power entering the thruster, P_f, is partially absorbed by the plasma, P_a, and partially reflected back through the coaxial cable, P_r, thus

$$\begin{equation}{P}_{\mathrm{f}}={P}_{\mathrm{a}}+{P}_{\mathrm{r}}.\end{equation} \tag{ 29 }$$

Notice that in the present simulations there are no free-space radiation losses; this is representative of the operation of the device inside a laboratory vacuum chamber, where the radiated power is confined.

The reflected power is obtained directly from the voltage standing wave ratio (VSWR), computed as the ratio of the maximum and minimum values of the radial electric field in the coaxial cable, and is directly related to the power coupling efficiency η_p as

$$\begin{equation}{\eta }_{\mathrm{p}}=1-\frac{{P}_{\mathrm{r}}}{{P}_{\mathrm{f}}}=1-{\left(\frac{\text{VSWR}-1}{\text{VSWR}+1}\right)}^{2}.\end{equation} \tag{ 30 }$$

Observe that the reflected power cannot be considered a priori a loss mechanism, as this power can be recirculated back to the thruster quite efficiently by impedance matching in the microwave circuit design. Nevertheless, a large reflection ratio would be symptomatic of an inadequate power coupling to the plasma.

Following [47], the power balance for all plasma species can be expressed as

$$\begin{equation}{P}_{\mathrm{a}}={P}_{\text{exc}}+{P}_{\text{ion}}+{P}_{\text{wall}}+{P}_{\mathrm{p}},\end{equation} \tag{ 31 }$$

where P_exc and P_ion are the power spent in excitation and ionization of the propellant; P_wall is the kinetic, thermal, and heat flux power lost to the walls; and P_p is the power flux of all species through the free boundaries. Consequently, we can define the loss ratios _exc = P_exc/P_a, _ion = P_ion/P_a, and _wall = P_wall/P_a.

The (overall) thruster efficiency is computed as

$$\begin{equation}{\eta }_{\mathrm{F}}=\frac{{F}^{2}}{2\dot {m}{P}_{\mathrm{a}}}.\end{equation} \tag{ 32 }$$

This efficiency can be approximately decomposed as

$$\begin{equation}{\eta }_{\mathrm{F}}\approx {\eta }_{\mathrm{u}}{\eta }_{\mathrm{e}}{\eta }_{\mathrm{c}}{\eta }_{\mathrm{d}}=\frac{{\dot {m}}_{\mathrm{i}\infty }}{\dot {m}}\frac{{P}_{\mathrm{p}}}{{P}_{\mathrm{a}}}\frac{{P}_{\mathrm{i}}}{{P}_{\mathrm{p}}}\frac{{P}_{z\mathrm{i}}}{{P}_{\mathrm{i}}},\end{equation} \tag{ 33 }$$

where ${\eta }_{\mathrm{u}}={\dot {m}}_{\mathrm{i}\infty }/\dot {m}$ is the utilization efficiency, i.e. the fraction of propellant mass flow rate that reaches the domain free boundaries as ions (${\dot {m}}_{\mathrm{i}\infty }$); η_e = P_p/P_a refers to the energy efficiency, i.e. the fraction of absorbed power that becomes plume kinetic and thermal power (P_p); η_c = P_i/P_p is the conversion efficiency, i.e. the portion of the plume power in the form of kinetic ion power (P_i); finally, η_d = P_zi/P_i is the divergence efficiency, i.e. the fraction of ion kinetic power which in the axial direction (P_zi), thus generating thrust. The difference in equation (33) between η_F and the product η_u η_e η_c η_d is due to the difference between F_i and F, which in a completely developed expansion is expected to be small. However, in the present simulation with finite domain, the electron contribution to thrust accounts for up to a 25% (see figure 4) at the free downstream boundary. This suggests that the expansion is incomplete by z = 4 cm, and that additional thrust is to be gained further downstream as the MN continues to convert electron power into directed ion kinetic power.

In addition to these partial efficiencies, another relevant quantity is the production efficiency, defined as the ratio of ion flow rate at the free boundaries with respect to that at all the simulation boundaries including thruster walls,

$$\begin{equation}{\eta }_{\text{prod}}=\frac{{\dot {m}}_{\mathrm{i}\infty }}{{\dot {m}}_{\mathrm{i}}},\end{equation} \tag{ 34 }$$

which is closely related to the number of times each atom undergoes reionization, and characterizes the plasma losses to the walls.

Table 2 displays the performance figures of the simulation. For comparison purposes, the available experimental data from ONERA thruster are also shown. The total thrust produced by the thruster in the simulation is 0.84 mN, and the specific impulse 428 s, which show a great agreement with the experimental measurements. The fraction of magnetic thrust to total thrust amounts to approximately 62%, which again agrees well with the available experimental data [72]: 60% for the device operating at 0.2 mg s⁻¹ and 40 W.

Table 2. Thruster performances. Values at the thruster (^*) are computed from the data reported at reference [11] (before the vacuum chamber feedthrough) with the estimated 2 dB losses between that point and the thruster, as indicated in that reference.

Parameter	Name	Units	UC3M	ONERA [11]
$\dot {m}$	Mass flow rate	mg s−1	0.2	0.2
Pf	Input power	W	30	53.7*
F	Thrust	mN	0.840	0.850
Fp	Pressure thrust	mN	0.322	—
Fm	Magnetic thrust	mN	0.518	—
Pa	Absorbed power	W	26.9	27.9*
Pr	Reflected power	W	3.1	4.28*
Isp	Specific impulse	s	428	429
Ii	Ion current	A	0.082	0.065
ηF	Thrust efficiency	%	6.6	6.5*
ηprod	Production efficiency	%	40.9	—
ηu	Utilization efficiency	%	50	45.1
ηe	Energy efficiency	%	25.3	—
ηc	Conversion efficiency	%	32.2	—
ηd	Divergence efficiency	%	91.2	—
exc	Excitation losses	%	4.8	—
ion	Ionization losses	%	7.0	—
wall	Material wall losses	%	63.0	—
VSWR	Voltage ratio	—	1.8	2.15*
ηp	Power coupling eff.	%	91.8	86.7*

In reference [11], efficiencies are computed with respect to the transmitted power (i.e. η_p P_f) and that this power includes the power losses in the cables, the feed-through, the DC block and the connectors/adapters, which are of the order of 2 dB. Table 2 shows the efficiencies and powers, taking into account the reported attenuation losses between the measurement location (prior to the tank) and the thruster. The thrust efficiency is very similar in the simulations and the experiments, around 6.6%.

The main inefficiency arises from plasma losses to the walls. First, the fraction of generated ions that become plume ions, given by η_prod, is quite low, about 41%, while the rest of them are recombined at the walls. Second, the energy efficiency is rather low, η_e ≈ 25%. While excitation and ionization losses (measured by _exc and _ion) are arguably small, about two-thirds of the absorbed power (_wall) is lost to the walls, transported by the ions and electrons that impact on them. In fact, the integrated energy fluxes (see section 3.3) over the inner thruster walls areas adds up to 16.9 W, where 51% is lost through the inner rod surface, 40% through the backplate, 9% through the lateral wall and 0.1% through the external wall. Large losses at the magnetically unshielded backwall are well-known from HPT analyses [21, 73] and the way to reduce them would be to shield that wall (without altering much the wave propagation). Magnetic shielding seems to operate rather well at the lateral wall; it does not behave well at the rod due to the large electron temperatures around it, leading to very large thermal loads in that thin element, which amounts for only a 6% of the total inner thruster walls area. To reduce this problem should be a central task of any design optimization of the thruster chamber.

The second main source of inefficiency is related to the meager utilization efficiency, a 50% in the simulation, and also in good agreement with the experimentally-obtained value. These values are also common in HPTs [74] and are likely due to the large wall recombination. The consequence is that a large fraction of the propellant is leaving the source as neutrals, generating essentially no thrust.

On a positive note, the divergence efficiency η_d of 91.2% shows that the kinetic power of ions leaving the domain is mainly axial and that good magnetic detachment takes place in the MN. On the other hand, and in line with the discussion above, the rather low conversion efficiency η_c evidences that the expansion is still incomplete in the simulation domain.

The experimentally-reported power coupling efficiency η_p is close to unity already in early designs [4] and in contemporary ones [18], where they report a 95%. These values are measured between the microwave generator and the microwave feedthrough on the tank, and therefore do not record the attenuation in the line from the tank feedthrough to the thruster. With the estimated power attenuation factor of 2 dB given in [11], the actual coupling efficiency at the thruster entrance for the prototype simulated here is 87%. This value is close to the computed 92% coupling efficiency.