Electrostatic Ion Thrusters - Towards Predictive Modeling

The development of electrostatic ion thrusters so far has mainly been based on empirical and qualitative know-how, and on evolutionary iteration steps. This resulted in considerable effort regarding prototype design, construction and testing and therefore in signiﬁcant development and qualiﬁcation costs and high time demands. For future developments it is anticipated to implement simulation tools which allow for quantitative prediction of ion thruster performance, long-term behavior and space craft interaction prior to hardware design and construction. Based on integrated numerical models combining self-consistent kinetic plasma models with plasma-wall interaction modules a new quality in the description of electrostatic thrusters can be reached. These open the perspective for predictive modeling in this ﬁeld. This paper reviews the application of a set of predictive numerical modeling tools on an ion thruster model of the HEMP-T (High Efﬁciency Multi-stage Plasma Thruster) type patented by Thales Electron Devices GmbH.


Introduction
Ion thrusters, where the propellant is ionized and the ions are accelerated by electric fields, are of increasing importance for scientific and commercial space missions. Compared to the commonly used chemical thrusters they have a 5 to 10 times higher specific impulse [1]. This results in a considerable reduced propellant budget such that a significant reduction of spacecraft launch mass by some 100 to 1000 kg can be achieved. The most straightforward concept of an electric propulsion system is the grid thruster. Here, gaseous propellant is ionized in a low pressure, low temperature gas discharge sustained by DC, radio-frequency or microwave fields, respectively. A biased grid at the thruster exit accelerates and extracts the produced ions. Since ionization and acceleration occurs independent of each other, it is possible to optimize these two processes separately. This yields a high efficiency of the device but also a huge system complexity. Therefore ion propulsion systems based on grid thrusters tend to be expensive and to exhibit reliability issues [2]. An alternative are grid-less ion thrusters which are based on magnetic confinement of the plasma electrons, where the trapped electrons both ionize the propellant and provide the potential drop for ion acceleration. Due to their low complexity in terms of system architecture they are becoming of increasing interest in particular for commercial satellites. In order to achieve reduced development and qualification effort, it is therefore needed to set up and apply a series of different modeling tools which can quantitatively describe the plasma physics within the thruster but also interactions of the thruster with the testing environment and finally the satellite. The integrated modeling strategy should include several modular components in a consistent way in order to provide the necessary complexity and accuracy depending on the problem. As an example, results for the grid-less High Efficient Multi-stage Plasma Thruster (HEMP-T) model DM3a are presented. The outline of the paper is the following: after a general introduction of the HEMP-T the general strategy for predictive modeling is presented. The characteristics of HEMP-T as deduced from 2D Particle-in-Cell (PIC) calculations are reviewed. A first example of 3D PIC results for self-consistent modeling of electrostatic turbulence in a HEMP-like geometry is given. Strategies and results for the simulation of the interaction of the thruster with the satellite in the plume and results on erosion are discussed afterwards, before the paper is summarized in the conclusions.

