Kinetic simulation of ion thruster plume neutralization in a vacuum chamber

The electrical environment of a ground vacuum testing chamber creates facility effects for gridded ion thrusters. For example, it is well known that the plume from the thruster generates current paths that are very different from what occurs in space, and the neutralization of this plume is also different. For reasons such as this, it is important to clarify how the experimental testing environment affects plasma flows, but understanding this effect solely through ground experiments is difficult. To that end, this study utilizes particle-in-cell and direct simulation Monte Carlo methods to simulate xenon beam ions and electrons emitted from a neutralizer. First, we compare simulations conducted within the chamber to those conducted in space, demonstrating that grounded chamber walls increase the electric potential and electron temperature. Next, we investigate the impact of the neutralizer’s position and the background pressure on the plume in the vacuum chamber. We find that as the neutralizer position moves closer to the location of maximum potential, more electrons are extracted, resulting in increased neutralization of the plume. We also observe that high background pressure generates slow charge-exchange ions, creating ion sheaths on the side walls that alter ion current paths. Finally, we discuss how the potential at the thruster and neutralizer exits affects the plume. The relative potential of the neutralizer to the vacuum chamber wall is observed to significantly influence the behavior of the electrons, thereby altering the degree of plume neutralization. These findings are shown to be consistent with experimental results in the literature and demonstrate the promise of high-performance simulation.


