Fluid-Kinetic-hybrid simulation for ion thruster using polymorphic particles

A hybrid Fluid-Kinetic model based on a polymorphic Particle-in-Cell method is developed for the simulation of plasmas in ion thruster. The numerical scheme is suited for scenarios in which thermal ions or electrons undergo strong acceleration. The polymorphic PIC model is tested for the extraction region of an ion thruster. In this paper, we report the results achieved by implementing the scheme on a single extraction orifice and compare our results with the simulation from known IGUN software for ion extraction.


Introduction
A Particle-in-Cell (PIC) method is extensively used for the simulation of plasmas since its early development in 1961 by Dawson [1,2].The PIC methods are sub-divided into kinetic PIC involving kinetic equation for collisionless plasma and the fluid-PIC that approximates plasma as fluid and quantities like mass, velocity and temperature are solved on Eulerian mesh using hydrodynamic equations [3].Both the kinetic-PIC and fluid-PIC uses computational particle to represent plasma or particle streams.In kinetic-PIC one defines "super-particles" representing ensemble of physical particle whereas in fluid-PIC there are computational-particles along with grid quantities.The kinetic description provides a more accurate treatment of many local and quasi-local processes at a cost of intensive computation.On the other hand, a fluid approach is suitable for treating large scale properties of plasmas involving mass, momentum, and flux transport with certain assumptions concerning the charge neutrality and the distribution function.A hybrid model tries to achieve the balance between two methods by treating one of the plasma components as fluid typically the electrons and other components with a kinetic description [4,5].The MHD-EPIC (MHD with Embedded PIC) is another example where kinetic-PIC was coupled with fluid-PIC.In this case results of Magnetohydrodynamics (MHD) results were fed to kinetic-PIC in specific domain of interests [6].
The description of plasma in ion thrusters, or in general, radio-frequency discharges pose serious challenges considering heavy computational efforts [7].In the earlier publication, we demonstrated the simulation of radio-frequency plasma employing a kinetic PIC method in 3 dimensions using parallelised MPI code [8].But even these simulations are restricted to small scale thrusters up to a cubic capacity of about 10000 mm 3 with plasma density up to 10 16 /m 3 .A commercial version of thrusters requires a more efficient and faster approach for plasma simulation.The hybrid simulation is frequently being used for thruster simulations [9,10].The treatment of the fluid model implies that particles follow the Maxwell distribution in order to define the temperature of the plasma and that the plasma is quasi-neutral [11,12].However, the transition from quasi-neutral plasma in the discharge chamber to a fully ionized extraction beam, forming the thruster plume cannot be treated with a single-fluid model.It has been also shown that the electrons in such cases do not entirely follow the Maxwell-Boltzmann distribution [8].The evolution of the non-Maxwellian distribution of electrons in such plasma entices us to look beyond a hybrid model to incorporate a non-thermal phase-space distribution.
The current work is based on the framework proposed as polymorphic nature of computational particle that allows changing fluid particle properties to the kinetic depending on the velocity distribution criteria [13].The main motivation comes from the fact that the plasma inside small scale ion thrusters can not be diagoned directly.The temperature and densities can only be estimated from extracted ions from plasma sheath near extraction grids.Unfortunately for kinetic simulation this poses bigger challenge as the velocity distribution inside and that of extracted ions is vastely different, leading to computational challenges in terms of different time and spatial scale in those regions.This is also true while simulating electrons and ions in any accelerating region or near vessel where electrons have higher velocities than ions.The non homogeneous grids and different times scale demands for higher computational cost and efforts.The proposed method follows a smooth transition from fluid-PIC approach to Kinetic-PIC approach depending on the changing velocity distribution function.In general this method can be extended to boundary plasmas or electron clouds in accelerator physics.
In this paper, we consider a plasma simulation in electrostatic fields near the extraction region of the ion thruster.In this region, a plasma sheath is formed, and the thermalized ions are accelerated and extracted through a hole with high longitudinal velocity.This scheme can also be adapted for multiple grids and adaptive time scale for faster implementation of ion extraction simulation.

Fluid-Kinetic model and particle polymorphism
We consider a plasma composed of two species, electrons and ions.The fluid equations for plasma specie α in terms of mass density ρ m , charge density ρ c , fluid bulk velocity u , internal energy I and pressure p, subjected to an electrostatic field, are written as (1) The first three equations describe the evolution of plasma through conservation of momentum and energy, whereas the last equation is called Equation of State (EoS).
The initial conditions for these quantities on grid points are achieved by homogeneous distributions on species in space and Maxwellian velocity distribution.Using a first-order weighting scheme the fluid bulk velocity, temperature, pressure and internal energy are calculated on grid points as, where m p , q p , v p and e p correspond to particle mass, charge, velocity and internal energy at position x p .The weighting function W can be defined using zeroth-order nearest grid point or volumic weighting in a 3-dimensional case.
The electric field over the domain is calculated using Poisson's equation, where all the quantities, potential, charge densities and electric fields ( φ, ρ c , E ) are defined on grid points.The solutions to these equations are obtained using Gauss-Seidel iterative method.The crux of the fluid-PIC method is using computational particles to calculate the convective term u • ∇ of the governing equation.In PIC method individual particles or super particles are moved using Lagrangian mover to track their individual positions and velocities in the continuous phase space.Where as the evolution of grid quantities are described by Eulerian mover e.g.fluid bulk velocity, energy and pressure etc on the stationary mesh.
For Eulerian mover the convective terms are removed from the momentum and energy equations and then rewritten as, These equations are discretized in 3-dimensional space and time as (4)

