2.1. Particle-in-cell method

In the PIC method, the plasma is treated as an ensemble of computational particles (Dawson Reference Dawson1983; Birdsall & Langdon Reference Birdsall and Langdon2018; Hockney & Eastwood Reference Hockney and Eastwood2021). The computational particles move on the simulation domain, where, at every time iteration, the charge of each particle is distributed to a dedicated mesh. The charge is distributed according to particle positions by a linear interpolation scheme. Once the charge density $\rho$ is computed on the nodes of the mesh, the electric field $\boldsymbol {E}$ can be integrated from the Poisson's equation. The full set of Maxwell equations has not been solved since the amount of RF power deposited outside the plasma source is a small fraction (usually less than 1 %) of the total power coupled to the plasma (Magarotto & Pavarin Reference Magarotto and Pavarin2020). This means that the propagation of RF waves in the plume is expected to have a minor influence on its dynamics. At the same time, currents (mainly azimuthal; Takahashi Reference Takahashi2019) induced by the presence of the magnetic nozzle are sufficiently low that the self-induced magnetic field produced can be neglected (errors less than 1 % are expected; Takahashi Reference Takahashi2019). Therefore the wave propagation in the plasma plume has been neglected and the magnetic field $\boldsymbol {B}$ is static and equal to that imposed by the presence of the magnetic nozzle (Magarotto et al. Reference Magarotto, Melazzi and Pavarin2020b). Consequently, the Maxwell equations reduce to the Poisson equation:

(2.1)\begin{equation} \tilde{\varepsilon}\Delta \phi ={-}\rho , \end{equation}

with $\phi$ the electrostatic potential and $\boldsymbol {E} = -\boldsymbol {\nabla } \phi$. The Poisson equation is solved via a finite element method with a preconditioned conjugate algorithm to obtain $\phi$ and consequently the field $\boldsymbol {E}$ (Rogier & Volpert Reference Rogier and Volpert2021).

The force $\boldsymbol {F_p}$ acting on every particle is given by Lorenz force:

(2.2)\begin{equation} \boldsymbol{F_p} = q_p (\boldsymbol{E} +\boldsymbol{v_p}\times\boldsymbol{B}), \end{equation}

where $\boldsymbol {v_p}$ is the speed of the particle and $q_p$ its charge. Force $\boldsymbol {F_p}$ is computed interpolating fields on the position of each particle, and the associated trajectory is then integrated via a Runge–Kutta Cash–Karp method (Sarrailh et al. Reference Sarrailh, Mateo-Velez, Hess, Roussel, Thiebault, Forest, Jeanty-Ruard, Hilgers, Rodgers and Cipriani2015). This sequence of operations is repeated at each time iteration (see figure 2), self-consistently evolving the particles and the electric field states. The simulations are performed until steady state is reached.

Figure 2. Time loop in PIC code.

It is worth noting that the Runge–Kutta Cash–Karp method shows a considerable propagation of numerical error integrating the trajectory of a charged particle (Qin et al. Reference Qin, Zhang, Xiao, Liu, Sun and Tang2013). A correction factor in the pusher algorithm has been added to reduce the energy numerical drift.

Finally, the effect of collisions has been neglected. This approach is proven to provide a valuable preliminary insight on the dynamics of the plume for typical operation conditions of MEPTs (Martinez-Sanchez et al. Reference Martinez-Sanchez, Navarro-Cavallé and Ahedo2015). Nonetheless collisions will be implemented in future works to estimate quantitatively the influence that phenomena such as the doubly trapped electrons have on the plume. The latter in fact are produced by collisional effects (Zhou et al. Reference Zhou, Sánchez-Arriaga and Ahedo2021).

2.2. Boundary conditions

The simulation domain and the boundary surfaces are represented in figure 3. The boundary conditions that have been imposed for the closure of (2.1) are the following.

1. (i) Dirichlet condition on the surfaces of the spacecraft, i.e.

   (2.3)\begin{equation} \phi = \begin{cases} \phi_o, & \text{ on the thruster outlet surface (particle source)},\\ 0, & \text{ on other surfaces}. \end{cases} \end{equation}