Introduction
Gridded ion thrusters (GITs), an electric propulsion device commonly used due to its high-specific impulse (> 3000 s) [1], generate thrust by selectively extracting ions using ion optics (grids).However, the performance of GITs relies not only on the ion optics but also on the neutralizer attached near the ion source location.Thermal electrons are released from the neutralizer to mitigate spacecraft charging caused by the space-charge effect of excess ion emission from the thruster, and the coupling between the two sources greatly affects the performance of the neutralizer.Research on this coupling of neutralizers dates back to the SERT II satellite [2], which demonstrated that the operation of the neutralizer reduced the electric potential of the ion plume and suppressed the decrease in spacecraft floating potential in space [3].Nakayama et al. [4] also investigated the ion current in a ground vacuum chamber by changing the neutralizer electron current and revealed that more ions returned to the thruster exit and body when the neutralizer electron current was insufficient.
Past GIT plume experiments have mainly focused on beam ion current and energy, which directly relate to thrust, unlike Hall effect thrusters, where the electron transport from the external neutralizer is important for plasma generation [5,6].In addition to current and energy, Polansky et al. [7] studied ion and electron number densities downstream in the radial direction, and Conde et al. [8] investigated the energy distribution of ions and electrons in a GIT plume.However, to date, few studies have experimentally investigated the detailed electron motion.
Table 1: Previous study of the gridded ion thruster plume using a fully kinetic PIC approach.
(co-located position), with the plume immediately neutralized just downstream from the thruster exit.In actual GITs, however, a large potential gradient occurs in principle because the ions and electrons are supplied from separate positions.Other work listed in Table 1, Refs.[19][20][21][22][23], studied plume neutralization for such shifted electron source position cases as described "external" or "internal" in the table.It was observed that the plume potential significantly changed when comparing co-located and external source cases, even for the same electron source size and density [20,23].The external source cases have relatively smaller R 0 /λ D0 than the co-located cases, where R 0 and λ D0 are the thruster exit radius and the initial Debye length, respectively.This means that they target smaller thrusters or lower-density plasmas because the electron density at the shifted neutralizer is normally much larger than the ion beam density.Hence, the minimum mesh size becomes smaller than the co-located cases, increasing computational costs.
Table 1 also shows that many simulations model space operations [9,10,[12][13][14][15][16][17][19][20][21][22][23].Although the in-space environment has almost an infinite volume for plume expansion, an otherwise-trapped electron may reach the end of the computational domain before it can reverse direction due to finite computational resources.Thus, Refs.[14,21,23] have focused on developing electron boundary conditions for fully kinetic simulations.As can be seen from the table, however, few studies have been conducted on ground chamber tests.
Ground tests of electric propulsion systems can introduce significant uncertainty in the onorbit performance prediction because of facility effects [26].For example, Nishii et al. [18] simulated the contamination caused by carbon backsputtering for different sputter models and plume conditions and demonstrated the importance of combining backsputtering and PIC plume simulations.Hu et al. [12] studied electrical effects by simulating a proton ion beam with different beam radii and showed that the beam width became 20 times larger than the thruster radius where ions reached their maximum velocity.Their work suggested that the vacuum chamber facility may significantly affect the plasma beam by prematurely terminating its expansion if the chamber size is less than 20 times the beam radius.However, since xenon has a higher mass than that of a proton (their assumed ion mass), the plume diffusion in their simulation would be smaller than that of an actual GIT.In addition, since the ion and electron sources were in the same position, which is not realistically possible, the effect of coupling of distant sources was neglected.Therefore, the electrical effects that affect the prediction of the neutralizer coupling voltage, the neutralizer bias voltage from the thruster common, and the plume divergence have not been fully investigated.
The primary objective of this study is to investigate electrical facility effects related to GITs using a fully kinetic simulation.We use an in-house developed 3-D PIC-DSMC solver, Cuda-based Hybrid Approach for Octree Simulations (CHAOS) [20], which enables us to solve for a steadystate GIT plume with relatively large R 0 /λ D0 for the true mass ratio between xenon ions and electrons, as shown in Table 1 (Refs.[16-18, 20, 23]).This study specifically seeks to understand the effect of the vacuum chamber boundary condition, the neutralizer exit location, and the finite background neutral density on the GIT plume, including beam ions and neutralizer electrons.
Additionally, the biased (non-zero) potential that appears in the ground chamber during GIT operation has not previously been investigated in numerical simulations.For instance, in a twogrid GIT system, there is no deceleration grid (decel grid), which is used to reduce the CEX ion backflow that can cause acceleration grid (accel grid) erosion [27].In such cases, the accel grid at a negative voltage is exposed directly to the plasma plume.Second, the neutralizer coupling voltage is normally altered to increase the neutralizer current [2,4,7].Taking into account the keeper positive voltage [2], which is crucial to the extraction of electrons from the hollow cathode, the neutralizer exit voltage should differ from the thruster common voltage.Therefore, we also investigate how such electrical potential boundaries affect the GIT plume in the vacuum chamber.
The outline of the remainder of this article is as follows.Section 2 reviews our plasma modeling approach and describes the boundary conditions implemented in CHAOS.Section 3 explains the geometry, species, and numerical conditions of the study and selection of variable parameters.Finally, we present and discuss comparisons between the simulated cases, including the effect of the simulation boundary condition (Section 4), neutralizer locations (Section 5.1), background neutral particles (Section 5.2), and electric potential at the thruster (Section 6.1) and neutralizer (Section 6.2) exits.

PIC and DSMC Modules and Their Coupling
In this section, we briefly discuss the computational framework implemented in CHAOS to couple the PIC and DSMC approaches in order to calculate the self-consistent electric field, taking into account the reaction between ions and neutral particles.In the EP plume, collisions and electric fields operate on significantly different time and length scales, differing by at least two orders of magnitude.To mitigate the impact of these differences, CHAOS employs several computational techniques as previously described in our earlier papers [16-18, 20, 23].
The DSMC module models three types of collisions: momentum exchange (MEX) collisions between Xe-Xe and Xe-Xe + , and CEX collisions between Xe-Xe + .The collision cross sections for MEX between neutral particles and MEX and CEX collisions between neutral particles and ions are obtained from Refs.[28], and [29].The no-time-counter collision scheme proposed by Serikov et al. [30] is used in this study since it accounts for the disparate timesteps and weighting factors of ions and neutral particles.The neutral particles move only when the DSMC module is executed, every 100 PIC timesteps, while the ions and electrons move every iteration.
In the PIC module, the electric potential is calculated using an explicit PIC technique.In the fully kinetic approach, the electric field, E, is self-consistently solved by: where ρ is the charge density, e is the elementary charge, n i and n e are the number density of ions and electrons, φ is the electric potential, and 0 is the permittivity of free space.A finite volume approach based on an unstructured octree grid is used to solve Eq. ( 2).
CHAOS has a number of PIC and DSMC coupling algorithms that save computational effort.First, using a Morton Z-curve, CHAOS constructs two separate grids with a linearized forest of octrees (FOT) in the PIC and DSMC modules, respectively, because the mean free path, λ, and Debye length, λ D , differ by at least three orders of magnitude.The FOT for DSMC (C-FOT) is constructed to resolve the local mean free path, while the FOT for PIC (E-FOT) is constructed to resolve the local Debye length, where the refined cell size x < λ D .We apply an adaptive mesh refinement method since the number density can vary widely in the computational domain.Both C-E-FOTs are reconstructed every 20,000 iterations before sampling starts.Second, weighting factors, W , are utilized to increase the number of charged computational particles compared to the neutral particles due to the disparate length scales of the C-and E-FOTs and disparate number densities of the neutral particles and CEX ions.In this study, the ratio of neutral and ion (W n /W i ) is set at 200,000.Third, time-slicing of the DSMC and PIC modules and speciesdependent timesteps are implemented due to the different timescales for collision and plasma frequencies.The positions of the neutral particles, ions, and electrons are updated with timesteps of ∆t n ∆t i = ∆t e to reconcile these disparate timescales.In this study, we use a timestep of ∆t n = 1.12 × 10 −4 s for neutral particles and ∆t i,e = 2.8 × 10 −10 s for both ions and electrons.

Boundary Conditions
To satisfy the objectives of this study, we use both in-space and in-chamber boundary conditions (BCs) for the outer edge domain boundary.For the in-space simulation, the charge-conserving energy-based BC (CCE BC) developed by Jambunathan and Levin [23] is used for the downstream boundary (z = 0.8 m), and the buffer BC is used for the other boundaries to simulate the infinite expansion of the thruster plume, similar to our previous calculations (Refs [16,23]).The buffer BC simulates the inflow of electrons from outside by placing a buffer region outside the computational boundary and copying the particles inside the boundary out to a distance of λ D0 beyond it.The CCE BC specularly reflects some electrons arriving at the edge of the computational domain and eliminates others using the following approach.The baseline total charge, Q 0 , and the average electron kinetic energy in the computational domain, E e , are obtained at the time step just before the beam ions reach the downstream domain boundary.In subsequent timesteps, when the total charge in the domain is less than Q 0 , electrons with energies less than E e are specularly reflected.In Ref. [23], it was verified that these BCs satisfy the requirements needed for the plume modeling by changing the size of the domain (see Table II in Ref. [23]).In terms of electrical boundary conditions, the inhomogeneous Neumann BC for the electric potential is implemented on all domain boundaries in the in-space simulations.For each boundary, the normal potential gradient (∂φ/∂n) bc is computed based on the current density through the boundary as follows; where N i,bc and N e,bc are the number of ions and electrons that cross the boundary, and A bc is the area of the boundary.
For the in-chamber boundary conditions, we use a fully diffuse reflection condition with a 300 K accommodation for ions and neutral particles on the chamber walls at the edge of the computational domain.The walls absorb all electrons by removing them from the domain and neutralize all incident ions by returning them into the domain as neutral particles with a temperature of 300 K.The 0 V Dirichlet BC is implemented on every boundary surface in the PIC module.The plasma screen, which is the housing of the thruster assembly and is normally electrically grounded, has the same BC as the vacuum chamber walls.We also assign a chargeabsorbing BC for particles impinging on the thruster and neutralizer exits and a Dirichlet potential BC for the electric potential with a baseline value of 0 V.The detailed settings about the potential boundaries are described in Section 3.2.

Calculation Geometry and Species
Figures 1 and 2 show the three-and two-dimensional schematics of the in-chamber cases investigated in this study.Note that only the geometry of the neutralizer differs in Figs. 1 and  2. The GIT is placed in a cubic vacuum chamber with a length of 0.8 m per side.In this study, only a half of the domain is simulated due to symmetry to save computational effort, i.e., a specular reflection BC and a Neumann BC (∂φ/∂n = 0) are implemented on the x = 0 m plane.Numerical pumps, shown as green volumes, are installed at all corners of the downstream face to remove heavy neutral particles from the vacuum chamber.Computational particles entering the numerical pump volume are deleted from the calculation, utilizing the same method employed in our previous studies [18,31].The cross-sectional area of the numerical pump is assumed to be 1.25 × 1.25 cm 2 , which produces a typical vacuum chamber pressure of 10 −6 Torr.The dimension of the thruster, which emits Xe neutrals and Xe + ions, is the same as that used in our previous calculations of an ion thruster system [16][17][18]23].The center position of the cylindrical thruster exit with a radius, R 0 , of 0.0625 m is located at (x 0 , y 0 , z 0 ) = (0.0, 0.4, 0.1) m.The domain contains the plasma screen offset from the inlet plane at z = 0.1 m.This study examines two neutralizer sizes and number densities.The first configuration is designated as "A," which has the same characteristics as in our previous work [17], as illustrated in Fig 1; (a): three-dimensional view, (b): two-dimensional view seen in the x-y plane, and (c): in the z-y plane.This neutralizer has the same radius as the thruster and is placed side by side on the same plane.In addition to the type A configuration, this study analyzes a more realistic configuration designated as "B."This neutralizer has a higher electron density and a smaller outlet downstream from the thruster, as shown in Fig 2 ; (a): three-dimensional view, (b): two-dimensional view seen in the x-y plane, and (c): in the z-y plane.The neutralizer exit radius, R e0 , is set to 1 cm based on the size of a typical neutralizer [32] and the exit electron number density is 32 times greater than type A to obtain the same current.Unlike the type A cases, the neutralizer exit shifts far from the thruster exit by a value of d 0 in the z-direction to investigate the effect of the electron exit position on the plume.
Table 2 summarizes the conditions of each species at the thruster and neutralizer exits.Similar to previous mesothermal studies [11,20,23], we chose a ratio of the initial ion temperature, T i0 , to the initial electron temperature, T e0 , of 0.01.All species are initialized at their sources with a stationary half-Maxwellian distribution in the streamwise and full-Maxwellian in the cross-stream directions.The reference plasma number density is considered to be the same as the ion number density at the thruster exit of n 0 = 1.0 × 10 13 m −3 .These selected values for T e0 and n 0 result in an initial Debye length, λ D0 = 3.32 × 10 −3 m, initial electron plasma frequency, ω pe = 1.78 × 10 8 rad/s, and initial electron thermal velocity, v te0 = 592, 892 m/s.In the kinetic simulations, the ion and electron timesteps should follow the requirements of ∆t < 0.1ω pe .In this work, we use a timestep of ∆t = 0.05ω −1 peo = 2.8 × 10 −10 s for ions and electrons.The superparticle parameter, F num , is set at 2,500 for all simulations such that there are at least 15 particles per species per cell for both the C-FOT and E-FOT.As a result, the total number of computational ion and electron particles is about 20 million at steady state.Table 2 also gives the Xe neutral particle parameters for the cases with a background pressure; otherwise, the plume is assumed to be collisionless.The exit number density of neutrals is n n0 = 1.0 × 10 17 m −3 , which is the typical order for actual GITs giving a total number of computational neutral particles of about five million at steady state.The source centers are (x 0 , y 0 , z 0 ) for Xe and Xe + and (x e0 , y e0 , z e0 ) for e − .
In this study, the following assumptions are made to simplify the model.Since the holes on the thruster optics are not modeled, ions are uniformly emitted from the thruster wall surface.In reality, the holes cause variations in ion density, and electrons present to some extent inside the hole [33,34].Unlike a hollow cathode, electrons are emitted from the neutralizer wall surface similar to a filament neutralizer [7,8], a direct emission neutralizer [35], or a diode mode neutralizer [36].Finally, collisions between Xe and e − , Coulomb collisions between Xe + and e − , the presence of multiply-charged ions (i.e., Xe ++ or more), and electron emission due to secondary electron emission (SEE) and ion-induced electron emission (IIEE) are neglected.

