Numerical simulation of an expanding magnetic field plasma thruster: a comparative study for argon, xenon and iodine fuel gases

Abstract

During a space mission, switching to an electric propulsion system from chemical propulsion, once the spacecraft is out of the Earth's gravity, significantly reduces the mission's overall cost. In electric propulsion, the Hall thruster and gridded ion thruster are established technologies. These thrusters compromise mission longevity due to continuous erosion of the device electrode material. To overcome this issue, an electrode less expanding magnetic field plasma thruster or helicon plasma thruster (HPT), was proposed and research is on going worldwide. The HPT shows scaling of thrust with input power while Hall thrusters and ion thrusters do not. Typically, an inert xenon gas is used as a fuel in HPT devices due to a low ionization potential and non-hazardous nature. Xenon is not easily available in nature and during a space mission it needs to be stored in high pressure tanks. Recently, iodine has been proposed as an alternate to xenon as it is easily available and does not have any storage issues. In most of the numerical simulations, argon is used as a fuel gas to reduce the simulation cost. Using a 1D3V particle-in-cell Monte Carlo collision code, we present here a net thrust generation for different fuel gases such as argon, xenon and iodine. We compare plasma flow rates and directed ion beam velocity for different fuel gases having identical inputs. Thrust and plasma flow are investigated for different magnetic field gradients in the plasma expansion region for unidirectional and bidirectional HPT and is reported here. Using iodine fuel, a significant increase in net thrust is obtained for higher magnetic field divergence for identical simulation input parameters while comparing with xenon fueled cases.

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 2008). 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 2008). 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 2003). 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. 2013b). 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 2004). Xenon is not easily available in nature and so it is extracted from liquefied air, which makes it more costly (Bailey 2001). 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 1970). 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. 2011, 2015; Mazouffre 2016; Martinez Martinez, Rafalskyi & Aanesland 2018; Lucken et al. 2019; Manente et al. 2019; Tirila, Demairé & Ryan 2023). Iodine is easily available in nature and is found both in atomic and molecular forms (Mazouffre 2016). 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. (2011) 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. (2016) 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. (2021) demonstrated a 0.6 mN thrust and 600 sec specific impulse for REGULAS, a complete propulsion device using the iodine fuel. Rafalskyi et al. (2021) 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 2005; Baalrud et al. 2011; Takahashi et al. 2018a; Saini & Ganesh 2020, 2022).

A 1D3V particle-in-cell Monte Carlo collision (PIC-MCC) (Hockney & Eastwood 1981; Chen 1990; Jackson 1991; Birdsall & Langdon 1992; Cartwright, Verboncoeur & Birdsall 2000) code, called EPPIC (Saini et al. 2018; Saini & Ganesh 2020, 2022), 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 1970; Qin et al. 2013; Ebersohn, Sheehan & Gallimore 2015). 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. 2015).

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. 2015; Saini & Ganesh 2020), 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. (2015)

(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 (2022).

To energize the stray electrons and set up the plasma discharge, an asymmetric inductive excitation scheme has been used (Meige et al. 2005; Ebersohn et al. 2015; Saini & Ganesh 2020). 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 1992; Verboncoeur et al. 1993; Vahedi & Surendra 1995; Frignani & Grasso 2006; Hamilton 2015; Carbone et al. 2021). 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 1992; Frignani & Grasso 2006; Hamilton 2015; Carbone et al. 2021). 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. 2018a). 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. 2011; Ebersohn et al. 2015).

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 2015; Takahashi, Takao & Ando 2016; Takahashi 2021).

2. Numerical simulation

To start with, in the presence of background neutrals of a particular fuel type, 1 ev seed electrons and ions following the Maxwellian velocity distribution are uniformly distributed spatially. These seed electrons are further energized using an externally applied RF electric field that evolves in time and the simulation attains a steady state after a sufficiently long time. We have compared net plasma thrust and charged particle flow profiles for different fuel gases such as argon, xenon and iodine for identical input parameters. Here we have simulated both a unidirectional and bidirectional expanding magnetic field plasma thruster. In the unidirectional case, the left half of the simulation domain is the heating region and the right half is the particle expansion region. For the bidirectional system, the heating region is situated in the centre and the plasma expands on both sides of the heating region.

