1 Introduction

In recent times, human exploitation of space has enormously increased, resulting in an increased number of flying pieces of debris or space junk which can collide with space stations, space shuttles and live satellites and may result in heavy damage. This increase in the amount with time of natural and man-made debris has posed a serious challenge to the safety of human activities in space, as well as for smooth operation of man-made satellites (Laurance & Brownlee Reference Laurance and Brownlee1986; Peter Reference Peter1988; Liou & Johnson Reference Liou and Johnson2006; Dov Reference Dov2020). A space debris orbits around Earth with an angular momentum such that the centrifugal force balances the Earth's gravitational force and maintains a constant altitude. For safe human activity in space, de-orbiting space junk is an emergent technological challenge. Due to the Earth's gravitational force, the natural orbit decay rate of a space debris is very small and, left to itself, space debris will take an enormous amount of time to re-enter into the Earth's atmosphere depending upon the current altitude of the space debris. An external drag force acting on space debris can decelerate/accelerate the debris and can result in the re-entry of debris into the Earth's atmosphere where it can burn out (Shan, Guo & Gill Reference Shan, Guo and Gill2015).

Several techniques have been proposed in the past to manoeuvre the re-entry of space debris into the Earth's atmosphere. For example, a ground based powerful laser pulse system can be used to ablate debris material. This will change the debris orbital velocity and make it possible for re-entry into the Earth's atmosphere where burn out can take place (Rubenchik, Fedoruk & Turitsyn Reference Rubenchik, Fedoruk and Turitsyn2014). A low divergence ion beam can be used for contactless removal of space junk by direct momentum transfer (Bombardelli & Peláez Reference Bombardelli and Peláez2011; Kitamura, Hayakawa & Kawamoto Reference Kitamura, Hayakawa and Kawamoto2014). Robotic arms and tether nets are potential methods to capture active space debris and de-orbit them (Castronuovo Reference Castronuovo2011; Shan, Guo & Gill Reference Shan, Guo and Gill2017). An electro-dynamic tether technology is investigated as it has an efficient orbital transfer system (Nishida et al. Reference Nishida, Kawamoto, Okawa, Terui and Kitamura2009).

More recently, a magnetic nozzle plasma thruster with both ends open or a bi-directional plasma thruster for active removal of space debris has been suggested by Takahashi et al. (Reference Takahashi, Charles, Boswell and Ando2018b) and Takahashi (Reference Takahashi2019). Figure 1(a) shows the concept of space debris removal via direct momentum transfer of an ejected beam, using a bi-directional plasma thruster. Here, $F_{1}$ and $F_{2}$ are the momentum fluxes ejected from two opposite ends of the device. A drag force $F_{1}$ decelerates the space debris and helps it to re-enter into the Earth's atmosphere where it is burnt out (Goebel & Katz Reference Goebel and Katz2008). Figure 1(b) shows a magnetic nozzle with both ends open, forming a single thruster device where the momentum flux is controlled via a magnet coil current and neutral fuel flow rate. This electric propulsion device can be mounted on a satellite and can act as a debris remover or accelerate–decelerate the satellite by altering the magnetic field profile using the coil current. Takahashi et al. experimentally demonstrated that, for a symmetric magnetic field expansion profile at both ends of the radio frequency (RF) heating region, a net thrust is imparted at a plate suspended into the experimental chamber. Here, the suspended plate works as a dummy space debris. A net zero thrust is imparted on a flying satellite working with a symmetric magnetic field expansion bi-directional plasma thruster; i.e. $F_{1}=F_{2}$. For an asymmetric magnetic field expansion at both ends of the heating region, a net thrust at a satellite can be imparted in a particular direction depending upon the magnetic field profile; i.e. $F_{1}\neq F_{2}$. Hence, a bi-directional plasma thruster can have three working modes viz., a satellite acceleration mode, a deceleration mode and a debris removal mode by altering the magnet coil current.

Figure 1.

(a) Working principle of a bi-directional plasma thruster while operating in space debris remover mode. (b) A single bi-directional plasma thruster where the momentum fluxes ejected at both ends are controlled by external operating nobes. (Picture credit: Kazunori Takahashi, Christine Charles, Rod W. Boswell & Akira Ando, Demonstrating a new technology for space debris removal using a bi-directional plasma thruster, Scientific Reports 8:14417 (2018), with approval of K. Takahashi.).

In this work, we have performed a set of numerical experiments for a doubly open-ended magnetic nozzle plasma thruster or a bi-directional plasma thruster using a standard particle in cell (PIC) simulation technique. Although this is an inherently three-dimensional problem, we have performed our simulations in one dimension, along the length of the nozzle's symmetric axis, as shown in figure 2. While averaging over two spatial directions (i.e. radial and angular directions), several physically relevant effects (e.g. $E \times B$ dynamics, radial transport; E, electric field and B, magnetic field) would naturally be neglected. To partially compensate, we have incorporated the effect of magnetic field expansion and compression on the plasma via a magnetic flux conservation model (Ebersohn, Sheehan & Gallimore Reference Ebersohn, Sheehan and Gallimore2015). We believe our reduced model can be useful in the understanding of experimental results for a bi-directional plasma thruster and to search the vast parameter space during optimization. In this work, variations in the directed ion beam velocity amplitude are scanned with a magnetic field expansion rate and neutral Ar gas density. A net thrust is generated for an asymmetric magnetic field configuration and qualitatively compared with recent experimental results for a bi-directional plasma thruster. To capture a more detailed picture of a bi-directional plasma thruster, a higher spatial dimensional code (e.g. 2D3V) may be required which can be expensive in terms of computational resources and time.