and
where n denotes the time level of the update.Note that i, j, k represent indices along x, y and z axes respectively and not written explicitly if only one of them changes.The fluid particle quantities are updated through ordinary differential equations whereas the kinetic particles follow The polymorphism of the PIC particle is defined through the velocity of the individual particle.As mentioned earlier the particles are initially defined with Maxwell's distribution of velocities in the domain under consideration.Thus at the beginning, all the particles are assigned the "Fluid" type which can be changed according to a predefined rule, which utilizes the velocity of the individual particle.
For example, if the velocity |v p | of a particle exceeding a multiple of the local ther- mal velocity v th is converted and assigned the "Kinetic" type where c 1 is some constant and v th is the thermal velocity.This formula provides gen- eral conditions to convert individual particles e.g."Fluid" to "Kinetic" type.One can define different rules in terms of velocity depending on the problem under study.But we have chosen here to switch from fluid to kinetic particles when fluid particles reach a threshold velocity or acceleration.In [13] authors found that a multiple of local thermal (11) velocity is a convenient threshold velocity.Thus corresponding equations of motion are switched from fluid to kinetic as per the rule dictates.

Fluid -Kinetic equivalency
We start with a simple case of plasma sheath formation in one dimension.A linear homogeneous ensemble of electrons and ions is defined with Maxwellian velocity distribution.The endpoints are fixed to ground potential.Under these conditions, we carry out simulations using three different models kinetic-only, fluid-only and polymorphic Fluid-Kinetic-PIC.Due to the higher mobility of electrons higher number of electrons are lost than ions leading to a positive potential with respect to the wall.The Debye shielding confines the potential variation to a layer called plasma sheath.Figure 1 shows the formation of a plasma sheath in terms of positive potential and the density profile of electrons using the three different models.Electrons were initialized with a temperature of 3 eV and Xe + ions of 1 eV.
In the case of the Fluid-Kinetic polymorphic PIC model, the threshold velocity was set at three times the thermal velocity for each species i.e. c 1 = 3 in Eq. 18.Initially, both ions and electrons are initialised with thermal velocities setting "Fluid" attribute.Since electrons are faster and have higher mobility, more electrons are converted from "Fluid" to "Kinetic".It can be seen that all three models give similar profiles in terms of density profile and potential distribution.This simulation provides the sanity check that polymorphic model still gives same results for given density and temperature of plasma as that of others even though execution time is different.

Simulation of ion thruster : RIT-4.0
Figure 2 shows the setup for the 3-dimensional simulation domain projected on an x-z plane.The model follows approximate geometry and parameters which are designed for ion thruster RIT-4.0.Homogeneous particle distribution is predefined in a plasma chamber with thermal velocities following Maxwell's distribution function.Typical geometric parameters are listed in Table 1.
All the geometric entities are defined in the 3D domain with grid points defined as structured cubic mesh.Declaration of particular node points as a material with fixed potential points also serve as boundary conditions.The structured cubic mesh consists of grid points defined at equi-distant locations, 0.2 mm apart in each direction.The homogeneous cubic grid model also know as sugar cube model leads to small aberrations where curved boundaries are designed.The smoothness in the electric field is achieved by calculation of potential at grid point using Shortley -Weller method that modifies the Poisson equation [14].The plasma chamber is defined with open boundary condition in longitudinal z-direction i.e. x-y plane at z = 0 .Radially the plasma is confined and held at the same potential as that of the screening electrode ( φ screen = 1200 V ).The screening grid is located at a distance of 2.0 mm from the origin.The accelerator grid is held at a potential of φ accel = −150 V located at a longitudinal distance of z = 3.25 mm from the origin.The ground electrode is located at a distance of 5.0 mm downstream.The detector plane is defined at a distance of 9.0 mm.
A quasi-neutral plasma with density 1 × 10 17 /m 3 is generated with homogeneous spatial distribution.A time-step of 1 × 10 −10 s is chosen for both ions and electrons.Initially, all the particles are initialized inside the plasma chamber with fluid properties.The velocities follow a Maxwellian distribution with thermal velocity v th = (3k b T /m) for each species.Along with position and velocity, each particle is defined with the "Fluid" attribute through an additional variable "Kflag" indicating its state.
Figure 3 shows ion distribution colour coded with velocities.As expected ions near the extraction aperture are subjected to an accelerating force, thereby increasing velocity in the longitudinal direction.Some of these particles fulfil the condition as stated in Eq. ( 18), their attribute is then changed to Kinetic flag, and consequently, they follow the kinetic equation of motion.
Figure 4 a) shows a qualitative comparison between ion distributions simulated using two different methods.The upper part shows the results of a simulation using a standard Particle-in-Cell code with kinetic equations alone, whereas the lower part shows the results using polymorphism for the ions.Both distributions show qualitatively the same density distribution in the simulation domain.
For a better comparison, we define Root Mean Square (RMS) emittance at the detector plate as (19) where defines the second central moment of the distribution, and x ′ = v x /v z denotes the angle between transverse to longitudinal velocity.The emittance is calculated only for vertical x-plane.Figure 4 b) compares emittance as a function of ion temperature in both cases.The functional variation in the emittance agrees quiet well between both methods.