Throughout the numerical simulation, for the unidirectional plasma expansion case, the right wall is grounded, i.e. the metal boundary condition (Dirichlet boundary condition), and the left boundary wall floats at a particular potential depending upon the amount of charge flux reaching the wall. For the bidirectional plasma expansion case, both ends of the simulation domain are grounded and particles crossing the simulation boundaries are removed.

To further justify the simulation boundary conditions, one takes into account that the typical length scale of a laboratory experiment is many times larger than the usual length scale of numerical simulations. Hence, performing numerical simulations for such large length scales becomes costly in terms of computation. However, we recently showed that bulk plasma physics and plasma sheath physics are clearly distinguished as the simulation length scale increases (Saini & Ganesh 2020). In the present work we have chosen the simulation length sufficiently large to delineate both the bulk plasma effect and sheath effect. In a HPT one of the major reasons for ion acceleration is the double layer potential formed in the bulk plasma region, which is far apart from the sheath potential in our simulation. Boundary conditions on both sides do affect the ion beam dynamics (Chen 2006; Takahashi, Shida & Fujiwara 2010; Boswell et al. 2015; Takahashi et al. 2018b, 2020); however, we have chosen a simulation length sufficiently large to minimize the sheath electron distribution effect over the bulk plasma. To reduce the computational cost, this approach is often used in numerical simulations (Meige et al. 2005; Ebersohn et al. 2015).

The potential and field quantities are time averaged over the last $10$ RF cycles. The detailed simulation input parameters are given in table 1 (Meige et al. 2005; Saini & Ganesh 2020).

Table 1. Standard simulation input parameters for an expanding magnetic field plasma thruster.

Figure 1 shows the schematic diagram of an expanding magnetic field plasma thruster. In the plasma heating chamber (inside red colour box), RF coils are wrapped around the cylinder to energize the electrons and generate an electron–ion pair via ionization collisions. Magnetic field coils are placed such that the magnetic field is constant inside the heating chamber and rapidly diverse outside the heating region. In a laboratory experiment, to control the net thrust, the fuel gas inflow rate and magnetic field gradient are the two main control parameters. Here our neutral gas pressure is constant and the only parameter to control the net thrust is the magnetic field gradient in the plasma expansion region.

Figure 1. Expanding the magnetic field plasma thruster schematic diagram.

Figures 2(a) and 2(b) show the magnetic field spatial profile for both a unidirectional and bidirectional expanding magnetic plasma thruster, respectively. As we have already discussed, the magnetic field gradient in the plasma expansion region increases with decreasing $lb$ values. For the bidirectional case, we can simultaneously control magnetic field expansion rates on both sides of the plasma heating region.

Figure 2. Spatial magnetic field profile for a (a) unidirectional and (b) bidirectional expanding plasma case with different $lb$ values. A steady magnetic field exists in the plasma heating region, while the magnetic field diverges in the plasma expansion region.

2.1. Unidirectional expanding magnetic field plasma thruster

For the unidirectional plasma expansion case, we performed simulations for different fuel gases such as argon, xenon and iodine for identical input parameters. Figure 3 shows the spatial plasma potential profile for the unidirectional plasma expansion case. The plasma potential amplitude in the heating region and gradient in the plasma expansion region increase with decreasing $lb$. For lower magnetic field gradient values, electrons have enough time to follow the ions in the plasma expansion region, which overall reduce the potential amplitude in the heating region and gradient in the expansion region. For the argon fuel case, the plasma potential shows a fall of around 20 V and the left wall charges around 7 V that helps to contain the plasma in the simulation region. For xenon and iodine, both the plasma potential amplitude in the heating region and the potential gradient in the expansion region decrease as compared with the argon case. The atomic mass of xenon and iodine is higher than argon, which causes the lower plasma potential amplitude in the heating region and the lower potential fall in the expansion region. For xenon and iodine, the maximum amplitude for the $lb=0.01$ case is around 20 V and the left wall attains a potential around 3 V.