The High Efficient Multi-stage Plasma Thruster (HEMP-T)
An ion thruster design which has the advantage of an electric field topology similar to grid thrusters at significant reduced thruster complexity is the High Efficient Multi-stage Plasma Thruster (HEMP), patented by the THALES group in 1998 [3]. As described in [1], HEMP-T's consist of a dielectric, rotationally symmetric discharge channel at the upstream end of which an anode is located. The anode is connected to the power supply and represents the only high voltage electrode in the thruster; it also serves as inlet for the propellant. The discharge channel is surrounded by a system of axially magnetized permanent magnet rings in opposite magnetization, the so-called PPM (Permanent Periodic Magnet) system. At the downstream end of the discharge channel, the thruster exit, a hollow cathode neutralizer is placed to provide the starter electrons for igniting the discharge and for neutralizing the ion beam emitted by the thruster. A schematic picture of the HEMP-T concept as described above is given in Fig. 1. The PPM system forms a linear magnetic multi-cusp structure inside the discharge channel. The magnetic field periodically changes between predominant axial and predominant radial direction. The level of magnetic induction B at any position within the thruster channel is chosen such that the Larmor radius of the electrons is much smaller than the geometrical dimensions of the discharge channel, whilst the propellant ions are hardly affected by the magnetic field due to their much higher mass. The gradients fulfill dBz dz > 0 in the axial zones and dBr dr > 0 in the radial zones, where in addition a strong gradient dBr dz is build up. In this way the plasma electrons are efficiently confined along the entire discharge channel and only few electrons are lost to the wall mostly at the cusps. Electron confinement due to cusp mirror oscillations is dominant compared to the Hall current which builds up due to the E × B drift. In this work an older prototype model named DM3a is discussed. It has implemented two cusps, the anode and the exit cusp and is described in [4]. In the HEMP thruster the ionization is particularly strong at the cusps regions. Due to the strongly reduced plasma wall contact in the cusp regions and mean energy of impinging ions at these locations below the sputter threshold, HEMP is characterized by a long lifetime of at least 10.000 h [5].
Potential, densities and plasma species temperatures of the DM3a were calculated with an electrostatic 2d3v Particle in Cell code with Monte Carlo collisions (PIC MCC) [6]. PIC MCC is a simulation method, used for low temperature plasma's [7], [8]. It gives a fully self-consistent microscopic description of a plasma and is able to involve complicated atomic and plasma-surface interactions. In the PIC-MCC simulation we follow the kinetics of so-called Super Particles (each of them representing many real particles), moving in the self consistent electric field calculated on a spatial grid from the Poisson equation. The particle collisions are handled by Monte-Carlo collision (MCC) routines, which randomly change particle velocities according to the actual collision dynamics. All relevant collisional processes are included in the model: electron-neutral elastic, ionization and excitation collisions, ion neutral momentum-transfer and charge exchange collisions. The dynamics of the background neutral gas is self-consistently resolved with direct simulation Monte Carlo. For a reliable plasma simulation, effects on the smallest length scales in a plasma, the Debye-scale λ D,e have to be resolved. Therefore, the smallest lateral grid size is Δx = 0.5λ D,e . The simulation has to resolve the fastest process in the system -the electron Langmuir oscillations. Therefore, the time step in the simulation is chosen dt = 0.2/ω pe , where ω pe is the electron plasma frequency. This dt guarantees a sufficient resolution of all time scales of the system. A more detailed description of PIC can be found in [8]. An important ingredient for speeding up the 2D PIC calculations is the use of a similarity scaling [6]. A typical value of the scaling factor is 0.1. Scaling is particularly well applicable for the HEMP-T concept because of the negligible wall effects due to the efficient magnetic plasma confinement. In Fig. 2 the potential as calculated with an electrostatic 2D PIC code is shown. For the HEMP thruster the potential in the plasma bulk is nearly constant with a steep drop at the thruster exit producing a narrow peak in the ion energy spectrum a high specific impulse. Close to the axis, the mainly axial magnetic field allows the electrons to flow along the electric field. A small perturbation of the electric potential is therefore quickly compensated by fast electrons. The inner surface of the channel walls are made of Boron Nitride based ceramics. Near the thruster walls, the potential is decreasing forming the electrostatic sheath. The maximum potential drops are located at the cusps, where magnetic field is perpendicular to the wall.
The influence of the plasma sheath can also be seen in the electron density in Fig. 3. The electron density is uniformly decreasing towards the channel walls and the cusp regions are clearly visible. The higher electron density in these regions due to electrons confined at the radial magnetic field at cusps can be seen. The transfer of energy from directed into thermal motion heats the electrons [9]. Thus, the density of the electrons in the cusp regions is high and the ionization can take place efficiently. For the HEMP thruster the electron source provides the primary electrons, which defines the operational point for the thruster due to the very strong amplification of it by ionization inside the thruster. The direct connection of anode and exit region close to the axis allows ignition of the thruster even without external source by some free electrons. These electrons can be created e.g. by cosmic radiation and will be accelerated to anode potential starting the ionization avalanche [10]. To characterize the basic physics of HEMP-T the velocity distribution functions resolved spatially along the thruster axis and averaged in radial direction were calculated. All velocity distribution functions (Fig. 5, 7 and 8) are temporally averaged over 10 6 time steps of a quasi steady-state run.
In Fig. 5 the axial velocity distribution function of electrons is plotted. Electron energy distribution functions at different locations are shown in Fig. 6. The distribution functions in the channel are non-Maxwellian, because the mean free paths for electron Coulomb collisions is much longer than the electron path in the system. The axial velocity distribution function shows that the electron axial energy is much higher in the regions about z= 20 mm www.cpp-journal.org and z= 51 mm, since electrons are heated at the cusps. Outside the cusps the velocity distributions are nearly identical within the channel.    The ion density profile is and the corresponding distribution function of the axial velocity are shown in Fig. 4 and 7. Within the channel all ions have rather low velocities due to the nearly constant potential profile. Only at the exit they are accelerated by the strong potential drop. Near the exit cusp and in the plume ions with high energies are existing. The first peak in the distribution at z= 48 mm with velocities up to 20000 m/s can be explained by the grounded wall after the dielectrics. The right side of the potential structure at this position is very steep and creates a region with large positive velocities. The boundary of this structure is a direct consequence of the potential structure in this region as shown in Fig. 2, because after the grounded wall the potential rises again and ions are slowed down by the counter-acting electric field. Close to the exit cusp inside the channel there are many ions with low velocities. At the exit cusp and afterwards they are accelerated, seen by the shift of the orange color to higher velocities. This acceleration is a consequence of the steep potential drop at the exit. In this acceleration region the ion density decreases. Most of them escape from the region due to their high velocity in radial direction. Further ions are accelerated until the grounded end of the computational domain, due to the continuous drop of the potential. Some of them reach a velocity of about 18000 m/s. This explains the trend to higher velocities in Fig. 7 beginning at the exit cusp at 51 mm up to the end of the domain. This structure looks like a quarter of a circle in the velocity distribution function. This is close to a motion with constant acceleration, where the velocity is proportional to the square-root of the distance. Ions are also scattered into the side part of the plume by charge-exchange collisions, where so-called wings can build up.
In Fig. 8 the axial velocity distribution function for neutrals is presented. One can see that charge exchange collisions have an impact on the velocity profile of neutrals. They couple the neutrals to the ions and are responsible for similarities in both distribution functions. Apart from the high-velocity wings the distribution looks quite more regular, since the neutrals are not affected by the magnetic or electric field. Most of the neutrals in the channel have velocities in the range of -2000 m/s . . . 2000 m/s, equally distributed in positive and in negative direction. This is again a result of the large neutral density there and their relative low temperature. In the plume the situation for neutrals is different. Neutrals are only flying out of the thruster and negative velocities are not detected.

