1. Introduction

A chemical propulsion engine is used to lift off the spacecraft and once it is out of the Earth's gravity, one switches to the electric propulsion engine due to a low thrust requirement. In a chemical propulsion engine, hot gas exhaust velocities are in the range of 3–4 Km sec$^{-1}$ and a specific impulse is generally less than 500 sec. An electric propulsion device has a specific impulse in the range of 1500–20 000 sec that helps to reduce the overall spacecraft mass for a mission and turns into a mission cost reduction (Goebel & Katz Reference Goebel and Katz2008). A gridded ion plasma thruster and Hall plasma thruster are established technologies and they are successfully deployed in several space missions for orbit correction and deep space missions (Goebel & Katz Reference Goebel and Katz2008). In these devices, continuous erosion of electrode material due to energetic charged particle bombardment, compromises the mission longevity. An electrode-less plasma thruster or helicon plasma thruster (HPT) is an alternative to overcome the issue of electrode material erosion (Charles & Boswell Reference Charles and Boswell2003). Helicon plasma thrusters have good thrust to power scaling. In a HPT, magnetic field coils have significant power consumption and make the overall device bulky. For a space mission, overall cost is directly propositional to total spacecraft mass. Recently, a laboratory experiment for the HPT has been performed with permanent magnets to obtain thrust (Takahashi et al. Reference Takahashi, Charles, Boswell and Ando2013b). This experiment shows a alternate approach for HPT devices to reduce the overall mass. A HPT operation with permanent magnets has several constraints such as one cannot do magnetic field steering and cannot change the magnetic field gradient for some specific purposes.

In a HPT device the most commonly used fuel gas is xenon. It is non-hazardous, non-reactive and stored in high pressure tanks during the operation. For a fuel gas, the first and second ionization energies should be far apart. For xenon, the first ionization energy is 12.12984 eV and the second ionization energy is 20.97503 eV (Saloman Reference Saloman2004). Xenon is not easily available in nature and so it is extracted from liquefied air, which makes it more costly (Bailey Reference Bailey2001). In principle, any atom can be used as a fuel because it can be ionized using required ionization potential energy and it can be accelerated using an external electromagnetic field. In the past, mercury and cesium were used as fuels due to their high mass and low ionization energy, but they are toxic and highly reactive (Brewer Reference Brewer1970). In recent times, iodine is considered as a very promising candidate and a good alternate of xenon for HPT fueling and several numerical simulation and laboratory experiments have been performed worldwide (Szabo et al. Reference Szabo, Pote, Paintal and Robin2011, Reference Szabo, Robin, Paintal, Pote, Hruby and Freeman2015; Mazouffre Reference Mazouffre2016; Martinez Martinez, Rafalskyi & Aanesland Reference Martinez Martinez, Rafalskyi and Aanesland2018; Lucken et al. Reference Lucken, Marmuse, Bourdon, Chabert and Tavant2019; Manente et al. Reference Manente, Trezzolani, Magarotto, Fantino, Selmo, Bellomo, Toson and Pavarin2019; Tirila, Demairé & Ryan Reference Tirila, Demairé and Ryan2023). Iodine is easily available in nature and is found both in atomic and molecular forms (Mazouffre Reference Mazouffre2016). The dissociation energy of molecular iodine $I_{2}$ is small (1.567 eV), which helps to form the atomic $I$ plasma easily. Iodine does not have any storage issue and during HPT operation, solid iodine is vapourised by sublimation and used as a fuel. As iodine is known to be relatively reactive, it is expected that proper precautions are taken during operation. Szabo et al. (Reference Szabo, Pote, Paintal and Robin2011) experimentally demonstrated that the Hall effect thruster with iodine vapour as a fuel showed similar performance to xenon fuel. These authors measured thrust, efficiency and specific impulse using both xenon and iodine fuel. Grondein et al. (Reference Grondein, Lafleur, Chabert and Aanesland2016) used iodine as fuel instead of xenon for a gridded ion thruster and showed that the overall performance is similar for both fuels. They showed that below a certain propellant mass flow rate (1.3 mg sec$^{-1}$), iodine showed higher efficiency than xenon. Bellomo et al. (Reference Bellomo, Magarotto, Manente and Trezzolani2021) demonstrated a 0.6 mN thrust and 600 sec specific impulse for REGULAS, a complete propulsion device using the iodine fuel. Rafalskyi et al. (Reference Rafalskyi, Martínez and Habl2021) reported an in-orbit demonstration of an iodine plasma thruster. They showed that by using an iodine fueled plasma thruster one could achieve more efficiency in comparison to a xenon plasm thruster. They performed small satellite maneuvering and confirmed using satellite tracking system. However, in numerical simulation and laboratory experiments argon is used most commonly as a fuel to cut down simulation cost and experiment cost, respectively.