Case Descriptions
We test eight conditions to investigate the effects of in-space versus in-chamber geometries, neutralizer locations, and the presence of background neutral particles, as defined in Table 3.The in-space geometry simulation, designated as "1," is compared to the in-chamber geometry simulation, designated as "2," for the neutralizer types A and B. The difference in the BC settings was described in Section 2.2.Then we discuss the ion beam and neutralizer coupling, changing the neutralizer exit positions at four locations: 2B-0, 2B-1/4, and 2B-1/2 cases indicating d 0 R 0 = 0, , respectively, where d 0 is the distance between the neutralizer exit and thruster exit, and R 0 is the thruster radius.The above-mentioned cases do not model neutrals to differentiate the electrical effect from the high-background pressure effect.Therefore, we also simulate a case with a finite background pressure designated as 2B-BP with d 0 R 0 = 1.An electrical schematic diagram for GIT ground operation modeled in this study is shown in Fig. 3. Ions with a potential of V dc generated in the discharge chamber inside the thruster are accelerated by ion optics and emitted from the external grid with a voltage of V th , where the beam ion kinetic energy in the axial direction corresponds to V dc − V th .Since this study uses 30,000 m/s for the beam ion velocity when V th = 0 V, V dc is set at 612 V, assuming that there is only ion motion outside the thruster, for simplicity.This study further investigates how the thruster plume is affected by using similar V dc and V th potential conditions to that in actual experiments.Table 4 shows these three additional conditions with respect to the 2B-BP case.The 2B-ACC case corresponds to the situation where there is no outermost decel grid so that the accel grid with V th = −200 V is exposed to the plasma plume.The exit density and velocity for the 2B-ACC case are corrected to 8.68 × 10 12 m −3 and 34,553 m/s, respectively, because the incoming ions are considered to have a 200 V higher axial energy compared to the baseline case 2B-BP.The 2B-NM and 2B-NP cases simulate the case where the neutralizer is negatively and positively biased relative to the chamber and plasma screen by V ne , which corresponds to the neutralizer coupling voltage applied in GIT experiments.: Electrical schematics of this study.I i0 and I e0 are the currents of ion emitted from the thruster and electrons emitted from the neutralizer, respectively.I ne is the current to the neutralizer exit, I th , I ps , and I vc are the currents to the thruster exit, plasma screen, and vacuum chamber, respectively.V dc , V th , and V ne are the biased voltages of the discharge chamber, thruster external gird, and neutralizer exit, respectively.
The ion bulk velocity is corrected to 34,553 m/s, and the ion exit density to 8.68 × 10 12 m −3 .All other conditions are the same as the 2B-BP case.
After the simulations reach a steady state, the current is calculated for ions and electrons by directly sampling the computational particles, as shown in Fig. 3. From the conservation laws, ion and electron currents are given by: CHAOS has multiple GPUs with MPI-Cuda parallelization strategies [20].This study uses 16 NVIDIA A100 GPUs on the Delta machine at the National Center for Supercomputing Applications for all cases.In the cases without neutral particles, we simulate 500,000 steps before sampling and then sample 200,000 steps to obtain the field macro-parameters and currents.In contrast, in the case with neutral particles, we simulate 5,000,000 steps prior to sampling due to the slow CEX particle motion and then sample 500,000 steps.The total simulation runtimes are; about 30 hours for the 1A and 1B cases, about 12 hours for the 2A, 2B, 2B-0, 2B-1/4, and 2B-1/2 cases, about 90 hours for the 2B-BP, 2B-ACC, 2B-NM, and 2B-NP cases.