Modeling strategy
The most complete model resolving all time scales would be a direct coupling of a kinetic plasma model with a molecular dynamics model for the walls. This would allow a fully self-consistent analysis of the complete system including plasma dynamics, possible erosion of thruster walls and interaction of the exhausted ions with surrounding satellite surfaces or, during test and qualification, with the testing environment, like the vacuum chamber walls, respectively. Such a solution is not possible due to the tremendous computational costs and high complexity of this combined model. Instead, we propose to use a hierarchical multi-scale set of models in which the parameterization for a lower hierarchy model can be deduced from a higher one (see Fig. 9). Fig. 9 The concept of a multi-scale modeling for combined thruster-plume models.
For example, a 3D particle-in-cell (PIC) model can deliver a parameterization of turbulence effects by appropriate anomalous transport coefficients. A first proof-of-principle example of a 3D PIC code with a simplified geometry will be presented in section 4. Transport coefficients based on such runs could then be used in a 2D PIC, which is more practical for production runs.
Furthermore, separating the analysis of the thruster from the plume is not adequate, because there exist a strong coupling between thruster and plume plasma. Due to strongly non-Maxwellian characteristics of the distribution www.cpp-journal.org functions of both electrons and ions it is not possible to keep any of them in a fluid representation. To get a correct description of both thruster and plume plasma one has to solve a kinetic problem for the whole region of interest including all significant physical processes. These are collisions, turbulence effects, surface driven sheath instabilities and breathing modes. That is why the choice of a PIC model is a natural one for such a problem. In addition, a similarity scaling is applied to further reduce the calculational costs [6].
However, despite the scaling, it is not possible to have a uniform grid for the whole region (thruster and plume) due to extremely different plasma conditions and also physical dimensions of both parts. Therefore, using an equidistant grid will lead to unacceptable grid sizes in terms of computational costs. A possible strategy for the plasma model is a 2D PIC model with a hierarchical grid (see section 5).
In order to describe erosion-redeposition processes one can use different approximation levels of the model. The most thorough description is given by the full molecular dynamic model, which would be far too timeconsuming because it resolves all individual atoms and their interactions. Next level can be represented by the binary collision cascade model assuming an amorphous target and the interaction of particles with the solid based on heavy particle collisions with ions and additional losses with electrons acting as viscous force. Such a model can use the detailed information about flux distributions provided by the PIC code and calculate then based on this the erosion response of the materials. The most crude approximation is given by a Monte-Carlo (MC) procedure simulating erosion-redeposition based on sputter yield tables calculated from the binary-collision cascade or molecular dynamics model together with the information about plasma fluxes. This model is particular useful due to its simplicity and flexibility allowing to quantify the life time of ion thrusters (see section 6).