In the current work, in our numerical simulations, we used argon, xenon and iodine as a fuel gas and compared the net thrust generation. We compare here, net particle flow and directed beam velocity in the plasma expansion region for different gases with different magnetic field gradients in the plasma expansion region.

We perform this study for both a unidirectional and bidirectional expanding magnetic field plasma thruster with different fuel gases. It is now known that a unidirectional plasma thruster may be used for satellite acceleration in a particular direction and a bidirectional plasma thruster can be used to accelerate/decelerate a space satellite or for debris removal via switching mode of operation (Meige, Boswell & Charles Reference Meige, Boswell and Charles2005; Baalrud et al. Reference Baalrud, Lafleur, Boswell and Charles2011; Takahashi et al. Reference Takahashi, Charles, Boswell and Ando2018a; Saini & Ganesh Reference Saini and Ganesh2020, Reference Saini and Ganesh2022).

A 1D3V particle-in-cell Monte Carlo collision (PIC-MCC) (Hockney & Eastwood Reference Hockney and Eastwood1981; Chen Reference Chen1990; Jackson Reference Jackson1991; Birdsall & Langdon Reference Birdsall and Langdon1992; Cartwright, Verboncoeur & Birdsall Reference Cartwright, Verboncoeur and Birdsall2000) code, called EPPIC (Saini et al. Reference Saini, Pandey, Trivedi and Ganesh2018; Saini & Ganesh Reference Saini and Ganesh2020, Reference Saini and Ganesh2022), is well tested and used here to study the dynamics of an expanding plasma following the diverging magnetic field lines such as a HPT or magnetic nozzle. This code resolves the magnetic nozzle (HPT axial direction) axial direction and all three velocity components. The geometrical forces (due to device geometry) acting on charged particles are considered and the Boris particle pusher algorithm has been used to solve the charged particle forces (Boris Reference Boris1970; Qin et al. Reference Qin, Zhang, Xiao, Liu, Sun and Tang2013; Ebersohn, Sheehan & Gallimore Reference Ebersohn, Sheehan and Gallimore2015). Throughout the simulation, we have considered that electrons are magnetized and are confined to a magnetic flux tube. A magnetic flux conservation model is used to incorporate the charge particle expansion physics. In this model the spatially varying magnetic flux tube cross-sectional area is used to include the magnetic field expansion effect on charged particles while following the diverging magnetic field lines (Ebersohn et al. Reference Ebersohn, Sheehan and Gallimore2015).

In a HPT device or magnetic nozzle, while plasma expands from the source region (heating region) to the expansion region, charged particle density typically varies 1–2 orders when comparing the source region with the plasma expansion region. This change in charged particle density along the axial direction is mimicked by variation of the magnetic flux tube cross-sectional area. For a known magnetic field axial profile and flux tube inlet, one can calculate the cross-sectional area throughout the nozzle axis, i.e.

(1.1)\begin{equation} A_{1}B_{1}=A_{2}B_{2}, \end{equation}

where $A_{1}$ and $A_{2}$ are magnetic flux tube cross-sectional areas. Here, one can easily calculate $A_{2}$ for an already given magnetic field axial profile. The variation of the magnetic flux tube cross-section area enables the Poisson equation to include the magnetic field expansion effect on charged particles while following the diverging filed line.

The magnetic nozzle is a cylindrical device and coordinate forces are incorporated here (Ebersohn et al. Reference Ebersohn, Sheehan and Gallimore2015; Saini & Ganesh Reference Saini and Ganesh2020), i.e.

(1.2)\begin{equation} \frac{{\rm d} \boldsymbol{v}}{{\rm d} t}=\frac{q}{m}(\boldsymbol{E}+\boldsymbol{v}{\times \boldsymbol{B}})+\frac{v^2_{\theta}}{r_{L}}\hat{r}- \frac{v_{\theta}v_{r}}{r_{L}}\hat{\theta}, \end{equation}