Effect of Space vs. Ground-Based Chamber Conditions
This section presents the outer edge boundary effect between the in-space and in-chamber cases.Figure 4 shows the volume charge density, ρ, of the 2B case to highlight the three-dimensional plume structure, where ρ is normalized by the reference charge density ρ 0 = en 0 .Ion beams with a velocity of 30,000 m/s and thermal electrons with a temperature of 2 eV are emitted from the thruster and neutralizer exit in the positive z-direction, resulting in a maximum and minimum ρ near the respective exit points.As the ion beam moves downstream, it is neutralized by coupling with the electrons but also expands due to the positive space charge.In the three-dimensional diagram, we indicated two planes that will be used for comparing subsequent cases.
Not only the type B, as shown in Fig. 4, but also the type A neutralizer are examined with different BCs at the edge of the domain, i.e., 1A versus 2A and 1B versus 2B. Figure 5 shows the charge density in the x = 0 m plane, with color contours indicating the degree of neutralization of the plume and green arrows indicating the electron streamlines in the plane.Regardless of the neutralizer type, the in-chamber results show positive charge density near the downstream wall (0.7 m < z < 0.8 m) due to the formation of a sheath by the charge-absorbing wall at 0 V.In addition, the green arrows indicate a higher density of electron streamlines that flow downstream compared to the in-space case.These are an example of facility effects caused by the presence of grounded potential walls.It is also evident that the same type of neutralizer as in the previous study [23] (type A) and the more realistic type of neutralizer (type B) create different spatial distributions of the volume charge density both in the in-space and in-chamber cases, even though the electron emission is set to the same level.This suggests that the type and location of the neutralizer significantly affect the neutralization of the plume, as discussed further in Section 5.1.
Figure 6 shows the potential, φ, in the y = 0.4 m plane, where φ is normalized by kT e0 /e = 2 V.The potential in the y = 0.4 m plane increases from the thruster exit downstream to a maximum value and then decreases further downstream.When comparing the in-space and in-chamber cases, higher potentials are observed in the chamber case for both neutralizer types.This is also quantitatively shown in the line plot of the potential along the thruster axis as given in Fig 7.
The difference in the number of electrons present in the computational domain explains the reason for this behavior.Table 5 displays the number of computational ions and electrons in the entire domain at steady-state.In all cases, the number of electrons is lower than that of ions, but it is more pronounced in the chamber case.This is because, at the edge of the computational domain, the chamber absorbs all electrons, whereas in the space case, the CCE BC simulates an actual space condition in which there is a backflow of electrons that are trapped in a potential well formed by the ion beam.A, space (1A) A, chamber (2A) B, space (1B) B, chamber (2B) Figure 7: Electric potential along the thruster axis for the in-space vs. in-chamber cases for neutralizer type A and type B, respectively.
Since the behavior of electrons is key to understanding the difference in the spatial variation of φ, it is important to evaluate the electron velocity distribution functions (EVDFs) as well as the macro parameters of the electron flow indicated by the green arrows in Fig. 5. Figure 8 shows the z-direction EVDFs obtained by sampling computational electrons at (x, y, z) = (0, 0.4, 0.1625) m, i.e., at a position R 0 downstream along the thruster axis.The dotted lines and markers indicate the fraction of electrons in a velocity bin, where the velocity and temperature values are normalized by the reference electron thermal velocity, v te0 , of 838,782 m/s, and the electron temperature at the neutralizer exit, T e0 , of 23,200 K (2 eV), respectively.In all cases, there are clearly two distributions for the electrons.One is a thermalized distribution with a high temperature peaking at nearly w e /v te0 = 0, and the other is a cold electron flow with a peak around w e /v te0 = 5.This flow can be understood from the streamlines moving toward the lower right at (z, y) = (0.2, 0.4) m in Fig. 5.In other words, there are two types of electrons in the steady-state ion thruster plume: thermal electrons trapped around the plume core and electrons flowing downstream driven by the potential gradient.
The solid lines in Fig. 8 are fitting curves obtained by a one-dimensional Maxwellian distribution function of: where w e is the electron z-velocity, m e is electron mass, k is Boltzmann's constant, and T e is the electron temperature.The fitting range is -6 < w e /v te0 <3 in the z-direction since only bulk electrons are considered.The obtained normalized fitting temperatures T e /T e0 are of the order of about 10 and are larger for the in-chamber cases for both types of neutralizers.This is due to the difference in the maximum potential and the sheath at the chamber downstream wall of the plume creating a steep potential gradient that attracts and accelerates more electrons, as shown in Fig 7.
The obtained potentials and electron temperatures are a few times larger than those typically obtained by experiments [2, 7, 8, 37] because some actual geometries and physics are neglected in the model of this study as described in the last paragraph of Section 3.1.First, SEE may be an important factor, especially on the chamber downstream wall at z = 0.8 m, where electrons incident with an energy of 6.37 eV on average are observed in this study.In this energy region, according to Ref. [38], the total emission yield from carbon is between 0.1-0.5, which means more than 10% of electrons are recovered when they hit the wall, although this estimate requires a large extrapolation of their data to the much lower energies of our case.The small energy secondary electrons emitted from the grounded walls would be trapped inside the high potential plume if they are born with energies of 2 eV, another unknown and could reduce the plume potential.Second, since an actual GIT has holes on the exterior grid, neutralization is initiated closer to or inside the exterior grid by electrons inside the grid holes.The third possibility is the axial location of the neutralizer exit.According to Ref. [37], the angle between the normal beam axis and the neutralizer axis can affect the plume potential.Finally, another past study [39] suggested that the neutralizer-ion beam coupling was enhanced by ions generated in electron-neutral collisions, which is also not modeled in this study.