Anomalous diffusion of electrons
In Hall effect thrusters a radial magnetic field is applied to confine the electrons in the thruster channel and to create the accelerating potential for the ions. The axial electron current towards the anode is possible only because of presence of some diffusion process. There are two major diffusion processes responsible for the electron current across the magnetic field. Classical diffusion due to collisions with neutrals and anomalous diffusion due to azimuthal fluctuations of the electric field.
In classical diffusion, the diffusion coefficient D is proportional to the square of the mean free path λ mf p and one over the average time τ between two collisions Due to the gyration of the electrons, the mean free path in a plasma with an external magnetic field is its Larmor radius r L,e = | v ⊥,e | ωc,e . For a strong magnetic field with ω c,e τ >> 1, the diffusion coefficient due to neutral scattering is given in [11] as While the Larmor radius is proportional to 1/B, the diffusion across the magnetic field scales with 1/B 2 . While for diffusion along magnetic field lines collisions are decreasing the mean free path and thus the collision coefficient, here collisions are necessary, therefore D ⊥ is proportional to the collision frequency ν = 1/τ . But in many experiments on the magnetized plasma's, the contribution of classical diffusion is not sufficient to explain the measured electron transport across the magnetic field lines [12]. Another kind of diffusion is the so-called "anomalous" or Bohm diffusion, which is due to electron scattering on the fluctuations of the perpendicular electric field. While in many E ×B experiments D ⊥ scales with B −1 rather than with B −2 , a semi empirical formula was proposed by D. Bohm, E. Burhop, and H. Massey, with a diffusion coefficient [13]: A general derivation for D ⊥ ∝ B −1 for this kind of diffusion was given by L. Spitzer in 1960 [12], with K 1 , K 2 , K 3 are empirical constants, determining the strength of the electron perturbation perpendicular to the magnetic field.
To study anomalous diffusion in the HEMP thruster, a three dimensional electrostatic PIC code is used. Because the calculation is extremely time consuming, only a proof-of-principle for studying anomalous transport in a HEMP-like geometry was possible. The HEMP thruster has a cylindrical volume. But in a PIC simulation a cylindrical grid causes self forces. To avoid that, a Cartesian grid is chosen. The HEMP model is set up as a cuboid, based on a Cartesian grid. To be able to obtain a solution in a reasonable run time (one week), the size of the system is scaled down by a factor of 50. In order to preserve the ratio of the charged particles mean free paths and the gyro-radii to the system length, the neutral Xenon density and the magnetic field are increased by the same factor 50 [6]. The simulation domain has a height and width of L x = L y = 23 mm, a length of L z = 130 mm and includes the thruster channel as well as the plume region, as shown in Fig. 10. The thruster exit is positioned at z = 51 mm and the two cusps at z = 20.5 mm (anode cusp) and z = 37 mm (exit cusp). The channel walls consist of boron nitride based ceramics. The external domain boundaries consist of metal walls with a potential of φ = 0. The anode voltage, applied at z=0 mm is set to 500 V. Particles which are hitting the metal as well as the dielectric boundaries are absorbed. In the thruster channel, these particles are contributing to the local surface charge. To identify anomalous diffusion driven by turbulence, no secondary electron emission was implemented. For the neutral xenon atoms, the density profile was prescribed by an exponential profile n n = n n exp{−Z/L}, with L = Z max /2.5 and Z max = 65 mm. Such density profile is close to the one obtained from self consistent models in 2D. In the following, the three dimensional data is shown in two dimensional cuts. In Fig. 11 a profile of the plasma potential along the y-axis, for y = L y /2, is shown. One can see that the potential has a step-like shape in the z-direction, which is caused by the two cusps. The potential shape is in overall agreement with the 2D results: the main potential drop occurs at the exit cusp, as shown in Fig. 2. At z = 20.5 mm, the potential gets negative at the wall. This is the anode cusp region, where the magnetic field lines are perpendicular to the wall. The second cusp region can be identified by a potential plateau in the magnetic bottle between z = 23-37 mm. In Fig. 12 the potential in this cusp is shown in the x-y-plane. The variation of the potential is not purely radial. One can identify a n=1 azimuthal mode at about R= 3 mm and a n=2 mode at about R= 7 mm with an amplitude of electric field E ≈ 100 V/cm. These fluctuations of the azimuthal electric field are responsible for the anomalous electron transport across the cusps in the simulation. One should note that this mode appears also for larger domains and seems not to be a consequence of the non-cylindrical geometry.
The electron density is shown in Fig. 13 and 14. The density profile for y = L y /2, Fig. 13, shows a large electron density at the anode (z = 0) and at the anode cusp (z = 20.5). The density drop between the two regions is due to the electron transport along magnetic field lines. At the exit cusp, practically no electrons can be found, because it is not yet filled with charged plasma particles due to run-time limits of the code. Due to that, no negative wall potential for the exit cusp can be seen in Figure 11. Estimation of the anomalous diffusion coefficient, caused by the fluctuations of the azimuthal electric field following the approach suggested in [12] gives which is in the order of magnitude of the diffusion coefficient given by Equ. (3). One has to mention, that the presented calculation is only a rough estimate, because of the large scaling factor used in the simulation. Due to the geometrical characteristics of the HEMP-T with large parts where the magnetic field is essentially axial and only small zones of radial fields at the cusps, the anomalous transport is only important at the cusps.
Here, the anomalous diffusion allows the electrons to overcome the radial magnetic field and to fill the region between cusps. In the region of axial fields the parallel transport along field-lines dominates. This is different to Hall thrusters, where in the whole acceleration channel the radial magnetic field is applied and the thruster characteristics is determined by anomalous electron transport. As a consequence modeling optimization with a 2D code of HEMP-T is possible, because classical transport dominates.