2. (ii) Robin condition on the external boundary (Thiebault et al. Reference Thiebault, Jeanty-Ruard, Souquet, Forest, Mateo-Velez, Sarrailh, Rodgers, Hilger, Cipriani and Payan2015):

   (2.4)\begin{equation} \frac{\textrm{d}\phi}{\textrm{d}\boldsymbol{k}} +\frac{\boldsymbol{r} \cdot\boldsymbol{k}}{r^2}\phi = 0, \end{equation}
   with $\boldsymbol{k}$ the external normal to the surface and $\boldsymbol{r}$ the vector going from the geometric centre of the particle source to the external boundary surface. Equation (2.4) derives from the assumption that at infinity $\phi =0$ and that $\phi \propto 1/r$ (Thiebault et al. Reference Thiebault, Jeanty-Ruard, Souquet, Forest, Mateo-Velez, Sarrailh, Rodgers, Hilger, Cipriani and Payan2015).

Figure 3. Simulation domain: the boundary surfaces are coloured in brown, the particle source in green. The simulation domain corresponds to the grey zone.

As far as particle motion is concerned, the definition of an electron energy boundary condition is required to simulate properly the dynamics of the plasma plume (Li et al. Reference Li, Merino, Ahedo and Tang2019). Indeed, to maintain the current-free condition (Takahashi Reference Takahashi2019), a positive potential drop occurs between the thruster outlet and ‘infinity’. The simulation should take into account the electrons that do not have enough kinetic energy to escape this potential drop and therefore return to the plasma source or become trapped in the plume. This is fundamental also from a computational standpoint to avoid the so-called ‘numerical pump instability’ (Brieda Reference Brieda2005). Despite the finite computational domain dimensions, these trapped electrons can be taken into account by adding an energy-based boundary condition. Energy $E_{\textrm {tot}}$ denotes the total energy of the electrons at the external boundary:

(2.5)\begin{equation} E_{\textrm{tot}} = \tfrac{1}{2}m_p v_p^2 +q_p \phi_b , \end{equation}

with $m_p$ the mass of the particle and $\phi _b$ the external boundary potential. Electrons that reach the external boundary surfaces with a negative $E_{\textrm {tot}}$ are reflected specularly; in contrast, electrons with a positive $E_{\textrm {tot}}$ are absorbed. In fact the former do not have enough energy to achieve infinity with $|v_p|\geq 0$. Ions that reach the external boundary are absorbed. Ions and electrons that reach the surfaces of the spacecraft, including the thruster outlet surface, are absorbed too.

2.3. Particle sources

The computational particles are ejected from the surface of the thruster named ‘thruster outlet’ in figure 3. Hereafter, an $xyz$ Cartesian coordinate system with the $z$ axis normal to the thruster outlet is considered; $z$ is also the axis of symmetry for the computational domain. As detailed in §§ 2.3.1 and 2.3.2, a Maxwellian distribution function is assumed for all the species at the thruster outlet. This hypothesis is justified by the value of the plasma density which is orders of magnitude higher in the production stage than in the plume, so a thermodynamic equilibrium is more likely (Cichocki et al. Reference Cichocki, Domínguez-Vázquez, Merino and Ahedo2017; Li et al. Reference Li, Merino, Ahedo and Tang2019; Magarotto & Pavarin Reference Magarotto and Pavarin2020; Magarotto et al. Reference Magarotto, Manente, Trezzolani and Pavarin2020a). Moreover, the assumption of a Maxwellian distribution function at the thruster outlet is supported by experiments performed on helicon sources (Chen & Blackwell Reference Chen and Blackwell1999). Hereafter, random variables are written in upper case, and particular realizations of such variables are written in the corresponding lower case.

2.3.1. Electron source

The speed $V_x$ along the $x$ axis and the speed $V_y$ along the $y$ axis of the ejected electrons are random variables that follow a normal distribution of mean $\mu = 0$ and standard deviation $\sigma$ given by