Neutralizer Position
This section examines how the neutralizer position affects the coupling between ion and electron sources in ground-based testing.Specifically, we consider only the type B neutralizer and perform calculations for the in-chamber case, where the neutralizer position is the only variable, i.e., the cases 2B-0, 2B-1/4, and 2B-1/2.The charge density in the x = 0 m plane is shown in Fig. 9 for different distances between the neutralizer exit and the plasma screen wall, d 0 .In all cases, electrons enter slightly above the ion source, form a negative charge density area, and then move toward the lower-right direction, as shown by the green arrows.Electrons begin to slow down and accumulate when they pass through (y, z) = (0.4,0.2) m, where the potential is at its maximum.A negative charge density area around (y, z) = (0.3, 0.35) m is formed on the opposite side of the electron source from the point of maximum potential, which we refer to as an "electron pool" in this study.This electron pool has also been observed in previous studies [19,22] and is a unique phenomenon in ion thruster plume neutralization with neutralizers adjacent to the thruster.
When d 0 /R 0 = 0, the electron pool is formed near the plume center, but as d 0 increases, the negative region shifts towards the lower right.However, these characteristics of electron motion are mainly observed on the thruster center plane, including the neutralizer exit (the x = 0 m plane).Figure 10 shows the charge density distribution in the x-y plane at the z = 0.45 m.Although there is significant variation in charge density around (x, y) = (0, 0.25) m or (0, 0.6) m in Fig. 10, the plume is uniformly distributed in most radial directions.
Next, Fig. 11a shows the normalized potential along the thruster axis for each neutralizer-type position.The plume potential increases as d 0 decreases, and the peak potential is particularly high for d 0 /R 0 = 0.This result indicates that the best neutralizer-ion beam coupling is obtained where 1/2 < d 0 /R 0 < 1 because the position of the maximum potential is seen to occur approximately (z − z 0 ) = R 0 , and there is no significant difference between the d 0 /R 0 = 1/2 and 1 cases.
Figure 11b shows the normalized potential on the neutralizer axis near the neutralizer exit.In all cases, the potential decreases just after the neutralizer exit and increases after the potential reaches a minimum value at approximately (z − z e0 )/R e = 0.2, forming a virtual cathode.The electron space charge limits the low-energy electron transport in the virtual cathode region, as explained in a previous experimental study [37].As d 0 increases, the minimum potential in Fig. 11b increases.This is due to the relaxation of the space-charge limitation as the exit of the neutralizer approaches the high potential space.In addition, when the neutralizer exit is on the wall (d 0 /R 0 = 0), there is no path for electrons to travel upstream to the neutralizer exit and into the plume, thereby reducing the neutralizer-ion beam coupling and increasing the electric potential in the plume.The results shown in this section conclude that moving the neutralizer downstream contributes more to neutralizing the ion beam.This trend is consistent with an experiment in which the coupling voltage decreased as the neutralizer moved toward the downstream [37].