Figure 2.

Schematic diagram of a two-open-ended magnetic nozzle where our simulation domain is the nozzle axis, which is shown by the red coloured line.

A 1D3V PIC 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 (Expanding Plasma PIC Solver) (Saini et al. Reference Saini, Pandey, Trivedi and Ganesh2018; Saini & Ganesh Reference Saini and Ganesh2020), is well benchmarked and used here to study the plasma dynamics in an expanding magnetic field plasma device such as the helicon plasma thruster (HPT) or a magnetic nozzle (Saini et al. Reference Saini, Pandey, Trivedi and Ganesh2018; Saini & Ganesh Reference Saini and Ganesh2020). The code resolves one spatial dimension, which is the direction of the nozzle axis and all three velocity dimensions. The geometrical forces acting on charged particles are incorporated and solved using the Boris particle pusher algorithm (Boris Reference Boris1970; Qin et al. Reference Qin, Zhang, Xiao, Liu, Sun and Tang2013; Ebersohn et al. Reference Ebersohn, Sheehan and Gallimore2015). We have assumed implicitly that electrons are magnetized and trapped inside a magnetic flux tube. A flux tube cross-sectional area (magnetic flux conservation model) is used to incorporate the effect of magnetic field compression and expansion on charged particles (Ebersohn et al. Reference Ebersohn, Sheehan and Gallimore2015).

In a HPT or magnetic nozzle, plasma expands from source region into an expansion region following strong guided magnetic field lines. Typically, in the plasma expansion region, particle density varies 1–2 orders in comparison with the source region. The change in density can be mimicked by a variation of the cross-sectional area of the magnetic flux tube. As plasma expands into the diverging magnetic field, particles following the magnetic field lines or plasma will get expanded to follow the magnetic flux surfaces. The cross-sectional area of the magnetic flux tube can be calculated if one uses the approximation that the radial contribution of the flux is negligible in comparison with the total flux. For known flux tube inlet cross-sectional area and magnetic field strength, the cross-sectional area can be calculated throughout the domain at any location along the $z$-axis, i.e.

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

where $A_{1}$ and $A_{2}$ are the cross-sectional areas of the magnetic flux tube at different axial locations. Here, $A_{2}$ can be calculated if $B_{1}$, $B_{2}$ and $A_{1}$ are known, where $B_{1}$ and $B_{2}$ are magnetic field values at two different axial locations.

This cross-sectional area is then used in the Poisson equation to incorporate the axial density variation effect, which in turn helps to mimic the magnetic field compression and expansion effect on plasma particles.

A magnetic nozzle device or HPT is a cylindrical device. Inertial effects due to the coordinate system can be incorporated into the force equation as (Ebersohn et al. Reference Ebersohn, Sheehan and Gallimore2015; Saini & Ganesh Reference Saini and Ganesh2020)

(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})$ and $\boldsymbol {E}$ and $\boldsymbol {B}$ are the electric and magnetic fields, respectively, $r_{L}$ is the particle Larmor radius, q is charge of the particle, m is mass of the particle and $\theta$ is azimuthal direction in cylindrical co-ordinates. As the problem considered is non-relativistic, the electric field $\boldsymbol {E}$ is independent of the particle velocities. Hence, we first solve (1.2) without the $\boldsymbol {E}$-term and include the same while computing the full equations of motion, as shown below.

Let us start with the Maxwell equation $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {B} =0$, which can be written in simplified form using the approximation that magnetized particles are small current loops orbiting around magnetic field lines and $B_{z}$ is assumed to be uniform over the Larmor radius $r_{L}$, resulting in

(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) (without $\boldsymbol {E}$ term), 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.5)\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} \end{gather}

(1.6)\begin{gather}\frac{{{\rm d}}v_{r}}{{{\rm d}}t}=0 . \end{gather}

It is pertinent to note the following: (1.6) indicates that $v_{r}$ is a non-zero constant in our simulations, which implies that the magnetic field, $\boldsymbol {B}$, does not alter the radial velocity component. Also, it has been assumed implicitly that particles are in the frame of reference along the magnetic field lines. The force equations (now including the electric field) in the frame of the magnetic field lines are

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

(1.8)\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 as direction along the magnetic field line. The azimuthal velocity ($v_{\theta }$) is the velocity perpendicular to the magnetic field line ($v_{\perp }$). The forces (1.7) and (1.8) are solved using the Boris particle pusher algorithm. The initial half-push is performed with an electric field. This is followed by the magnetic field push and the final half-push is performed by the electric field. To solve the magnetic field force on the particles, the energy conservation principle is used because the work done by the magnetic field on charged particles is zero. Because of energy conservation, one can write

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