where $\boldsymbol {v}=(v_{r},v_{\theta },v_{z})$.

We have assumed that charged particles are magnetized and they are displaced over one Larmor radius $r_{L}$ from magnetic field lines. Here, the axial magnetic field $B_{z}$ does not change over the particle orbit. Using the above approximations, Gauss's law of magnetism, $\triangledown \cdot B =0$, can be written in a simplified form:

(1.3)\begin{equation} B_{r}={-}\frac{r_{L}}{2} \frac{{\rm d} B_{z}}{{\rm d} z}. \end{equation}

After substituting $B_{r}$ into (1.2), the simplified model equations are

(1.4)\begin{gather} \frac{{\rm d} v_{z}}{{\rm d} t}={-}\frac{1}{2B_{z}} \frac{{\rm d} B_{z}}{{\rm d} z} v^{2}_{\theta}, \end{gather}

(1.5a,b)\begin{gather}\frac{{\rm d} v_{\theta}}{{\rm d} t}=\frac{1}{2B_{z}} \frac{{\rm d} B_{z}}{{\rm d} z} v_{\theta} v_{z},\quad \frac{{\rm d} v_{r}}{{\rm d} t}=0. \end{gather}

Throughout the simulations, force equations are solved in the frame of magnetic field lines. The force equations following the magnetic field frame of reference are

(1.6)\begin{gather} \frac{{\rm d} v_{{\parallel}}}{{\rm d} t}={-}\frac{1}{2B} \frac{{\rm d} B}{{\rm d} s} v^{2}_{{\perp}}, \end{gather}

(1.7)\begin{gather}\frac{{\rm d} v_{{\perp}}}{{\rm d} t}=\frac{1}{2B} \frac{{\rm d} B}{{\rm d} s} v_{{\perp}}v_{{\parallel}}, \end{gather}

where $s$ is defined in the direction along the magnetic field line and $v_{\perp } =\sqrt {v^2_{\theta }+v^2_{r}}$.

To solve the force equations, the Boris particle pusher algorithm has been used along with energy conservation to solve the magnetic field forces. Using the energy conservation principle for a given magnetic field, one can write

(1.8)\begin{equation} (v^{n}_\parallel)^2+(v^{n}_{{\perp}})^2=(v^{n+1}_\parallel)^2+(v^{n+1}_{{\perp}})^2=K. \end{equation}

After a small algebraic calculation, the parallel velocity component at the ${n+1}$th time step is Ebersohn et al. (Reference Ebersohn, Sheehan and Gallimore2015)

(1.9)\begin{equation} v^{n+1}_{{\parallel}}=\sqrt{K} \tanh\left[\tanh^{{-}1}\left(\frac{v^{n}_{{\parallel}}}{\sqrt{K}}\right)- \zeta\sqrt{K}\Delta t\right], \end{equation}

where $\zeta =({1}/{2B}) ({{\rm d}B}/{{\rm d}s})$ and $K$ are constant over $\Delta t$. One can find a detailed algebraic calculation in Saini & Ganesh (Reference Saini and Ganesh2022).

To energize the stray electrons and set up the plasma discharge, an asymmetric inductive excitation scheme has been used (Meige et al. Reference Meige, Boswell and Charles2005; Ebersohn et al. Reference Ebersohn, Sheehan and Gallimore2015; Saini & Ganesh Reference Saini and Ganesh2020). A time varying radio frequency electric field $E_{y}$ is applied normal to the axial direction to energize the plasma electrons. Here the $\hat {y}$ direction belongs to the $\hat {r}-\hat {\theta }$ plane. Using Maxwell's equation, the total current density is

(1.10)\begin{gather} J_{\rm total} = J_{\rm displacement} + J_{\rm conduction}, \end{gather}

(1.11)\begin{gather}J_{\rm total}=J_{0}\sin(\omega t)=\epsilon_{0} \frac{\partial E_{y}}{\partial t}+e\varGamma_{e,y}, \end{gather}

where $\varGamma _{e,y}$ is electron flux along the $y$ axis. Here $\varGamma _{e,y}$ is

(1.12)\begin{equation} \varGamma_{e,y}=\frac{N_{\sigma}}{L_{\rm source}}\sum_{i \in {\rm source}} v_{e,y,i}, \end{equation}