Background Pressure
To understand the effect of background pressure in the vacuum chamber, we use the 2B case as the baseline.Figure 12 shows the background pressure, p, in the x = 0 m plane, where p, is calculated from the ideal gas equation of p = n n kT n , and n n and T n = 300 K are the neutral number density and neutral particle temperature, respectively.As a result of the neutral number density at the thruster exit and the size of the numerical pump, the highest pressure is about 4 µTorr at the thruster exit, with a minimum pressure of about 1.6 µTorr around the center of the vacuum chamber.This is similar to the typical background pressures found during ground test experiments.CEX and MEX collisions with the background neutrals change the ion velocity, which, in turn, alters the plume.CEX ions, originally neutral particles, have much smaller velocities than the beam ions and are the main cause of facility effects due to the high-background pressure.Figure 13 shows the number density distribution of CEX ions in the x = 0 m plane.CEX ions are produced in all regions where beam ions exist, but many CEX ions are particularly produced immediately downstream of the thruster, where the neutral number density is a maximum.However, highdensity areas also appear outside the plume core, indicating an asymmetric structure.This occurs because the CEX ions produced with a very small velocity remain in the electron pool region (y, z) = (0.3, 0.35) m, as indicated in Fig. 9. Similarly, CEX ions are trapped in the virtual cathode near the neutralizer exit.
Figure 14a shows the potential on the thruster axis, with the maximum potential decreasing in the presence of background neutrals.The decrease is due to two reasons.First, more electrons are present in the plume core.When collisions with neutral particles are modeled, the charge density on the thruster axis for ions and electrons in Fig. 14b indicates that the presence of CEX ion increases the ion density of the plume by up to 14%, while the electron density increased by nearly 35%.The second reason for the decrease in electric potential is the decrease in electron temperature.The electron temperature at (x, y, z) = (0, 0.4, 0.1625) m in the 2B-BP case obtained by a Maxwellian fitting (Eq.( 7)) of 15 eV is nearly half that obtained in the 2B case due to the increase in the minimum potential of the virtual cathode.The simulations show that the accumulation of CEX ions in the virtual cathode increases its minimum potential of the virtual cathode from φ = −3.2 to −2.7 V.The smaller decrease in voltage at the virtual cathode means that fewer electrons return to the neutralizer, and the kinetic energy of the electrons that can pass through the virtual cathode is lower, resulting in greater neutralization of the plume.Another interesting difference is the ion sheath formed on the side walls.Figure 15 shows the normalized volume charge density, ρ/ρ 0 , of the 2B and 2B-BP cases.Only when background neutral particles are present, the ρ/ρ 0 = 0 line appears between the chamber side wall and the plume core region.This is due to slow CEX ions that create an ion sheath in front of the chamber wall, where we define the ion sheath as the volume near the wall where ρ/ρ 0 > 0. The thickness of the ionic sheath is approximately 0.1 m.A 1-D analytical expression for the Child-Langmuir sheath thickness, s, can be calculated as [40] where α i = 0.61 is the number density factor [40]. λ D,s , φ s , and T e,s are the local Debye length, the electric potential, and the electron temperature at the sheath edge (ρ/ρ 0 = 0 line), respectively.Using (x, y, z) = (0.3, 0.4, 0.3) m as the reference point, we obtain from our simulations: n i = n e ∼ 1.0 × 10 11 m −3 , T e,s ∼ 2 eV, and φ s ∼ 7 V.Thus, the analytical sheath thickness is calculated as s ∼ 0.086 m, which is close to our simulated sheath thickness.Having demonstrated that not only background neutrals but also the difference between the space and ground chamber significantly influence electron transport, electric potential, and electron temperature, we next investigate electrical facility effects in terms of the current flow to each location in the simulation domain for cases 1B, 2B, and 2B-BP.Table 6 displays the currents for ions and electrons based on the current definitions shown in Fig. 3. Using the CCE BC (Case-1B),    Comparing the results with and without background neutral particles in Cases 2B vs. 2B-BP, we find that the current flowing to the plasma screen increased for both ions and electrons due to the backflow of CEX ions.Additionally, we demonstrated that the presence of CEX ions increases the potential of the virtual cathode, which results in a decrease in the current flowing back to the neutralizer.The current flow to the downstream and the side walls (I vc,end and I vc,side , respectively) exhibits an interesting behavior due to the presence of the background neutral particles.With respect to ions, 24% of the ion currents change their direction toward the side wall (see the difference in I i,vc,side ) due to collisions with neutral particles, whereas only 3% of the electron currents (see the difference in I e,vc,side ) change their direction.Consequently, ion-neutral particle collisions have little effect on the destination of electrons.However, note that this study ignores electron-neutral particle collisions, which may affect the electron current.2. b See Fig. 3 for definitions of the current to each part.