(2.6)\begin{equation} \sigma = \sqrt\frac{k_BT_e}{m}, \end{equation}

with $T_e$ the electron temperature expressed in kelvin, $k_B$ the Boltzmann constant and $m$ the electron mass. The speed along the $z$ axis $V_z$ of the ejected electrons has the following probability density function (Li et al. Reference Li, Merino, Ahedo and Tang2019):

(2.7)\begin{equation} f_{V_z}(v_z)=\begin{cases} 0 & \text{ if } v_z<0,\\ \dfrac{2}{\sigma\sqrt{2{\rm \pi}}} e^-\dfrac{v_z^2}{2\sigma^2} & \text{ if } v_z\geq 0. \end{cases} \end{equation}

Equation (2.7) represents a flux of electrons ejected from the thruster outlet with a Maxwellian distribution function. Instead the flux of electrons returning at the same surface ($v_z<0$) is an output of the simulation that cannot be imposed a priori. It is easy to show that the norm of the speed of the ejected electrons, denoted as $V$, reads (Bartsch Reference Bartsch2014)

(2.8)\begin{equation} f_{V}(v)= \sqrt{\frac{2}{\rm \pi}}\frac{v^2}{\sigma^3} e^-\frac{v^2}{2\sigma^2}. \end{equation}

In particular, from the definition of electron current (Chen et al. Reference Chen1984), it is possible to compute the electron density $n_e$ at the thruster outlet as

(2.9)\begin{equation} n_e = \frac{\dfrac{|I_e|}{|e|S}}{|E[\widetilde{V_z}]|}, \end{equation}

with $e$ the electron charge, $S$ the surface area of the thruster outlet, $I_e$ the electron net current at the thruster outlet and $|E[\widetilde {V_z}]|$ the expected value of the electron speed $\widetilde {V_z}$ at the thruster outlet. The electron speed $\widetilde {V_z}$ at the thruster outlet corresponds to electrons both ejected and absorbed at this surface.

2.3.2. Ion source

Ions are ejected from the thruster outlet with the following velocities:

(2.10)\begin{equation} \left. \begin{aligned} V_x & = V_n +u_B \sin{\varTheta}\cos{\varLambda},\\ V_y & = V_n +u_B \sin{\varTheta}\sin{\varLambda},\\ V_z & = V_n +u_B \cos{\varTheta}, \end{aligned} \right\} \end{equation}

where $V_n$ is a random variable that follows a normal distribution of mean $\mu = 0$ and standard deviation $\sigma$ given by

(2.11)\begin{equation} \sigma = \sqrt\frac{k_BT_i}{M}, \end{equation}

with $T_i$ the ion temperature, $M$ the ion mass and $u_B$ the Bohm speed:

(2.12)\begin{equation} u_B = \sqrt{\frac{k_BT_e}{M}}. \end{equation}

Here $\varTheta$ is a random variable that follows a continuous uniform distribution in the interval $[0, \theta _t ]$. Angle $\theta _t$ corresponds to the nozzle divergence angle (see figure 3). It is associated with the presence of a conical source (Charles, Boswell & Takahashi Reference Charles, Boswell and Takahashi2013) or an expansion nozzle (Manente et al. Reference Manente, Trezzolani, Magarotto, Fantino, Selmo, Bellomo, Toson and Pavarin2019). Considering typical MEPTs, in this study the assumption $\theta _t < 45^{\circ }$ has been made (Bellomo et al. Reference Bellomo, Di Fede, Magarotto, Manente, Pavarin, Trezzolani, Selmo, De Carlo, Toson and Minute2021). Parameter $\varLambda$ is a random variable that follows a continuous uniform distribution in the interval $[0, 2{\rm \pi} ]$. The expected value $E[V_z]$ and the variance $\textrm {Var}[V_z]$ of $V_z$ are

(2.13)\begin{gather} E[V_z] = \frac{u_B\sin{\theta_t}}{\theta_t}, \end{gather}