After simplification of the parallel velocity component, $v_{\parallel }$ at the $n+1$th time step can be written as

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

where $\zeta =({1}/{2B})({{{\rm d}}B}/{{{\rm d}}s})$ and $K$ are constant over $\triangle t$ (Appendix A).

An asymmetric inductive excitation scheme is used to energize the electrons (Meige, Boswell & Charles Reference Meige, Boswell and Charles2005; Ebersohn et al. Reference Ebersohn, Sheehan and Gallimore2015; Saini & Ganesh Reference Saini and Ganesh2020). An RF electric field $E_{y}$ of frequency $f={\omega }/{2 {\rm \pi}}=13.56$ MHz is applied normal to the axial direction to heat the electrons. The $\hat {y}$-direction is in the $\hat {r}-\hat {\theta }$ plane, which is perpendicular to the axial direction. The amplitude of this RF electric field is chosen to avoid overheating of the electrons. Hence, the RF current density $J_{0}$ is used as the control parameter for the RF electric field amplitude $E_{y}$. The total current density in the source region is

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

(1.12)\begin{gather}J_{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 the electron flux along the $y$-axis. Here, $\varGamma _{e,y}$ is

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

using Maxwell's equations, which results in

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

where $J_{0}$ is the RF current density, $L_{source}$ is the length of the heating region and $N_{\sigma }$ is the weight of a macro-particle.

An energetic electron collides with a neutral Ar fuel gas atom and further generates an electron–ion pair via ionization collision (Verboncoeur et al. Reference Verboncoeur, Alves, Vahedi and Birdsall1993; Vahedi & Surendra Reference Vahedi and Surendra1995; Frignani & Grasso Reference Frignani and Grasso2006). These particles evolve in time, inducing further ionization, or remain as bulk species and, after a sufficiently large simulation time, the system attains a steady state. Once steady state is achieved, desired quantities are stored to understand the late time behaviour of the system. It is important to note that, in our PIC-MCC simulations, electron–neutral and ion–neutral collisions comprising of elastic, ionization and excitation collisions as well as charge exchange collisions have been included, whereas electron–electron and electron–ion and ion–ion collisions have been neglected. For low energy electrons, scattering is more isotropic and energetic electrons witness forward scattering (small scattering angles). Here, we have used direct analytical formulae for argon collision cross-sections following the strategy used in the XOOPIC code (Frignani & Grasso Reference Frignani and Grasso2006). These formulae are known to accurately reproduce a database for argon collision cross-sections and are, at the same time, easy to implement.

The spatial magnetic field profile used in our simulations 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}}{l_{b}^{2}}\right)^{3/2}} & \mbox{if } 0.6\times L < z < 0.4\times L \end{cases}, \end{equation}

where $L$ is the total axial length of the system, $l_{b}$ is a control parameter for the magnetic field expansion rate and $B_{0}$ is the magnetic field amplitude. We have chosen $z_{0}=0.6\times L$ on the right-hand side of the heating region and $z_{0}=0.4\times L$ on the left-hand side of the heating region. With decreasing $l_{b}$, the magnetic field gradient increases in the expansion chamber and this results in rapid expansion of plasma charged particles in an experimental chamber (Baalrud et al. Reference Baalrud, Lafleur, Boswell and Charles2011; Ebersohn et al. Reference Ebersohn, Sheehan and Gallimore2015). Figure 3 shows the spatial magnetic field profile for different $l_{b}$ values having a symmetric magnetic field expansion profile on both sides of the plasma heating region. For the magnetic field strength and electron temperatures considered in this work, a typical average Larmor radius is approximately $1.9 \times 10^{-4}$ m. Here, in our simulation, electrons are regarded as magnetized and the strength of the bulk potential depends upon the magnetic field gradient in the expansion region.

Figure 3.

The spatial magnetic field profile with different $l_{b}$ for a symmetric magnetic field expansion. Here, the heating region starts from $z=0.08$ and ends at $z=0.12$.

This model does not include the neutral flow dynamics, although it is a critical parameter in a working device (Takao & Takahashi Reference Takao and Takahashi2015; Takahashi, Takao & Ando Reference Takahashi, Takao and Ando2016; Takahashi Reference Takahashi2021). Here, we also do not include the plasma plume interaction with the spacecraft body and charge particle detachment from the magnetic field lines, which are the critical phenomena for thruster performance. As indicated earlier, our 1D3V PIC-MCC solver is written in a Cartesian geometry. The HPT is a cylindrical device and geometrical forces (inertial forces) are included indirectly using model equations (Delzanno & Camporeale Reference Delzanno and Camporeale2013; Ebersohn et al. Reference Ebersohn, Sheehan, Gallimore and Shebalin2017). It may be noted that, in the magnetic field frame of reference, while we need only two velocity components, namely $v_{\perp }$ and $v_\parallel$, 3 velocity components are required for calculating the particle collision dynamics. We have used the energy conservation principle to calculate the magnetic field force on the charged particle and any change in $v_{\parallel }$ would be reflected in both perpendicular directions, $v_{\perp }=\sqrt {v_{x}^{2}+v_{y}^{2}}$.