Sensitivity of Ion Plume Due to Thruster and Neutralizer Potential
This section investigates the changes in the plume due to different electric potential BCs, which can vary with each thruster potential condition.Figure 16 shows the calculated potentials for the four cases shown in Table 4.As can be seen from this figure, the differences significantly affect the plume.A detailed discussion follows.

Accel Grid Potential
First, we present the results of changing the potential at the thruster exit (V th ), simulating the case when a negatively biased accel grid is exposed to the plasma without a decel grid Case 2B-ACC.As shown in Fig. 16, low potentials are observed near the negatively biased outlet, while a higher potential region is seen downstream.In Fig. 17, the potential plotted on the thruster axis shows that in the 2B-ACC case, it is initially -200 V but reaches a larger maximum potential than that in the 2B-BP case at z ∼ 0.2 m or 0.05 m downstream from the peak of the 2B-BP case.The reason for the larger maximum plume potential despite the lower potential of the thruster may be attributed to the degradation of the coupling between the neutralizer and the thruster.In the 2B-ACC case, the ion inflow velocity increases by an amount proportional to the square root of the accel grid potential of 200 V because there is no deceleration by a decel grid (2B-ACC).The axial ion velocity, w i , in the y = 0.4 m plane is shown in Fig. 18a, where w i is normalized by the reference beam ion velocity of w i0 = 30, 0000 m/s (612 eV).While the 2B-ACC case shows larger velocities near the thruster exit, there is almost no difference between the two cases downstream.The line plots in Fig. 18b also show that the velocities are almost identical, but the beam width is narrower in the 2B-ACC case.This is due to the convergence of the plume near the thruster caused by the external electric field induced by the potentials between the plasma screen of 0 V and the thruster exit of -200 V.This was also observed in an experiment (Ref.[41]), where the grounded thruster cover around the accel grid reduced the beam divergence angle.

Neutralizer Potential
Next, we investigate the cases where the neutralizer exit potential (V ne ) is 5 V lower (2B-NM) and 5 V higher (2B-NP) with respect to the ground voltage of 0 V of the 2B-BP case.As shown in Fig 16, the potential of the plasma plume changes significantly even though the exit potential is only changed by ±5 V. Figure 19a shows the potential along the thruster axis where it can be seen that the maximum potential is increased by eφ/kT e0 ∼ 50 for the 2B-NM case, while for 2B-NP it is decreased by eφ/kT e0 ∼ 50.Additionally, as shown in Fig 19b, a comparison of the potentials near the neutralizer indicates that as the neutralizer potential is lowered, the downstream potential well becomes smaller, and it disappears completely in the 2B-NM case.The reason that this occurs is related to the current distributions, as we discuss next.Table 6 shows the currents to different parts of the GIT for three cases (2B-BP, 2B-NM, and 2B-NP).Regarding ion currents, the large plume potential for Case 2B-NM shown in Fig 16 causes more divergence of the ion beam, and a larger number of CEX ions return to the chamber side walls and the thruster, compared to case 2B-BP.However, a more significant difference is observed in the electron currents than in ion currents.As shown in Fig 19b, when V ne is small, all electrons are allowed to flow downstream, while when V ne is large, more electrons return to the neutralizer, as a result of virtual cathode formation.Here, I e0 − I e,ne can be considered an effective emission current from the neutralizer, I e,eff , which increases as V ne decreases.This trend is consistent with experimental results showing that lowering the coupling voltage of the neutralizer releases more electron current [37].Furthermore, there is a large difference in I vc,end , but a small difference in I vc,side because electrons do not return to the neutralizer.
Despite the increase in I e,eff , the plume has a very high potential in the 2B-NM case.The electric potential where electrons are emitted with respect to the vacuum chamber, V ne , changes the electron motion in the vacuum chamber.Electrons produced at lower potentials than the vacuum chamber (V ne = −5 V, 2B-NM) do not return to the plume because of the sheath in front of the vacuum chamber walls and are all absorbed by the wall, resulting in a shortage of electrons for neutralization.On the other hand, the sheath reflects almost all electrons produced at higher potentials (V ne = 5 V) except for those with high kinetic energy.As a result, a sufficient number of electrons remain in the chamber when V ne is large (2B-NP), as shown in Fig. 20.In Table 6, I vc decreases as V ne increases, which confirms this electron confinement effect inside the vacuum chamber from the viewpoint of the currents.Figure 21 shows the EVDFs in the z-directions obtained by sampling computational electrons at R 0 downstream on the thruster axis for the 2B-BP, 2B-NM, and 2B-NP cases.Two distributions are clearly seen in the cases 2B-BP and 2B-NM, but almost all electrons are thermalized in the 2B-NP case.Fitting the electron temperature to the distribution for the population centered around zero velocity reveals that T e decreases as V ne increases.This is because the plume potential is reduced by the above-mentioned change in electron motion (see Figs. 16 and 19a), allowing the otherwise-trapped electrons to move with less energy.