Plume modeling
Up to now, either fluid or hybrid models (with kinetic ions) are used for the plume modeling. However, the plasma in the acceleration channel is already non-Maxwellian [14]. The mean free paths of particles in the plume of about ∼10 m does not allow a relaxation of the distributions to a Maxwellian in the plume region. Therefore, a correct model of the plume plasma has to be kinetic. Unfortunately, the well-established Particle-In Cell (PIC) method is momentum conserving and free of artificial self forces only in the case of equidistant meshes [8]. This means that a self-consistent PIC model for an ion thruster including the full plume plasma has to resolve the smallest length scale (usually the Debye-scale at the largest electron density inside the thruster) and the fastest time scale (usually the plasma frequency at the largest electron density inside the thruster). A typical thruster size is about 10 cm whereas the typical size of interest for a plume is 1 m∼10 m. Electron density inside the thruster is about ∼ 10 13 1/cm 3 and for such a density λ D = 7.43·10 −4 cm for T e = 10 eV. In the plume plasma density drops exponentially, such that on a distance of ∼5 cm from a thruster nozzle it is around ∼ 10 10 1/cm 3 which results in λ D ∼ 10 −2 cm. This allows only small domains for the plume due to computational time restrictions [6]. Such small domains suffer from the strong influence of the solution from the boundary conditions at the end of the plume, which can lead to results that are not consistent with experimental data, especially in terms of angular ion distribution functions. Fig. 15 demonstrates the influence of boundary conditions on a potential solution in the near-plume region. In addition, non-linear coupling of a plasma inside the thruster with the plume part may influence the operational regime of the thruster, including currents, density distributions and thrust. The presented method uses a "matryoshka-like" hierarchy of equidistant grids with different cell sizes. The hierarchy is constructed such, that the most dense grid covers only the thruster and the near-plume region; the next level with a cell size of 2 ∼ 4 times larger extends further to capture more of the plume and so on. Thus, the thruster and the near-plume region is covered with all grids and the most distant from the thruster regions are covered only by the most coarse grid. Furthermore, the density of the charged particles is gathered in all grids independently. The Poisson equation is solved one by one starting from the most coarse and finishing the most dense grid such, that the boundary values for the next level are taken from the coarser mesh obtained on the previous step. Such an approach appears to be not only accurate enough, but also remarkably fast compared with the solution for a single non-equidistant mesh. Such a fast solver is possible due to two reasons. First, the matrix which one gets after the discretization of the Poisson equation for an equidistant grid has a block structure. Such a structuring can be used by specialized solvers, for example [15] or [16]. And second, for such a single mesh one has a larger matrix to solve, rather then for the case of "matryoshka-like" grids. However, non-equidistant grids suffer from artifacts [8], like self-forces, and corrections to minimize such errors are needed. To overcome this, a modified two point central difference scheme for calculation of the electric field on non-equidistant grids is implemented [17].