Figure 3. Plasma potential spatial profile with different $lb$ values for argon, xenon and iodine, where the fuel neutral pressure is fixed at 1.23 mTorr. Potential profiles for (a) argon (b) xenon and (c) iodine.

We have shown the ion phase space plots for argon, xenon and iodine in figure 4 for the $lb=0.01$ case. The phase space plots show a clear signature of the directed ion beam in the plasma expansion region. The argon ion phase space plot shows the beam velocity around 10 000 m sec$^{-1}$ in the expansion region. At both ends of the simulation walls, ions are accelerated due to the plasma sheath. Xenon and iodine ion phase space plots show directed ion beam velocities around 4000 m sec$^{-1}$. For argon, the potential gradient in the plasma expansion region is higher compared with xenon and iodine, which causes a high electric field in the bulk plasma and results in a larger directed ion beam amplitude.

Figure 4. At the end of simulation time, the ion phase space plot for (a) argon, (b) xenon and (c) iodine, where the neutral fuel gas pressure is fixed at 1.23 mTorr. Here we have shown the $lb=0.01$ case.

Figure 5 shows the ion velocity distribution function for different magnetic field gradients in the plasma expansion region. Here we have sampled the ion velocity distribution around $z=0.07$ m. For the argon case, the maximum directed beam velocity is around 9000 m sec$^{-1}$ for $lb=0.01$. Apart from a clear signature of the directed ion beam, it is observed that ion-neutral elastic and charge exchange collisions cause a broadened ion velocity distribution function. As we increase $lb$, the directed ion beam velocity decreases in the plasma expansion region and, for $lb=0.04$, the distribution function shows a shifted Maxwellian kind behaviour. For the xenon and iodine plasma case, the maximum ion beam velocity is around 4000 m sec$^{-1}$. Here also while increasing the $lb$ value, the directed ion beam velocity amplitude decreases and, for the $lb=0.04$ case, the ion velocity distribution is near Maxwellian.

Figure 5. Ion velocity distribution function at $z=0.07$ m for (a) argon, (b) xenon and (c) iodine, where the fuel gas neutral pressure is fixed at 1.23 mTorr. A directed ion beam having a velocity around 9 km s$^{-1}$ for the argon plasma and around 4 km s$^{-1}$ for the xenon and iodine plasma cases are observed.

Figure 6 shows the mean ion velocity having different $lb$ values for argon, xenon and iodine. As we previously explained for the ion velocity distribution function, the mean ion velocity also follows a similar signature that is quite obvious. In the beginning of plasma expansion region, a sharp bulk electric field (plasma double layer) occurs that accelerates the ions and turns into the beam. The strength of the plasma double layer is low for the xenon and iodine cases, which causes a lower mean ion velocity compared with the argon plasma case.

Figure 6. Spatial ion mean velocity profile for (a) argon, (b) xenon and (c) iodine, where the fuel gas neutral pressure is fixed at 1.23 mTorr. Ions get accelerated due to the plasma double layer electric field at the start of the expansion region.

In figure 7 we show the ion flow rate profiles for argon, xenon and iodine plasma expansion cases. Here the ion flow rate is defined as $n_{i}v_{i}A$, where $n_{i}$ is the ion spatial density, $v_{i}$ is the ion mean velocity and $A$ is the simulation cross-sectional area. In the plasma heating region, the ion flow profile shows a constant gradient due to a source term in the plasma steady state. The ion flow rate is nearly constant in the plasma expansion region because this region does not have any plasma source term. For the xenon and iodine plasmas, the ion flow rates are smaller compared with the argon plasma case due to a smaller mean ion velocity.

Figure 7. Ion flow profile for (a) argon, (b) xenon and (c) iodine for different $lb$ values, where the neutral fuel gas pressure is fixed at 1.23 mTorr.