Comparison with IGUN
IGUN is a well-known software developed at Goethe University of Frankfurt am Main, which has been used multiple times for designing ion extraction systems in the field of accelerator physics and within the ion thruster community [15,16].Simulations were performed using IGUN with the above listed (Table 1) geometric and plasma parameters to compare the results to our polymorphic-PIC method.In the case of IGUN software, the input parameter includes electron temperature, plasma density, ion species and electrode voltages.In this simulation code, Debye lengths are matched with our simulation parameter.About 10000 test particles i.e. ions are extracted and tracked to calculate the emittances and transmission factor.Figure 5 compares transmission factor ( β ) defined as a ratio of extracted current to Bohm current.Figure 5 a) compares the transmission factor as a function of ion density using two models.Except at a few points, the values simulated using polymorphic-PIC lie within the error range.Figure 5 b) compares the emittance as a function of temperature.Although the IGUN values are a little higher for the lower plasma temperatures than calculated by polymorphic PIC, the behaviour with respect to increasing temperature agrees with each other.The discrepancies can be attributed to different method of charge deposition at grid points and the implementation of sharp electrodes at grid points; the methodology is not known in the case of licensed IGUN software.

Conclusion
We developed a fluid-kinetic hybrid code with polymorphism for plasma simulation in full 3-dimensional cartesian coordinate system.This method is well suited for the simulation of ion thrusters.On the one hand, kinetic description of plasmas is computationally heavy requiring large number of particles, while on the other hand fluid equations rely on assumed velocity distribution failing to capture important physical phenomena.The fluid PIC require relatively low number of particles compared to the kinetic PIC.
The electric potential distribution near debye sheath can be analytically calculated in one dimensional case.We simulated the potential using two methods and compared the average deviation in the electric potential as a function of number of macro particles.The system is evolved in steady state already in each case.Figure 6 shows that in the case of poly-Fluid calculation lower number of macro particles suffice to achieve acceptable accuracy predicted analytically.This value in almost 60% to that used in Kinetic simulation.When calculations for extractions inclusive MCC were performed on a highperformance cluster using a MPI (Message-Passing-Interface) domain decomposition on Intel Xeon Skylake Gold 6148, with 24 cores on a single node it required 8 hours as compared 12hrs for Kinetic only simulation.
The polymorphic PIC model put Fluid and Kinetic descriptions in the same framework.In this work, we have implemented polymorphism based on the velocity of individual particles.In the systems like ion thrusters or ion sources, it is intuitive that a large number of particles in the plasma will be in a thermal state and only a small number near extraction will have non-Maxwellian distribution.This fact can be effectively used to further ease the simulation strategy.
Figure 7 shows ion distribution colour coded with a "Kinetic" flag.In this example, polymorphism was defined as a spatial function.A spatial grid-based approach is adapted in such a way that for region z < 1.8 mm , all the particles are attributed fluid properties whereas the region near extraction carries the kinetic attribute.This approach is still under investigation but could be useful in implementing multi grid systems with different resolutions and field equations.With this approach plasma inside gridded ion thruster can be calculated with using entirely fluid-PIC and computational particle near grid are morphed into super-particles for kinetic-PIC and consequently extracted for further calculation without need for separate coupling of to method as in MHD-EPIC.

Fig. 1 a
Fig. 1 a Density Profile in one dimensional plasma and b Plasma potential distribution for three different models (colour-online)

Fig. 2 a
Fig. 2 a Simulation model for ion thruster showing three-gridded extraction system and location of detection plane.b Potential distribution projected on the x-z plane (colour-online)

Fig. 3 Fig. 4 a
Fig. 3 Spatial homogeneous distribution of ions after a few simulation steps showing accelerated ions near extraction aperture colour coded with respect to speeds (colour-online)

Fig. 5 aFig. 6
Fig. 5 a Transmission factor as a function of density and b Emittance as a function of temperature comparing polymorphic model with IGUN (colour-online)

Fig. 7
Fig. 7 Ion distribution colour coded with Kinetic flag as a spatial function instead of thermal velocity (colour-online)

Table 1
Parameters for simulation of ion extractiona Although the radius of the plasma chamber in RIT-4.0 is about 21 mm we have restricted the plasma chamber radially to simulate single aperture extraction and avoided large mesh in the transverse direction