using Maxwell's equations, which results in

(1.13)\begin{equation} \frac{\partial E_{y}}{\partial t}=\frac{1}{\epsilon_{0}}\left[ J_{0}\sin(\omega t)-\frac{eN_{\sigma}}{L_{\rm source}}\sum_{i \in {\rm source}} v_{e,y,i} \right], \end{equation}

where $J_{0}$ refers to the amplitude of the radio frequency (RF) current density, $L_{source}$ is the plasma source region length and $N_{\sigma }$ is the simulation particle weight.

To set a plasma discharge, stray electrons are accelerated by an external RF electric field, creating a electron–ion pair via ionization that evolves further in time until the system reaches a steady state (Sakabe & Izawa Reference Sakabe and Izawa1992; Verboncoeur et al. Reference Verboncoeur, Alves, Vahedi and Birdsall1993; Vahedi & Surendra Reference Vahedi and Surendra1995; Frignani & Grasso Reference Frignani and Grasso2006; Hamilton Reference Hamilton2015; Carbone et al. Reference Carbone, Graef, Hagelaar, Boer, Hopkins, Stephens, Yee, Pancheshnyi, van Dijk and Pitchford2021). Once the system attains a steady state, several diagnostics are used to understand the plasma dynamics. As already discussed, three different fuel gases (argon, xenon and iodine) are used. Here we have considered elastic, excitation and ionization collisions for neutral-electron interaction and for neutral-ion interaction, elastic and charge exchange collisions are also taken into account (Sakabe & Izawa Reference Sakabe and Izawa1992; Frignani & Grasso Reference Frignani and Grasso2006; Hamilton Reference Hamilton2015; Carbone et al. Reference Carbone, Graef, Hagelaar, Boer, Hopkins, Stephens, Yee, Pancheshnyi, van Dijk and Pitchford2021). In the present study, only atomic iodine chemistry and its dynamics is considered in net plasma thrust generation.

Here both a unidirectional and bidirectional expanding magnetic field plasma thruster are simulated. The spatial magnetic field profile for a unidirectional plasma thruster is

(1.14)\begin{equation} B(z) = \begin{cases} B_{0} & \mbox{if } z\leq 0.5 \times L, \\ \dfrac{B_{0}}{\left(1+\dfrac{(z-z_{0})^2}{lb^{2}}\right)^{3/2}} & \mbox{if } z\geq 0.5\times L. \end{cases} \end{equation}

The magnetic field profile for a bidirectional expanding magnetic field plasma thruster is

(1.15)\begin{equation} B(z) = \begin{cases} B_{0} & \mbox{if } 0.4\times L \leq z\leq 0.6 \times L, \\ \dfrac{B_{0}}{\left(1+\dfrac{(z-z_{0})^2}{lb^{2}}\right)^{3/2}} & \mbox{if } 0.6\times L < z < 0.4\times L, \end{cases} \end{equation}

where $L$ represents the total simulation length, $B_{0}$ is the magnetic field strength and $lb$ is a control parameter for the gradient length scale of the magnetic field. In a realistic experimental device this length scale can be controlled using the magnetic field coil current and its geometrical location (Takahashi et al. Reference Takahashi, Charles, Boswell and Ando2018a). In simulations, $lb$ is a control parameter to include the information of how fast magnetic field lines spatially diverge in an expanding magnetic field plasma thruster. For the unidirectional case, we have chosen $z_{0}=0.5\times L$ and, for the bidirectional case, $z_{0}=0.6\times L$ is chosen on the right-hand side and $z_{0}=0.4\times L$ is chosen on the left-hand side of the plasma source region. The magnetic field divergence in the plasma expansion region increases with decreasing $lb$ values and the $lb=0.01$ case has the largest magnetic field divergence rate (Baalrud et al. Reference Baalrud, Lafleur, Boswell and Charles2011; Ebersohn et al. Reference Ebersohn, Sheehan and Gallimore2015).

In this model, neutral flow dynamics is not included, although, it is a critical parameter in a working device for net thrust generation (Takao & Takahashi Reference Takao and Takahashi2015; Takahashi, Takao & Ando Reference Takahashi, Takao and Ando2016; Takahashi Reference Takahashi2021).