Figure 8 show the net thrust generation for argon, xenon and iodine having different magnetic field divergence rates in the plasma expansion region. The net ion thrust has been measured at axial location $z=0.07$ m for the magnetic nozzle configuration. The thrust measurement diagnostic is $T= \dot m v_{\rm exh}$, where $\dot m$ (ion flow rate times fuel atomic mass) represents the ion mass flow rate and $v_{\rm exh}$ (see figure 5) refers to the ion exhaust velocity. In an expanding magnetic plasma thruster, the radial electron pressure is converted in an axial momentum flux and the electron temperature is reduced in the axial direction by losing their internal energy while ions get accelerated. The plasma is adiabatic in the diverging section of a magnetic nozzle, so any energy lost by the electrons must be transferred to the ions via the electric field. The total force over a bounded plasma system of the magnetic nozzle is the sum of the static electron pressure force exerted on the axial boundary, the Lorentz force onto a magnetic nozzle (para axial approximation for the one-dimensional case) and the force caused by the bulk electric field that arises due to sudden plasma expansion (Fruchtman 2006; Takahashi et al. 2011; Takahashi, Charles & Boswell 2013a; Takahashi 2019; Takahashi, Sugawara, & Ando 2020). The $lb=0.01$ case, which has the largest magnetic field gradient, shows the maximum thrust generation that reduces with increasing $lb$ values. Here, low $lb$ values show significant change in the net thrust generation for different fuel gases. Iodine shows a $1.5$ times jump in the thrust compared with argon for identical input parameters, which is basically due to a high atomic mass.

Figure 8. Net plasma thrust profile with different magnetic field divergence for argon, xenon and iodine fuels.

2.2. Bidirectional expanding magnetic field plasma thruster

We have also investigated plasma dynamics in an expanding magnetic field for a bidirectional plasma thruster configuration. Here we have chosen one specific value of $lb$ throughout the simulation that is 0.02. The simulation length and grid points are doubled here to have a consistency with the unidirectional case. Simulation macro particle numbers are also increased two times such that the number of particles per cell is identical to the unidirectional case. The rest of the input parameters are the same as in table 1.

Figure 9 shows the spatial potential profile for argon, xenon an iodine for fixed $lb=0.02$. The peak potential in the plasma heating region for argon is around $32$ V, which is significantly larger compared with xenon and iodine fuel plasmas. Similarly, this shows a larger potential gradient for the argon plasma compared with the xenon and iodine cases, which turns into a large electric field in the bulk of the plasma. The peak potential in the heating region and potential gradient in the plasma expansion region depends upon both the fuel gas atomic mass and its collisional probability with the neutral fuel gas. Here argon has rapidly expanded compared with xenon and iodine, which result in a larger potential gradient or bulk electric field.

Figure 9. Spatial plasma potential profile for argon, xenon and iodine for the fixed $lb=0.02$ case.

In figure 10 we show instantaneous ion phase space plots for argon, xenon and iodine at $t=t_{\rm end}$. In all three cases, the maximum ion density lies in the plasma source region that is situated at the centre for a bidirectional plasma expansion simulation. A bulk electric field at the start of the expansion region accelerates the ions and finally it forms a directed ion beam. Here near the simulation end walls, ions are accelerated due to the sheath electric field. A directed ion beam having amplitude 9000 m sec$^{-1}$ is observed for the argon plasma expansion and the beam amplitude for xenon and iodine is nearly similar at around 3500 m sec$^{-1}$.

Figure 10. Instantaneous ion phase space plot for (a) argon, (b) xenon and (c) iodine, where the neutral fuel gas pressure is fixed at 1.23 mTorr, for a bidirectional plasma expansion at $t=t_{\rm end}$. Here we show the $lb=0.02$ case.