Effects due to impinging ions
Ions created in the thruster discharge may impinge surrounding surfaces which can induce sputter erosion and redeposition of eroded material. Depending on the surface region this may affect operational and performance characteristics of the thruster itself, of the ion thruster module or even of the whole satellite, respectively. For the simulation one can distinguish: a) Impact on inner thruster surface by ions generated in the inner thruster discharge. b) Impact on exit-sided surface of the thruster and the neutralizing electron source by ions generated in the plasma plume downstream the thruster exit.
c) Impact on satellite surface producing erosion and redeposition. d) Impact on vacuum chamber walls during testing and life-time qualification creating redeposition onto thruster and thruster module surface.
The proposed multi-scale modeling strategy is well suited to address the above mentioned ion impingement effects. A good approximation of a plasma surface interaction are binary collisions between the impinging ion, which gets neutralized next to the surface, and the target atoms, producing collisional cascades in the solid. When a part of the collided particles get enough energy to leave the surface, the target emits them as sputtered particles. Sputtered particles are impurities in a plasma, values of sputter yields are important for plasma experiments and simulations. A tool for simulating binary collisions in matter is the SD.Trim.SP (Stationary/Dynamic Transport of Ions in Matter, with the calculation mode Serial or Parallel) code. The SD.Trim.SP computer program simulates sputtering, backscattering and transmission effects of ion bombarded material and can additionally take the modification of the target into account, when it runs in the dynamic mode. It applies the Monte-Carlo Binary Collision Approximation (BCA) and assumes therefore an amorphous (randomized) material with a infinite lattice size and a temperature of OK. In SD.Trim.SP the particle movement in matter is approximated as a series of inelastic binary collisions between atoms, the BCA and a continuous friction, to simulate the interactions of moving atoms with electrons. For additional information about the use of SD.Trim.SP, see [18]. The domain of SD.Trim.SP, is a one dimensional simulation space, where the Cartesian x-component is perpendicular to the surface. Also two dimensional simulations are possible A negative x-component indicates the space above the surface, while a positive one shows the position in the solid. Also layers of different materials can be implemented. SD.Trim.SP in static mode proceeds in the following way. At first a projectile is initialized with the kinetic energy E 0 and the direction r 0 . After a distance of λ, a collision partner is determined by the stochastic choice of an impact parameter p. While SD.Trim.SP assumes an amorphous structure of the material, no lattice structure has to be taken in account and therefore λ and p are determined by their distribution functions given by the BCA. Both are implemented with inverse Monte-Carlo sampling. The azimuthal angles between two collisions are chosen randomly between [0; 2π]. The BCA gives the energies of the particles after the collision and the scattering angle ϑ 1 as well as the recoil angel ϑ 2 , which are determining the new direction of the projectile and the target atom. The energy loss of atoms traveling through matter, due to interactions with electrons, is simulated as a continuous friction in between two collisions. Three scenarios are possible for each particle. If the energy is smaller than the binding energy of the matter E < E b the particle sticks and is not followed any more. If E > E b and the particle is close enough to the surface, it gets emitted as a sputtered atom and is also not followed any more. In the third case the particle moves through the matter and produces a collision cascade through several collisions, proceeded as described above. Reflection at the surface is realized with different binding energies for particles coming from inside or outside the target. To determine the dynamics of the target thickness, SD.Trim.SP has a dynamic mode. Here, the material is resolved one dimensional and the target is segmented into slabs. These slabs have an initial thickness, which changes during the calculation due to collisional transport. For many particles, the calculation as well as the memory occupation of every collisional cascade becomes very costly. Therefore, for large fluency pseudo particles which are representing a number of real particles are introduced to minimize the numerical costs. For an entire dose of Φ 0 the material should be exposed with, pseudo particles with a differential fluency of ΔΦ = Φ 0 /N d are followed in N d simulation steps. Moreover, the physical sputter yield rapidly decreases for energies lower than a threshold energy E thr .
Although the coupling of the BCA and the PIC models is promising in terms of analyzing erosion of a thruster during its operation, it is in practice inapplicable for other above mentioned tasks (c and d). That is why a Monte-Carlo model, which uses sputter yield tables pre-calculated with a binary collision cascade model is developed [18]. Coupling the plasma model with this erosion module an integrated model can be set up. The plasma fluxes impinging on the walls from PIC are used in a Monte-Carlo procedure for erosion redeposition simulations, where the erosion fluxes are determined from the tables. Application of this model is discussed in [19].
A terrestrial qualification of a thruster has a significant difference from outer space exploitation in that it is held in a limited vessel, which can create different artifacts on the measured thruster properties. For example, the back scattered flux from vessel walls can be deposited on the walls of the thrusters and by that create a conducting layer influencing the thruster operating regimes. The quantitative characterization of such an influence is possible by means of a self-consistent coupling between the PIC code modeling of the plasma and the Monte-Carlo (MC) www.cpp-journal.org erosion-deposition code modeling of erosion of the thruster walls due to plasma-wall interaction and of deposition of the eroded particles both from vessel and thruster walls. Due to large size difference of a thruster and a vessel it is possible to parameterize the back-scattered flux from vessel walls as an effective source for the MC erosiondeposition code. The primary distribution of ions with respect to energies, angles and species is specified and pseudo-particles are followed interacting with the vessel walls. Hitting a wall, based on sputter rates calculated by a binary collision cascade code, the back-flow of eroded particles from vessel walls towards the ion thruster acceleration channel is calculated [20]. In case of metal walls large erosion is appearing, whereas for carbon walls much smaller physical sputtering happens. However, in the case of carbon the release of hydrocarbons is a major problem linked to the sponge-like characteristics of carbon with respect to its interaction with hydrogen. This means that every time the vessel is opened a large amount of hydro-carbons are created due to the interaction of air with carbon. Due to the porous structure of graphite air molecules can diffuse quite deep into the bulk of such graphite tiles and produce there hydro-carbons. The ions bombarding these tiles release this large reservoir and create back-flow to the thruster. Co-deposited layers of hydro-carbons are created in the acceleration channel. These layers are getting conductive hence changing the potentials and produce subsequent problems. The result is a different performance of the thruster in the vessel compared to the one in space. Using instead of carbon metal walls the rates of physical sputtering are larger, but the evaporation at hot channel parts will prevent deposition inside the thruster.
Strategies to overcome this limitations by additional baffles are studied with the help of the Monte-Carlo erosion code. In Fig. 17 the influence of a baffle on the back flow towards the thruster channel is sketched. Particles, accelerated by the thruster, impinge on the vessel walls and sputter its surface. Due to micro roughness of the surface, the distribution function of the sputtered particles follows a cosine function [21]. Its orientation is mainly vertical to the surface, but nearly independent on the incident angel of the accelerated particles. For baffles tilted away from the thruster, this cosine distribution results in a reduced amount of back-scatted particle towards the thruster channel. The calculated angular distribution of the back scattered particles by the Monte-Carlo simulation and from an analytical model are shown in Fig. 16. Due to the larger intersection of the solid angle with the thruster exit plane, the particle flux is larger for small angles of incidence. The returned particle flux consists of two parts, due to the vessel geometry. Sputtered particles from the spherical end are giving the contribution of small angles between 0 and 15 degrees, while for larger angles the sputtered particles are from the cylinder walls. The Monte-Carlo simulation (full line) agrees to the analytical calculations (dashed lines).
In further simulations these results can be used to verify structure of co-deposited layers inside the thruster channel and its impact on the measured thrust. The advantage of the shown simulation is its flexibility to analyze different vessel geometries.
In spite of the low back-flowing fraction of emitted particles (< 1%), these particles are producing codeposited layers during the long test runs causing artifacts in the thrust. A strategy for avoiding this is to implement baffles in the vessel. While the cos(θ) distribution of the sputtered particles does not depend on the angle of incidence of the impinging ion, tilted baffles are turn the direction of the distribution and reduce the back-flow to the thruster exit, as sketched in Fig. 17. A implementation of baffles in the whole vessel, as sketched in Fig. 18 can reduce the back flow by 40% − 66%. The results of this analysis have entered directly the preparation of the life-time qualification set-up for the HEMP-T system foreseen for the SGEO satellite [5]. The vacuum vessel was equipped with baffles tilted by about 20 • to reduce artifacts due to sputtered particles from the vessel walls.