Conclusions
In this article, we have simulated ion thruster plumes in a vacuum chamber and space configurations using the in-house multi-GPU CHAOS solver.A fully kinetic PIC-DSMC approach is applied to model electron motion emitted from two types of external neutralizers separated from the thruster.This study qualitatively confirms important aspects of electrical facility effects and improves our understanding of how electrons move in a vacuum test chamber for practical GIT configurations.
First, we have shown that the plume potential and electron temperature are larger in the ground test than in the space environment because the former absorbs all the electrons, reducing the number of electrons trapped in the ion beam potential.This causes a large change in the electron current.The backflow current from the plume to the propulsion system, including the neutralizer, obtained under vacuum chamber conditions is less than half that under space conditions, indicating that most of the current flows to the downstream wall.Another facility effect that we investigated  ne = 0 V, (2B-BP)  ne = -5 V, (2B-NM)  ne = 5 V, (2B-NP) Figure 21: EVDF in z-direction for the 2B-BP, 2B-NM, and 2B-NP cases.Computational electrons are sampled at R 0 downstream from the thruster exit center.The solid lines are Maxwellian distribution fitting curves, based on Eq (7), and the normalized electron temperature obtained by fitting, T e /T e0 , are 15 eV (2B-BP), 34 eV (2B-NM), and 2 eV (2B-NP), where T e0 = 2 eV. is how background pressure affects the plume.In particular, it was shown that the electron flow is also changed when slow CEX ions accumulate in areas with negative charge density.This causes the maximum potential of the plume and electron temperature to decrease.As for the ion current, the reverse current to the thruster and the current to the side walls of the chamber increase significantly.The electron current also changes but at a smaller rate than the ion current.
For further practical insight, we investigated the coupling between the ion beam and electrons by changing the position of the neutralizer while keeping the emitted electron current constant.As the neutralizer is moved downstream from the thruster wall, the maximum potential of the plume decreases, and the minimum potential of the virtual cathode created at the neutralizer exit increases.Electrons emitted from the neutralizer flow downstream in the central plane of the thruster, moving back and forth across the center of the plume where the potential is high, creating an asymmetric charge density distribution.
In addition, we have also simulated different electric potentials at the ion and electron source exits.In the case of no decel grid, where the -200 V accel grid is exposed to the plasma, the maximum potential increased due to the difference in coupling with the electron source, regardless of the lower exit potential.From a performance standpoint, the ion beam is more focused, but the final ion velocity remains the same.When the neutralizer potential was lower than the chamber wall, the plume potential increased significantly because more electrons were extracted and absorbed toward the high potential chamber wall, resulting in insufficient neutralization.Conversely, when the potential of the neutralizer is higher than the chamber, the chamber acts as a cage for electrons, and the plume is more neutralized.

Figure 1 :
Figure 1: Computational domain setups for Type-A electron source cases, where the thruster exit radius, R 0 , and the neutralizer exit radius, R e0 , are the same as R 0 = R e0 = 6.25 cm.

Figure 2 :
Figure 2: Computational domain setups for Type-B electron source cases, where the neutralizer exit radius, R e0 = 1 cm, is smaller than the thruster exit radius, R 0 = 6.25 cm.The distance between the neutralizer exit and thruster body, d 0 , changes for four cases as d 0 = 0, 1 4 R 0 , 1 2 R 0 or, R 0 .
Figure3: Electrical schematics of this study.I i0 and I e0 are the currents of ion emitted from the thruster and electrons emitted from the neutralizer, respectively.I ne is the current to the neutralizer exit, I th , I ps , and I vc are the currents to the thruster exit, plasma screen, and vacuum chamber, respectively.V dc , V th , and V ne are the biased voltages of the discharge chamber, thruster external gird, and neutralizer exit, respectively.

Figure 5 :
Figure 5: Volume charge density in the x = 0 m plane for the in-space vs. in-chamber cases for neutralizer types A and B, respectively.The green lines show electron streamlines based on electron z-and y-velocities.

Figure 6 :
Figure 6: Electric potential in the y = 0.4 m plane for the in-space vs. in-chamber cases for neutralizer type A and type B, respectively.Here, and in similar subsequent figures, kT e0 =2 eV.

2 )Figure 10 :
Figure 10: Volume charge density in the x-y plane for different d 0 /R 0 values of the type B neutralizer in the vacuum chamber in the z = 0.45 m plane.

Figure 11 :
Figure 11: Electric potential along the z-direction for different neutralizer locations for the neutralizer Type: B in the vacuum chamber.

Figure 12 :
Figure 12: Background pressure distribution for the 2B-BP case in the x = 0 m plane.

Figure 13 :
Figure 13: Number density of CEX ions for the 2B-BP case in the x = 0 m plane.

Figure 14 :
Figure 14: Electric potential and volume charge density along the thruster axis for w/ vs. w/o background neutral for the type B neutralizer in the vacuum chamber.

Figure 16 :
Figure 16: Electric potential in the x = 0 m plane for various electric potential BCs.

Figure 17 :
Figure 17: Electric potential along the thruster axis for different thruster exit potential BCs.

Figure 19 :
Figure 19: Electric potential in z-direction for different neutralizer exit potential BCs.

Figure 20 :
Figure 20: Volume charge density along the thruster axis for different neutralizer exit potential BCs.

Table 2 :
Parameters of the species at the thruster and neutralizer exits.

Table 3 :
Test conditions of each case ID for facility effect and neutralizer position study.The neutralizer types, A and B, are described in Table2. *

Table 4 :
Test conditions of each case ID for electric potential boundary study.V th Neutralizer exit, V ne

Table 5 :
Number of computational particles in the entire computational domain at the steady state.

Table 6 :
Ion and electron current from thruster plume to different locations.ps I vc (= I vc,end +I vc,side ) I ne I th I ps I vc (= I vc,end +I vc,side ) All current values are normalized by the emitted ion or electron current (I i0 , I e0 ), given in Table a