Figure 11 shows the ion velocity distribution function for argon, xenon and iodine for a fixed $lb=0.02$. This ion distribution is measured at $z=0.14$ m on the simulation axis. Peaks in the velocity distribution is evident the formation of accelerated ion beam which impart the net thrust in the opposite direction. Here argon shows the maximum amplitude for the directed ion beam. These ion distribution functions have broadening that is basically due to ion-neutral charge exchange collisions and elastic collisions.

Figure 11. Ion velocity distribution measured at axial location $z=0.14$ m for a fixed $lb=0.02$ case.

In figure 12 we show the ion mean velocity spatial profile for argon, xenon and iodine for the fixed $lb=0.02$ case. This shows two jumps in the ion mean velocity profile on either side. The first jump in the ion velocity right after the heating region is due to a bulk electric field that occurs due to plasma expansion dynamics and the second time velocities are picking up near the simulation walls, which is due to the sheath electric field. Similarly, argon also shows larger numbers of ion mean velocity compared with xenon and iodine.

Figure 12. Ion mean velocity profile for argon, xenon and iodine for the fixed $lb=0.02$ case.

Figure 13 shows the ion flow profile for argon, xenon and iodine plasma for a fixed $lb=0.02$ case. Ion flow is governed by both the mean velocity profile and the ion density spatial profile. The net plasma thrust is measured using this ion flow profile that is 2.0 mN, 2.67 mN and 3.1 mN for argon, xenon and iodine, respectively. This also follows a similar trend to the unidirectional plasma expansion case and net thrust is increased due to the fuel atomic mass for xenon and iodine.

Figure 13. Ion spatial flow profile for argon, xenon and iodine for the fixed $lb=0.02$ case.

3. Summary and conclusions

We simulated expanding plasma dynamics along the thruster axial direction using argon, xenon and iodine as a fuel gas. Plasma expansion physics are incorporated using the magnetic flux conservation model. We simulated plasma thrust generation for a unidirectional and bidirectional expansion case using different magnetic field divergence rates in the plasma expansion region. For unidirectional plasma expansion, we performed a detailed investigation for plasma dynamics using different $lb$ values. Potential gradient in the plasma expansion region reduces with increasing $lb$ value for argon, xenon and iodine fuel cases. A maximum potential gradient is observed for the $lb=0.01$ case with argon fuel, which results in a larger bulk electric field and causes a high velocity directed ion beam. Clear emergence of a directed beam is observed in the ion phase space for argon, xenon and iodine fuel cases. For the argon plasma case, the directed beam velocity is around 9000 m sec$^{-1}$ and a similar directed beam velocity amplitude is observed for xenon and iodine plasma cases at around 4000 m sec$^{-1}$. We observed maximum thrust for the $lb=0.01$ value and it decreases as we increase the $lb$ values. Due to a high atomic mass number, xenon and iodine show larger thrust numbers. For iodine, maximum thrust is observed at around 3.4 mN. For the same input parameters, the thrust observed for argon and xenon for the $lb=0.01$ case is $2.2$ mN and $2.7$ mN, respectively. For larger $lb$ values, the net thrust variation between different fuel gases reduces and the $lb=0.04$ case shows similar thrust numbers.

We also simulated a bidirectional plasma expansion dynamics for different fuel gases such as argon, xenon and iodine. Here we presented simulation results only for the $lb=0.02$ case. Similar to the unidirectional case, peak plasma potential and potential gradient in the expansion region are also maximum for the argon plasma. For the bidirectional case, we showed symmetric plasma expansion, which means both sides of the heating region $lb$ values are the same. Here we also obtained maximum thrust for the iodine fuel plasma at axial location $z=0.14$ m, which is $3.1$ mN.

As we already discussed, to capture the plasma expansion dynamics, we used the magnetic flux conservation model that holds good only if the radial component of the magnetic field is very small compared with the axial component, i.e. $B_{r} \ll B_{Z}$. In laboratory conditions, this approximation may not hold. To alleviate this situation, a minimum 2D3V PIC-MCC code is essential. We have developed a cylindrical 2D3V PIC-MCC code that is in the testing phase and a full scale simulation will be presented soon.