Conclusions
An integrated, modular approach is suggested to address the multi-scale problem of combined thruster-plume models. This approach covers ion-thruster plasma, plume and plasma wall interaction. A hierarchical multiscale set of models in which the parameterization for a lower hierarchy model is deduced from a higher one is proposed. In the frame of such an approach the 3D PIC model could be used to parameterize turbulence effects on the electron mobility in the 2D PIC model. Due to the non-Maxwellian characteristics of plasma the 2D PIC is chosen as the core of the approach for analyzing both in-thruster and plume plasma's. The idea of the "matryoshka-like" set of grids is utilized to reach acceptable length scales with the capability to resolve even the plume kinetically. Furthermore, the modified central different scheme is used to obtain the electric field in order to minimize the error in the momentum conservation introduced by the use of the non-equidistant grid.
For the erosion, redeposition analysis a direct coupling of the kinetic PIC model with binary-collision codes allows a detailed analysis of sputtering inside the thruster. To address longer scales including the plume and its interaction with the satellite Monte-Carlo approaches offer the best perspective, allowing even studies of the interaction of ion thrusters with the walls of the testing facilities and the development of strategies to minimize back-flows from vessel walls to the thruster by using appropriate baffle designs.
For carbon walls the release of hydrocarbons is a major problem linked to the sponge-like characteristics of carbon with respect to its interaction with hydrogen. This means that every time the vessel is opened a large amount of hydro-carbons are created due to the interaction of air with carbon, which can be released due to the interaction with energetic particles.