(2.14)\begin{gather} \textrm{Var}[V_z] = \sigma^2+\frac{u_B^2}{2}\left(\frac{\sin{2\theta}}{2\theta_t}+1 -2\left(\frac{\sin\theta_t}{\theta_t}\right)^2\right). \end{gather}

In typical MEPT configurations the plasma is mesothermal (i.e. $T_e \gg T_i$) (Balkey et al. Reference Balkey, Boivin, Kline and Scime2001); therefore $u_B$ is higher than $\sigma$ and it is reasonable to assume that

(2.15)\begin{equation} P(V_z < 0) \approx 0, \end{equation}

where $P(V_z < 0)$ is the probability that ions are ejected with a negative $v_z$. The ion density $n_i$ at the thruster outlet can be computed as

(2.16)\begin{equation} n_i = \frac{\dfrac{|I_i|}{q_iS}}{|E[\widetilde{V_z}]|}, \end{equation}

with $q_i$ the ion charge, $I_i$ the ion net current at the thruster outlet and $E[\widetilde {V_z}|$ the expected value of ion speed $\widetilde {V_z}$ at the thruster outlet. The ion speed $\widetilde {V_z}$ at the thruster outlet is the speed corresponding to ions both ejected and absorbed by the thruster outlet. The outlet potential $\phi _o$ is larger than the boundary potential $\phi _b$; therefore a negligible amount of ions return to the thruster outlet and it is reasonable to assume that $\widetilde {V_z} \approx V_z$.

2.4. Control loop

For the sake of simplicity hereafter, it has been considered that the ion population is composed of only one species. Considering typical operating conditions of MEPTs, this assumption is reasonable where the thruster is operated with noble gases (Chabert et al. Reference Chabert, Arancibia Monreal, Bredin, Popelier and Aanesland2012; Lafleur Reference Lafleur2014). Regardless, the following control strategy can be easily generalized if several ion species are considered. With the focus of this work on MEPTs, the current-free condition has been imposed at the thruster outlet (Manente et al. Reference Manente, Trezzolani, Magarotto, Fantino, Selmo, Bellomo, Toson and Pavarin2019). It is worth noting that this condition is truly respected only ‘on average’ since RF wave propagation can induce instantaneous non-null currents. However, this effect is expected to influence mildly the plume dynamics and, in particular, the computation of macroscopic quantities for the propulsive performance. In fact, the RF power deposited outside the source is expected to be only 1 % of the total (Magarotto & Pavarin Reference Magarotto and Pavarin2020). The current-free condition at the thruster outlet can be expressed as

(2.17)\begin{equation} |I_i| = |I_e|. \end{equation}

If quasi-neutrality holds true, from (2.9) and (2.16):

(2.18)\begin{equation} n_i = \frac{\dfrac{|I_i|}{q_iS}}{|E[\widetilde{V_{z,i}}]|} \approx n_e = \frac{\dfrac{|I_e|}{|e|S}}{|E[\widetilde{V_{z,e}}]|}. \end{equation}

Combining (2.17) and (2.18):

(2.19)\begin{equation} |E[\widetilde{V_{z,i}}]| \approx \frac{|e|}{q_i} |E[\widetilde{V_{z,e}}| . \end{equation}

Assuming $\widetilde {V_z} \approx V_z$, from (2.13) and (2.19):

(2.20)\begin{equation} \frac{u_B\sin{\theta_t}}{\theta_t} \approx \frac{|e|}{q_i} |E[\widetilde{V_{z,e}}]|. \end{equation}

The term $|E[\widetilde {V_{z,e}}]|$ being highly nonlinear dependent on the boundary conditions at the thruster outlet (Martinez-Sanchez et al. Reference Martinez-Sanchez, Navarro-Cavallé and Ahedo2015), any noise or deviation from analytical behaviour could lead to a condition that strongly departs from current-free and quasi-neutrality. Consequently, a control loop is needed to reach the steady state while ensuring these conditions (Li et al. Reference Li, Merino, Ahedo and Tang2019).

The loop controls the parameters that mainly drive the simulation, namely the potential at the thruster outlet $\phi _o$, along with the electron and ion currents injected into the domain ($I_e^+$ and $I_i^+$ respectively). From (2.15), $I_i^+$ is assumed constant during the simulation and is given by

(2.21)\begin{equation} I_i^+{=} \frac{q_i\dot{m}}{M}, \end{equation}

with $\dot {m}$ the propellant mass flow rate. To achieve the current-free condition at the steady state, for every time step $t$ the potential at the thruster outlet ($\phi _o^t$) is updated according to the explicit relation

(2.22)\begin{equation} \phi_o^t = \phi_o^{t-1}-k_1\frac{I_e^{t-1}+I_i^{t-1}}{|I_i^{t-1}|}, \end{equation}

with $I_e^{t-1}$ and $I_i^{t-1}$ respectively the electron and ion net currents at the thruster outlet at the previous time step and $k_1$ a positive control coefficient. To achieve quasi-neutrality at the steady state, for every time step the emitted electron current (${I_e^+}^t$) is updated according to the explicit relation

(2.23)\begin{equation} |{I_e^+}^t| = |{I_e^+}^{t-1}|+k_2\frac{n_i}{n_e}(n_i-n_e), \end{equation}

where ion and electron densities are computed at the thruster outlet and $k_2$ is a positive control coefficient. It is worth noting that, according to (2.9), $n_e$ depends on both $I_e^+$ and $I_e^-$, the latter being the absorbed electron current that cannot be imposed directly as an input of the simulation since it is driven by $\phi _o$. Namely (2.22) and (2.23) impose a self-consistent relation between $\phi _o$, $I_e^+$ and $I_i^+$ to enforce current-free and quasi-neutrality conditions.

At the first time step, the quantities $\phi _o$ and $I_e^+$ are imposed according to the analytical values derived in Appendix A for a non-magnetized case, and then updated according to (2.22) and (2.23) up to the achievement of the steady-state condition (see § 2.5). As a result, the simulation provides as an output also the self-consistent value of the total potential drop across the plume (i.e. $\phi _o$) which depends on the geometry of the system and the imposed flow of propellant.

2.5. Time loop termination

A simulation is assumed to be at the steady state if the number of particles leaving through the external boundaries matches the number of new particles introduced into the system by the source. This condition is met for a sufficiently high number of time steps ($\alpha$) in order to avoid a ‘false positive’ induced by PIC noise. Specifically, the steady-state condition is implemented by requiring that (Gallina et al. Reference Gallina, Magarotto, Manente and Pavarin2019)

(2.24)\begin{equation} \frac{1}{\alpha+1}\sum_{t = h-\alpha}^{h}\frac{|N^t_{\textrm{sp}}-N^{t-1}_{\textrm{sp}}|}{N^t_{\textrm{sp}}}<\nu, \end{equation}

where $\nu$ is a user-defined number that represents a steady-state condition (usually $\nu \sim 10^{-4}$), $N^t_{\textrm {sp}}$ is the number of super-particles in the system at time step $t$ and $h$ is the current time step. As a rule of thumb, $\alpha \sim 50$ has been assumed.

2.6. Performance indicators

The thrust $F$ produced by the thruster is obtained by summing, for each species $p$, the contribution due to the momentum exchange and to the pressure. Specifically it is computed with (Zhou et al. Reference Zhou, Pérez-Grande, Fajardo and Ahedo2019)

(2.25)\begin{equation} F =\int_{S_b} \sum_{p}(n_p m_p v_{p,z} \boldsymbol{v_p} \boldsymbol{\cdot} \boldsymbol{k} + n_p k_B T_p \boldsymbol{\hat{z}} \boldsymbol{\cdot} \boldsymbol{k})\,\textrm{d}S_b , \end{equation}

with $S_b$ the external boundary surface, $\boldsymbol {k}$ the external normal to the surface and $\boldsymbol {\hat {z}}$ the unity vector parallel to the $z$ axis. The specific impulse is obtained by

(2.26)\begin{equation} I_{\textrm{sp}} = \frac{F }{g_0 \dot{m}}, \end{equation}

with $g_0$ the standard gravity.