Global models for radio-frequency ion thrusters

The emerging “new space” age strengthens the importance of rapid development and qualification procedures of electric engines and their peripheral devices. A key element is the reliable simulation of the thrusters and their supply units on short time scales. Global models seem to be well suited for this purpose. In this article, three variants of global models are presented and validated by comparison with experimental results. All models show excellent agreement with experiment, illustrating the strength of this modeling approach. Future developments of radio-frequency ion thrusters can be significantly accelerated with the help of these global models.


Introduction
As of 2021, electric spacecraft propulsion (EP) has become an essential part of modern spaceflight. It enables highly mass-efficient and, hence, economically beneficial missions with prolonged durability in contrast to conventional chemical propulsion. Two prominent applications come to mind for which EP had caused a disrupting influence: geostationary/geosynchronous station keeping (GSK) and electric orbit raising (EOR) of telecommunication satellites [1]. Typically, small chemical orbital propulsion devices such as cold gas or mono-propellant thrusters had been used to compensate for ecliptic effects and solar winds during GSK missions, which led to a high wet-to-dry-mass ratio of the satellites due to low specific impulse I sp of said propulsion devices. Additionally, to bring the geostationary telecommunication satellites (geo-sats) into their orbits of operation, the launcher would take the satellites directly to a geostationary elliptical transfer orbit (GTO). Typically, chemical mono-or bi-propellant thrusters had been used to initiate a kick-burn at the apogee to ultimately bring the telecommunication satellites to their final orbits. Those kick-burns consumed a large amount of propellant due to comparably low I sp of the so-called "apogee engines". The propellant required contributed significantly to the launch mass and volume of the satellite. However, the complete transfer only took about one week until the satellites would be fully operational in orbit. Within the last recent years, substantial financial savings could be achieved using EP for GSK and EOR [2]. Due to the much higher I sp (up to a factor of 10) in comparison to typical chemical devices, considerable economic advantages occur because the amount of propellant to be launched could significantly be reduced. Since the total cost of a mission relates strongly to the absolute mass to be launched, either considerable financial savings or a substantial increase of dry mass might be achieved. Nowadays, most geosats are brought to low earth orbit (LEO) with the launcher. With the use of continuously operating electric thrusters, a constant (electric) orbit raising (EOR) is achieved, which consumes only fractions of the propellant mass typically associated with an orbit raise. However, compared with chemical engines, EOR requires a significantly longer period of time for spiral-up due to the low thrust of EP devices (up to 6 months). Given the advantages listed above (most prominently the fact that the geo-sats can use the same thrusters for EOR and GSK, ridding the need of different systems with different piping, valves, tanks, etc.), many satellite operators favor this option nowadays.
There are different types of EP devices available and ready for those applications [3,4]. Hall effect thrusters (HET) have been flight-proven since 1971 [5] and have, up to now, been the first choice for EP demands. However, radio-frequency (RF) ion thrusters represent a viable alternative. This is mainly due to their typically almost twice as high I sp , which resides in the area 3, 000 − 4, 000 s. Additionally, plasma generation and extraction can be regarded as decoupled mechanisms which facilitates scalability of RF ion thrusters. This makes them, in principle, the optimal technology for EOR and GSK. Unfortunately, RF ion thrusters often struggle with lifetime issues which result from extraction grid erosion. In contrast to HET where the plasma ions are extracted through an annular orifice (for a detailed description of HET refer to Ref. [3] and references therein), in RF ion thrusters ions are extracted from the plasma using a gridded extraction system; hence their classification as gridded ion engine (GIE). An inherent issue of GIEs is grid erosion effects caused by high-energy ions bombarding the grids directly (direct impingement) or charge-exchange collisions (CEX) between fast ions and slow neutrals which exchange momentum of those species. The latter results in fast neutrals and slow (positive) ions which are attracted by the negative potential of the accelerator grid, which in turn manifests itself in so-called "barrel erosion" of the accelerator grid [6].
To overcome those drawbacks, numerical analyses are often performed starting at a very early stage in the thruster development process. Since timescales are limited in industry and detailed knowledge of the internal microscopic processes is often not needed, global and hybrid models which are trimmed to speed are the preferred choice. Therefore, a precise knowledge of details such as the electron energy distribution function, which can be determined by PIC-MCC approaches [7], is not necessarily required, as long as the accuracy of the determination of quantities such as the required input power of the RFG is not diminished significantly. Furthermore, global models can be rather easily used to simulate plasmas fed by molecular propellants, such as an iodine [8] or oxygen [9]. For both iodine and oxygen, the discharge can be electronegative, which shows that global models can even be applied to such plasmas. This a remarkable advantage of modeling ion sources used for materials processing using global models, since those often use molecular gases, which exhibit a much more complicated plasma chemistry in comparison with noble gases. Electronegative gases, such as iodine or sulfur hexafluoride are also investigated as propellant for a neutralizer free electric propulsion system, where both negatively and positively charged ions are extracted in order to avoid a charging of the spacecraft. This principle has been introduced in the frame of the project PEGASES (plasma propulsion with electronegative gases) [10].
In this paper, we want to analyze different approaches of global (and hybrid) RF ion thruster models published in literature, discuss their pros and cons as engineering tools during thruster development and try to give recommendations how to integrate global and hybrid models into the development steps of a real thruster. For this purpose, we briefly introduce the most relevant physics of RF ion thrusters in "Physics of RF ion thrusters" section, then introduce three different, but nevertheless somewhat similar models in "Global RF ion thruster modeling" section: pure global (furthermore referred to as 0D model), 0D/2D (referred to as 2D model) and 0D/3D (referred to as 3D model) hybrid models. The nomenclature is due to the specialized handling of the plasma and electromagnetic phenomena, respectively. The basic plasma processes are described in a global (0D) fashion, albeit with slight modifications depending on the attached electromagnetic model (0D, 2D, or 3D). We work out the theories used for each of the models in detail in the respective sections. A benchmark thruster (RIT-4), as depicted in Fig. 1, is analyzed with each model in "Comparison and discussion of the results" section and their respective outputs are compared and validated by experiments in "Experimental validation" section.

Physics of RF ion thrusters
The foundation of RF ion thrusters, which were formally invented by H. Löb in the 1960s [11], are inductively-coupled plasma (ICP) sources which have been and still are widely used for materials processing purposes [12][13][14]. The advantageous characteristics of ICP sources for materials processing purposes, e.g. scaleable ion beam energy and low divergence angle, are also what makes them interesting for space propulsion. Low divergence angles are correlated with high specific impulse which relates the effective axial exhaust velocity of the thrusters c eff to the gravitational acceleration on the Earth surface g 0 : Fig. 1 Left: RF ion thruster assembly including peripheral HV supplies (negative: NHV, positive: PHV), mass flow controller (MFC) and radio-frequency generator (RFG). The screen grid is floating, since the PHV-supply is connected to the gas inlet instead of the screen grid. Right: CAD drawing of the RIT-4 benchmark thruster used in this paper incl. coil current, electric and magnetic field High I sp implies high mass utilization efficiency which can be deduced from the equation for thrust: This relation is derived directly from Newton's laws (mainly Netwon's third law " actio=reactio" reflecting momentum conservation) and is an essential figure of merit for space propulsion systems. Substituting Eq. (1) in Eq. (2) results in which shows that for constant thrust F low mass consumptionṁ p corresponds to high I sp . From Eq. (1) it is evident that high I sp is based on high effective axial exhaust velocity which can be obtained by adjusting the extraction grid voltages to the ICP parameters accordingly. The grid voltages V screen and V accel correspond to the screen grid and the acceleration grid potentials, respectively. V screen determines the kinetic energy and hence the exhaust velocity of the extracted ions of charge-to-mass ratio q i /m i reads: This is an approximate relationship due to the negligible influences of the ion flux towards the walls induced by the voltage drop V p V screen between the plasma center and the walls as well as, secondly, neutral particles escaping the thruster with only thermal velocity. Furthermore, η denotes the divergence angle due to radial velocity components which will not contribute to thrust, if a perfect axi-symmetric ion beam is assumed.
The extraction of high-velocity particles will only be accomplished, if the plasma yield (energy and density of ions at the plasma sheath-edge) is high enough. In other words, there has to be a sufficient amount of flux of charge-carriers (electrons and ions) towards the extraction grid to meet the amount of thrust to be generated. The thrust can also be expressed as a function of the ion beam current I b and the effective exhaust velocity, as can be derived when substituting Eq. (4) in Eq. (2): Equal amounts of flux of electrons and ions characterize bounded plasmas. Since in all plasmas quasi-neutrality is demanded; i.e., positive charge must equal negative charge within the bulk volume on a macroscopic scale, and since many technological plasmas, such as those used in RF ion thrusters, consist of mainly electrons and singly positively charged ions, the number density of both species is approximately equal: This equation, however, does not hold true in the plasma sheath between the bulk of the plasma and structural walls of the discharge vessel. Due to their significantly lower mass, electrons will traverse the sheath layer much faster than ions, leading to negative charging of the walls with respect to the plasma bulk. This negative potential promotes further ion flux and prohibits further electron flux towards the walls and ultimately a steady state with equal fluxes of ions and electrons arises, resulting in stable plasma parameters with n + n − in the sheath layer [12]. The particle flux is a direct consequence of the plasma parameters and thus related to the input parameters coil current I coil , radio-frequency f, and propellant mass flowṁ p . The (inter-)dependence of plasma parameters, electromagnetic fields and input RF coil current can be described using a time-harmonic diffusion equation of the electrodynamic vector potential In this equation, J s denotes the source current density generating a vector potential A and ω = 2πf . Throughout this paper, the underline denotes a complex quantity.
In typical low-temperature plasmas, the permeability approximates to the vacuum permeability (μ ≈ μ 0 ) and is thus independent of space. The interconnection of plasma and EM (electromagnetic) fields is expressed by the conductivity κ, which can be correlated to plasma density and electron temperature in terms of collision parameters. For this purpose, Drude models are often used for global (RF) plasma modeling [13,15,16]: In this equation, ω pe = n e e 2 / (m e 0 ) 0.5 denotes the electron plasma frequency; i.e., the collective oscillation of the electron "cloud", 0 the vacuum permittivity and ν eff the effective collision frequency of electrons. Strictly for RF plasmas, the ions are often regarded as immobile with respect to the quickly alternating EM fields due to their inertia and thus are not taken into account as primary particles in Eq. (8). They are, however, considered as background particles with whom the electrons may collide, thus they contribute to the effective collision frequency: The first term in this equation is the elastic electron-neutral collision frequency which is derived from the collision cross section and the second term describes the above mentioned Coulomb collisions between electrons and ions. Detailed information about these processes and the underlying cross sections is given elsewhere [12][13][14][17][18][19][20]. Equation (8) determines the electric power needed to sustain the plasma, depending on the stationary conditions of the plasma discharge parameters, mainly electron temperature T e , neutral gas density n 0 and bulk ion density n + . Equation (7) can be solved for A and related to the electric field vector by assuming no scalar potential gradient within the plasma bulk area [21], leading to The absorbed RF power is then given by volume integration of the power density, with E = E . To correlate the RF power to the propulsive performance of a RF ion thruster, the generated beam current I b has to be further examined. The beam current is given by with the elementary charge e, the ion density at the plasma pre-sheath edge n +,s , the Bohm velocity of ions u B = (k B T e /m i ) 0.5 , the total cross section area of the extraction grid A g and its ion transmission factor β i [12,22]. Only singly positively charged ions are assumed here which holds true to a certain extent for most RF ion thruster operational regimes [22,23]. As can be deduced from Eqs. (11) and (12), the ion density as well as the electron temperature, in terms of u B , relate the ion beam current to the electromagnetic field effects inside the plasma chamber and thus the necessary RF power to sustain the corresponding state of the plasma.
An important quantity is the ion transmission factor β i which is self-consistently dependent on the plasma parameters (T e and n − ≈ n + as well as n 0 ). It is given by the ion beam current related to the current traversing the plasma (pre-)sheath. Thus, it can be obtained by analyzing the ion trajectories through the extraction system. If a strong defocusing of ions occurs, β i will be reduced, also reducing the beam current. Hence the thrust and overall propulsive performance of the thruster will suffer. Additionally, the beam divergence is affected, which is hence inherently coupled to the transmission factor. More information on ion optics in general can be found in [24][25][26][27][28][29][30] and in particular, in terms of transmission factor modeling, will be shown later in the paper.
To focus the ions and accelerate them to high exhaust velocities, screen voltages in the kV range are typically used. The power dissipated by the extraction grid system is given by where V screen is the screen grid voltage. The acceleration grid is negatively biased with V accel . Purposes of this approach are the above-mentioned focusing of ions and the repelling of electrons to avoid electron back-streaming (EBS) which changes operational conditions for the thruster and hence gives rise to adverse eroding mechanisms. However, since typically I accel I b and |V accel | < V screen , there is hardly any power deposited on the acceleration grid. The same applies for the ground grid in 3-grid-systems.
Altogether, the power to sustain a given mode of operation is given by the sum of Eqs. (11) and (13). The powering electronics can be divided into RF and DC components. In the remainder of the paper, only the RF electronics will be analyzed in greater detail, since it comprises an inextricable part of the plasma generation system. The DC part can easily be decoupled as DC electronics can be considered linear under most operation conditions.
To ensure efficient energy coupling from RF generator (RFG) to the thruster plasma, its output frequency is often matched to the resonant frequency of the thruster impedance, which is determined by the ohmic-inductive plasma load and a series resonant capacitor. Additionally, the influence of the feeding cable has to be accounted for as well. More information on the subject is given elsewhere [31,32].
The thruster impedance can be correlated to its propulsive performance and hence its input terminal's parameters (voltage, current, frequency) following Poynting's theorem: Radiative processes are neglected, which holds true due to comparably low frequency (i.e. large wavelengths of about 100 m) and short conducting structures which do not form any kind of waveguide. The impedance Z = R + iωL can be derived from Eq. (14) following the procedure given in Ref. [21].

Global RF ion thruster modeling
Global plasma and ion source models have been reported on since the early days of exceedingly increasing computing capabilities. Most of them can be traced back to Ref. [33] and references therein. Inductively-coupled plasma sources have been used most prominently for materials processing applications. Hence, a majority of existing global models focus on plasma sources of the inductive type [15,16]. They are regarded as basis for many published global (0D) RF ion thruster models [8,22,[34][35][36][37][38]. Those models have in common that they treat all of the underlying physical mechanisms as described in Eqs.
(2)- (14) in a global, volume-averaged way. This speeds up simulations and allows one to derive, even for large thrusters, results in a reasonable time frame. Thus, it provides a useful virtual prototyping tool for thruster development. There are, however, certain drawbacks and restrictions of global models. Especially, if strong geometrical features break symmetries. This is typically the case for RF ion thruster induction coils which are often electrically short. The magnetic field generated inside such coils cannot be considered purely axial; especially at the coil ends there will be a strong radial component which builds up approaching these limits. Those radial components will have a measurable effect on the coupling of RF energy into the plasma and may lead to significant errors if not properly taken care of (scaling parameters, fit functions, etc.).
Most of the problems induced by global assumptions have been solved by 2D axisymmetric approaches. In some models only the EM solver is capable to handle 2D geometries and the plasma is still regarded as homogeneous [34]. There are more consistent approaches which pre-define a density distribution [39][40][41], either derived from empirical studies [42], fluid models [43,44] or full particle-in-cell/Monte-Carlo-collision (PIC/MCC) approaches [7,45,46]. It is evident that fluid and especially PIC/MCC models aim at a totally different output (deeper understanding of microscopic plasma behavior, evolution of statistical distribution functions etc.) than global models, which, due to typically long simulation duration, makes them inappropriate for virtual prototyping applications.
To enable the most realistic representation of the EM fields inside actual thruster geometries, 3D models have to be used. To keep simulation duration at a minimum, even those models are coupled to homogeneous or pre-defined plasma density profiles [47,48]. This approach is considered as the most reasonable trade-off between accuracy and simulation duration. However, as will be shown in "Comparison and discussion of the results" section, for typical cylindrical discharge vessel and coil geometries, the 2D approach should be a good enough representation of the actual physics involved-even though simplified solvers are used. The models developed at our respective institutes are shown in the following sections. Table 1 lists the main working mechanisms and differences of the proposed models.

0D model
The global model examined here is based on a publication by Chabert et al. [36]. In addition, it has been extended by linking it to the software IGUN [49] to adequately model the ion optics of the thruster's grid system.  Individual treatment of each beamlet as a function of n + (r) and τ i (r) Individual treatment of each beamlet as a function of n + (r) and τ i (r) Peripheral electronics

--LTSpice simulation
The global model is composed of four coupled differential equations describing the time evolution of neutral gas density n 0 , ion or electron density n + , neutral gas temperature T 0 , and electron temperature T e . All rate coefficients are calculated using the cross sections provided by the "Biagi database" [50] assuming a Maxwellian electron energy distribution function. While the 3D model uses the same cross sections, the 2D axi-symmetric model uses those provided by the "Hayashi database" [51]. Assuming a Maxwellian distribution function is quite common in global modeling, due to simplicity and a lack of precise knowledge of the actual distribution function in the plasma. The actual distribution may, among other factors, even depend on the thruster's operational point. Multiply charged ions are not considered and due to quasi-neutrality of the main plasma, ion and electron density are described by the same equation in the model. The ion and neutral gas temperatures are assumed to be equal. The input parameters used to represent the geometry are length l DC and radius of the discharge chamber R in addition to the number of coil windings N Coil . The coil windings are considered to be evenly distributed over the entire length of the cylindrical discharge chamber. The transmission properties of the gridsystem are described by two dimensionless parameters for the neutral gas β n and the ions β i . In a first approach, the efficient transparency for neutral gas β n can be calculated as the product of the geometrical transparency of the screen grid, with respect to the cross section corresponding to the diameter of the discharge vessel, and the total transmission probability P Total through one extraction channel. P Total is determined as a serial connection of pipes of arbitrary length [52]: • A 1 : Cross section of a aperture in the screen grid.
• A i : Cross section of element with the index i. • P Total : Total transmission probability.
• P i : Transmission probability of element of index i.
The transmission probability for one element of the extraction channel P i , either an aperture, or a part of the inter space between the grids, is considered and calculated the with the empirical formula for cylindrical pipes of arbitrary length l and diameter d [52]: This approach typically estimates slightly smaller values, about 5 % to 10 %, of β n compared to the more advanced methods, such as described in "2D axi-symmetric model" section. For this reason, this approach was not used for the simulations of this publication. Instead, the value obtained by the more advanced method presented in "2D axi-symmetric model" section was used; yielding a better comparability of the results. Furthermore, the frequency of the RF current in the coil, the ohmic resistance of the coil at this frequency, the propellant mass flowṁ p injected into the discharge chamber and the RF power P RF are also input parameters of the model. From the electron temperature, ion density and neutral gas density, the ion current generated by the thruster can be calculated. In contrast to the works of Chabert et al. the neutral gas temperature was kept constant in order to be able to compare this modeling approach to the other two models.
Based on the electron temperature and ion density as well as the grid geometry and applied voltages, IGUN allows the extracted ion current to be calculated [49]. The ion density is weighted by the factor h l , which describes the edge-to-center plasma density ratio on the grid surface, before entered into IGUN. For calculating the factor h l as well as h r the solution for the intermediate pressure regime introduced by Godyak [53] is used. The ion current obtained that way serves as a comparison for the ion current obtained from the global model in order to be able to find a value for the ion optics parameter β i . This must be determined iteratively by the aforementioned comparison, since the ion optics parameter both depends on and affects the particle densities and electron temperature.
To determine the plasma parameters and the necessary RF power at a given mass flow and beam current set point I b,set , the procedure shown in the flow chart in Fig. 2 is applied.
With the given thruster parameters, the differential equations from Chabert et al., denoted as ODEs in Fig. 2, are solved. The values of the densities and temperatures in steady state are used to calculate the beam current I b . The RF power is adjusted until the beam current matches the beam current set point. Then the output values of the differential equations are used to perform a calculation with IGUN. The beam current output from IGUN (I IGUN ) is compared to the beam current from the differential equations. If the currents do not match, the ion optics parameter in the global model will be adjusted accordingly. If these currents are identical, it must be verified if the beam current still matches the beam current set point. If this is not the case, the RF power will be adjusted accordingly and the process will be repeated. If the beam current of the global model and the beam current set point match, the calculation at the specific operating point will

2D axi-symmetric model
A quite detailed description of the 2D model is given in Ref. [41]. Here, only the most important features should be addressed to give the reader an idea about the model's inner mechanisms and how to work out the main differences to the other presented models. It is noteworthy that within the 2D axi-symmetric simulation domain, which is described by cylindrical coordinates (r, φ, z), all quantities x are assumed to be constant with respect to the azimuth angle φ, leading to simplified equations. Hence, all quantities x that can be deduced from the vector-potential are independent of φ and, therefore, ∂x/∂φ = 0. This leads to simplified equations. A flow chart of the 2D model is given in Fig. 3.
Here, the grid transmission coefficient for ions β i is composed of several single transmission coefficients τ i ; one for each beamlet, depending on the plasma density adjacent to each simulated grid cell. Their values are stored as functions of T e and n + in a database and their values are interpolated when accessed by the global model. The ion transmissions coefficients are obtained by an in-house 2D trajectory tracker based on the PIC method [54,55], as discussed in "Ions" section. The transmission coefficient for neutrals β n is pre-determined using Sandia National Lab's SPARTA DSMC package [56,57]. It is an explicit function of the extraction grid geometry since no inter-molecular collisions are assumed. Further information can be found further below in this section. The set point parameters include the angular frequency ω, the beam current I b,set , the extraction voltages V screen and V accel , the propellant mass flowṁ p and the temperature of the propellant gas and all of the thruster components T 0 . A full thermal equilibrium of neutrals and ions is assumed. Additionally, initially best guesses of the bulk ion/electron density n + /n − and the Maxwellian electron temperature T e are used to start the iterative operation.
Before the iterations start, the neutral gas density, which is assumed constant throughout the volume, can be derived as a function of the set point beam current: Here, v th = (8k B T 0 / (πm)) 0.5 denotes the thermal Maxwellian velocity magnitude of the neutrals, with their mass m and A g the extraction grid area.

Charge conservation
The innermost loop is a charge conservation equation which relates the volume ionization rate to the charge lost due to flux out of the plasma: In this equation, K iz = K iz (n 0 , T e ) denotes the rate coefficient for ionization of neutrals by incident electrons and is strongly coupled to the underlying distribution function of the particles, which is here assumed to be Maxwellian (electrons, ions and neutrals are all assumed to obey a Maxwell distribution in phase space; however, since T e T i ≈ T 0 , only ions (T i ) and neutrals (T 0 ) are assumed to be in thermal equilibrium with the structural components). R furthermore denotes the upper radial, z min and z max the lower and upper axial boundaries of the plasma in the simulation domain, respectively. H is a heuristic density distribution function first introduced in Ref. [39]. It is based on analytic sheath-to-bulk density factors [33,42]. It follows a parabolic shape and can thus be used to examine most of the operational points of an RF ion thruster plasma [58]: Here, h r and h l denote the sheath-to-bulk density ratios which depend strongly on n 0 and z 0 denotes the axial center of the plasma. Equation (19) is solved by means of numerical integration and Newton's method for T e .

EM field solver / energy conservation
With T e and n 0 , the plasma conductivity from Eq. (8) can be formulated. The structural components of the thruster, such as gas inlet, induction coil, extraction grid and housing, are parameterized by temperature dependent conductivity models following first order approximations, with the electrical resistivity ρ 0 at temperature ϑ and an empirical fit parameter α which is provided for many metals by databases like NIST. A thorough overview is given in Ref. [59]. The inductively coupled RF power from Eq. (11) is obtained by a finite-difference representation of Eq. (7) which is numerically solved [47]. This absorbed power within the plasma region is computed against the power lost to collisions and wall flux and solved for n + : The collision and wall flux losses are given by (with the temperatures T e and T 0 given in V ): where K m , K ei , K iz , K ex , E iz and E ex denote the rate coefficients for momentum transfer (electron-neutral), Coulomb collisions (electron-ion), single ionization (electron-neutral), excitation (electron-neutral), the first ionization energy threshold and the excitation energy threshold. For the latter, about 50 states are incorporated in case of xenon as propellant gas. The voltage drop across the plasma pre-sheath and sheath region V p is furthermore given by with the mean Maxwellian electron velocity v e = (8eT e / (πm e )) 0.5 . Equation (22) can finally be solved for n + in an outer iterative loop.

Extraction and performance model
With T e and n + converged to stable values, a variation of Eq. (12) can be solved for the beam current. The corresponding ion transmission parameter β i (T e , n + ) is calculated from each τ i stored in a look-up table which is accessed during evaluation from several sub models. Therefore, it is advisable to obtain this table prior to the global model simulations which will be shown in the next section.
If the calculated beam current does not match the desired set point value (I b = I b,set ), the coil current is updated by a proportional-integral controller. After convergence, the thruster performance (Eqs. (3) and (5)) and impedance (derived from Eq. (14)) are calculated.

Transmission factor modeling
This section gives a brief introduction of the two sub models that determine the transmission factors for neutral flux β n and ion flux β i . Both models need the actual extraction grid geometry as input.

Ions
The ion transmission coefficient β i is determined by a trajectory tracking code loosely based on the PIC method. In this case, axi-symmetry is assumed, hence only the center beamlet can be evaluated. If the distance between adjacent beamlets is large enough, the error is negligibly small as can be demonstrated by comparison to full 3D PIC solutions [60,61]. A flow chart of the trajectory tracking code is shown in Fig. 4. The scheme follows approaches shown in Refs. [49,60,[62][63][64] with a special treatment of the differential operator for cylindrical axi-symmetry. In principle, a non-linear static electric potential equation is solved iteratively while the ions traverse the simulation domain until all trajectories have exited the modeled domain of the extraction channel. The numerical solution is rather complicated due to the non-linear term in the Poisson equation, which is given by the Boltzmann factor for the electron density at the sheath edge, with 0 denoting the potential of ions and electrons when entering the simulation domain at the plasma sheath-edge. It arises due to the negative wall potential with respect to the Maxwellian plasma bulk and hinders electrons to traverse the sheath unless they have a sufficiently high kinetic energy. It is however crucial to incorporate this term due to otherwise non-realistic ion trajectories, since all particles are extracted from a quasi-neutral plasma and the sheath is governed by space-charge effects.
After each successive potential update, a chain of operations is performed which a) calculates the space-charge of ions within the simulation domain ("Scatter"), b) maps the electric field vectors, which are stored on discrete grid points, onto the charges which can be located at arbitrary positions ("Gather"), c) uses those electric field values as driving forces in an equation of motion q p E = m p dv p /dt yielding updated velocities and d) pushes the ions to new coordinates with the updated velocity values ("Push"). For the mechanisms in c) and d) commonly the Leapfrog method is used where velocity and spatial vectors are staggered by a half time step size to increase numerical stability [65]. After all ions have either collided with structural parts or have left the simulation domain, the resulting space-charge map is checked for changes compared to its preceding iteration. If this is the case, a new ion density map is generated (as a function of said space-charge map) and fed back to the Poisson equation. The loop is repeated until all quantities converge within a given tolerance.
The grid ion transmission factor β i is then calculated from the transmission factors for each beamlet τ i . Those can be broken down to functions of the current entering the simulation domain I in and the beam current I b leaving the domain (the sum of all beamlet currents I bl ):

Neutrals
The transmission coefficient for neutrals is an explicit function of the grid geometry.
Here, a 2D representation induces non-negligible errors which will ultimately lead to very wrong power balances. Hence, a 3D solution for this particular problem is chosen. The problem is solved using Sandia's SPARTA (stochastic parallel rarefied-gas time-accurate analyzer). The simulated section is shown in Fig. 5 and manifests as a unit-cell with symmetric boundary conditions to all lateral boundaries (not in flow direction). As can be seen, with this approach, only a small fraction of the total grid system has to be simulated. Due to the symmetry boundary conditions a hexagonal aperture pattern arises which is typical for RF ion thruster grids.
The flow of particles is considered to be free-molecular, so no particle collisions are implemented. Instead, only diffuse reflections with structural components are taken into account.
The transmission coefficient is obtained in a similar manner than the one for ions, simply by relating the flux of particles at the downstream boundary γ ds to the fixed number at the inlet γ in and scaling to the full grid geometry: In this equation, N denotes the number of extraction beamlets and A SPARTA the rectangular simulation domain area in SPARTA representing half an extraction beamlet, as shown in Fig. 5.

3D model
The model described in this section is an extension of our previously published model [48]. An in-house 3D ion extraction code, published in Ref. [61], is used to simulate a unit cell of the extraction system with the geometry shown Fig. 5. The simulated region is chosen to represent one unit cell of an indefinite large array of extraction channels arranged in a hexagonal manner. This approach with an 3D implementation takes into account the effects of adjacent holes and beamlets on the ion optics behavior; in contrast to the axi-symmetric treatment presented in the last section. The same 3D approach is used to simulate the neutral gas transmission coefficient by utilizing the test particle Monte-Carlo (TPMC) method [66]. This method is based on simulating the trajectories of multiple test particles. The transmission coefficient is given by the amount of test particles leaving the ion thruster divided by the amount of generated test particles. The particles are generated at a boundary representing the plasma. Hereby it is important to use statistically evenly distributed positions and velocity vectors. An isotropic velocity of neutral particles inside the plasma is assumed. The trajectories are calculated by integrating the equation of motion and applying a mirror reflection at the boundaries towards neighboring simulation domains and an explicitly diffusive reflection when colliding with an extraction grid. The neutral gas density as well as the temperatures of neutrals, ions and electrons are assumed to he spatially homogeneous. As well as in the 2D model a density profile inside the plasma for the electrons and ions is assumed, which is considered consistently in each sub model. For the calculation of this profile h parameters describing the sheath-to-bulk density ratios are used, as described in Ref. [48]. The h parameters are based on analytic equations valid for different density regimes. We use a fitting equation between different regimes introduced in Ref. [53]. The profile already described by Eq. (20) was separated into two parts in order to better cope with non-cylindrical geometries. For this, the center of mass of the ionization chamber is calculated and assumed to be the center of the plasma. Starting from this point the distance R towards the most distant point in radial direction, the distance l p towards the most distant point in positive axial direction and the distance l n towards the most distant point in negative axial direction are calculated. Based on these distances the corresponding h parameters h r , h l p , and h l n are calculated. The profile is then calculated with where l is the coordinate in axial direction with l = 0 m at the position of the plasma center and r is the radial distance from the plasma center. Furthermore, the plasma conductivity model was updated. In the model used the electromagnetic fields are solved in 3D. The calculation of losses inside the plasma is done similarly to Eqs. (23)- (27) described for the 2D model. An integration method was implemented to enable the simulation of arbitrary geometries. Therefore, the 3D model also includes the gas inlet. While P iz , P ex and P w are calculated the same way as in the 2D model, we neglect P ei and use a different equation to calculate P m . The latter difference is not described in detail because the term is of negligible influence. The equations used to calculate the amount of generated and lost ions are also the same.

Solving the self-consistent set of balance equations
To reach a self consistent-solution we use the approach presented in Ref. [67] and updated it to also incorporate our new peripheral model, which is described in the following subsection. The corresponding flow chart is illustrated in Fig. 6. By using multiple conditions, the balance equations are decoupled to be able to calculate the three unknown state quantities (neutral gas density n 0 , electron temperature T e and ion density n + at the plasma center) sequentially. First, we define the set point values. In particular, we set the neutral gas temperature T 0 and the material temperatures T m to 423 K. Furthermore, we set the radio-frequency f, the grid voltages V ext and the geometry. The neutral gas density is then calculated in the same manner as already mentioned in "2D axi-symmetric model" section. In the 3D model, we iterate n 0 until the balance equation of neutral particles (denoted as mass conservation in Fig. 6 is satisfied. After n 0 has been obtained, we assume that the density profile of electrons and ions (Eq. (32)) does not depend on the value of the ion density in the plasma center n + . Accordingly, the solution of the charge conservation; i.e., the balance equation of ions generated inside the plasma and ions lost at the edge of the plasma, is independent of n + . Therefore, any value of n + can be chosen to solve the equation for T e . This is done by iterating the electron temperature T e until a value is found which satisfies the balance equation ( I ≈ 0 A).
After n 0 and T e have been obtained, the spatial ion extraction can be calculated with a given value of n + . This is done by iterating the value of n + until the simulated beam current I b matches the beam current set as boundary condition. After this, all state quantities (n 0 , T e and n + ) have been obtained. With these values the plasma losses P n as well as the spatial plasma conductivity κ p are calculated. The latter is used for the electromagnetic field solver to simulate the impedance of the thruster system Z th which includes the plasma resistance R p .
With R p and P n known, the coil current I coil is calculated with I coil = P n /R p 0.5 . Finally, the simulated impedances and I coil are used inside the peripheral model, in which die DC input voltage V DC of the RFG is iterated until the simulated coil current matches I coil .
Then the input power P DC , the input voltage V DC and the input current I DC are known.

Peripheral model
The utilized radio-frequency generator (RFG) was developed at the Technische Hochschule Mittelhessen -University of Applied Sciences (THM). It is based on a field programmable gate array (FPGA) processing a digital phase and frequency control algorithm presented in [31,68]. The algorithm is used to generate the driving signal for a half-bridge consisting of two GaN MOSFETs of type "GS66516T-E02-MR" from the manufacturer GaN Systems, with the hardware assembly shown in [32]. The rectangular output signal is used to excite a resonant circuit consisting of the coil inductance and a resonant capacitance with its resonant frequency. To be able to simulate the switching losses for this setup, a more sophisticated approach than the one presented in Ref. [48] was necessary. Instead of an analytic model a simulation with LTSpice is used. The standard approach would be to use a SPICE model from the manufacturer to describe the MOSFETs in LTSpice. Utilizing these models only the currents and voltages at the MOSFETs connection interfaces; i.e., drain (D), gate (G), and source (C) are visible in the simulation. Given this information, it is difficult to deduce the internal behavior and the resulting losses. For example, turning on the MOSFET while a voltage is applied between source and drain leads to the discharge of the parasitic drain-source capacitance within the MOSFET. This generates losses inside the MOSFET, but the responsible currents do not appear as external currents measurable at the terminals. In addition, the drain-source capacitance of the "GS66516T-E02-MR" is up to 1 nF, which may lead to significant displacement currents. Depending on the topology, these are, due to the non-linear behavior, hard to distinguish from the active current. To allow one to evaluate the switching process in detail as well as the occurring switching losses one has to simulate the internal mechanisms of the MOSFETs used. A schematic image of the MOSFET with this modeling approach is given in Fig. 7. The values of the capacitance C GS (between gate and source), C GD (between gate and drain) and C DS (between drain and source) are derived from data sheet values for C ISS (small signal input capacitance), C OSS (small signal output capacitance) and C RSS (small signal reverse transfer capacitance). While C GS is nearly linear (we use C GS = 510 pF) the values of C GD and C DS are set up as look up tables to describe their non-linear behavior. The value R G = 1.5 for the internal parasitic gate resistance is also given in the data sheet. All variables with a dash correspond to internal variables of the MOSFET.
The MOSFET's conductance is modeled by a current source for I DS based on fitting the data sheet parameters. The current I DS is calculated by with the threshold voltage V th , the drain-source voltage V DS and the internal gate-source voltage V GS . First we define based on the data sheet. V p and I p define the pinch-off point between linear and saturation region. For V GS > V th (linear region), the term I p 1 − V DS − V p 2 /V 2 p is used, which describes a parabolic behavior of I DS with respect to V DS . At V DS = 0 V we have I DS = 0 A and the vertex is at I DS = I p with V DS = V p . For V DS ≥ V p (saturation region), the current I DS is assumed to be constant over V DS . Since I DS = f (V DS ) both I p and V p change with V GS . This is taken into account by describing I p and V p as functions of V GS . For I p , the fit is used to describe the vertex of the parabola (coinciding with the current at the pinch-off point). For V p , the fit is used, so that the initial slope of the chosen parabola describes the switching resistance between drain and source R DS(on) which is approximately 40 m . The fitting is performed for an assumed junction operating temperature of 343 K with the values given in the data sheet (Rev. 161007).
The reverse conductivity is modeled by a diode. The complete peripheral model is shown in Fig. 8.
The two-wire cable is described by inductance L cab = 140 nH and resistance R cab = 43 m . The 0.15 m long transmission line of type "Belden 9220" is described by a -model with the transmission line's input and output capacitance C tl = 7.6 pF, its inductance L tl = 38 nH and its resistance R tl = 7 m . Furthermore the resonant capacitance inside the RFG C res = 2.725 nF is set and its equivalent series resistance R esr is calculated by its loss tangent tan(δ) = 0.001. The voltage sources for the gates are set to the given excitation frequency with the corresponding delay times.

Experimental validation
For validating the presented global models, a comparison with experimental results is mandatory. A meaningful comparison can only be drawn if several requirements are fulfilled. Most importantly, the pressure in the vacuum tank has to be low enough to reduce backflow of neutrals into the thruster to a minimum. Otherwise this effect is a source of error when comparing the results of the simulations and the experiment.

Experimental setup
The cylindrical vacuum chamber BigMac Evo of about 3 m in length and a diameter of 1.6 m was used for the experiments. Operating the thruster, xenon is mainly pumped by two cryogenic pumps with a total pumping speed of 34, 000 l s −1 . Therefore, the xenon partial pressure does not exceed 10 −4 Pa for flow rates up to 100 μg s −1 . Assuming a homogeneous pressure distribution inside the tank and taking the conductance of the thruster's grid of approximately 5 l s −1 into account, the backflow of xenon into the thruster is only about 0.01 μg s −1 and, therefore, negligible. The actual backflow might be slightly higher, since the pressure along the thruster's beam axis can be assumed to be slightly larger than the mean value. However, even if the actual backflow would be an order of magnitude larger, it would still be very small compared to the flow injected into the discharge vessel. A turbomolecular pump of 2, 350 l s −1 nominal pumping speed for nitrogen is used for the pumpdown and later on keeps the residual pressure of nitrogen and oxygen on an acceptable level; since the cryo pumps can not continuously be kept cold enough for pumping those gases. The pressure in the vacuum facility is monitored by a Pirani/cold-cathode combination gauge from by Pfeiffer Vacuum. The thruster is mounted radially in the center of the vacuum tank and is attached to an ISO-K-160 flange. A test power supply (TPS) unit is used for operating the thruster. The TPS contains the power supplies for controlling the RFG, the screen and acceleration grid of the thruster. For providing power to the RFG, a "TDK-Lambda GENH-60-12.5" is used. The screen and accel grids are supplied by "FuG MCP 700-2000" and "FuG MCP 35-2000" sources, respectively. The deceleration grid is directly connected to the facility ground. Using an ordinary computer, all power supplies as well as the mass flow controller, type Bronkhorst "El-Flow Select" with a maximum range of 140 μg s −1 xenon, are digitally controlled by a LabView-based control software. A beam current controller is implemented to keep the current flowing to the screen grid constant. This current is adjusted by tuning the RFG's input voltage. As requirement of stability for each data point of the performance-mapping, the largest tolerated drift of the RFG's input power is 1 % in a time interval of 30 min. In order to fulfill this condition, the overall time span operating the thruster was 6 h.
For further validation of the global models, the current through the coil is monitored with a Pearson probe 2878.

Measurement results
The recorded performance mapping is shown in Fig. 9 a). The input power of the RFG decreases strictly monotonously with an increase of the mass flow. The shape of the curve is typical for operating this kind of thrusters with noble gases [48,69,70]. The current through the coil exhibits a similar dependence on the mass flow.
Mass and electrical efficiency deduced from the experiment are plotted in Fig. 9 b) as function of the mass flow. The largest obtained mass efficiency η m , which is calculated as ratio of the mass flow of ions leaving the thruster to the mass flow injected into the thruster, is 71.8 %, the smallest mass efficiency is 43.4 %. Another measure for the thruster's efficiency is the electrical efficiency η e . In the case of a RF ion thruster, its common definition is the ratio of the power applied to the screen grid to the total input power required for operating the thruster. The total power is mainly determined by the sum of the RFG's input power and the power provided to the screen grid. Other contributions to the total power balance, such as the power supplying the acceleration grid (≤ 0.1 W) are negligible. In the covered range of the performance mapping, the obtained electrical efficiency is 39.4 % for the smallest and 48.2 % for the largest injected mass flow.

Comparison and discussion of the results
The performance curve obtained by the experiment serves as benchmark for the three global models. Table 2 lists the properties of the RIT 4 that are necessary to reproduce the results of the 0D model. A comparison of the results is shown in Fig. 10. The 3D model shows a very good agreement with the experimental data, as shown in Fig. 10 a). The 2D axi-symmetric as well as the 0D model do not take losses in the RFG into account. Therefore deviations to the RFG's input power are expected. Comparing the RFG's output power values derived with the three models at the same operational points a general trend is given over the entire range of the performance-mapping shown in Fig. 10 b). The 2D axi-symmetric model yielded the smallest output power of the RFG, while the 0D model estimated the largest output power. The output power values determined by the 3D model are situated between the estimates of the other models. However, at low mass flows the agreement is better between the 3D and the 0D model, whereas at large flows, a better agreement of the 3D model with the 2D axi-symmetric model is given.
In Fig. 10 c) the 2D model and the 3D model estimated coil currents close to the actual value, unless the mass flows were smaller than about 65 μg s −1 . In this range of lowest flows, both models slightly overestimate this current. The 0D model generally strongly underestimates the current transmitted by the coil. A likely reason is the very simplistic description of the coil using the approximation of long coils, which is quite inaccurate, since the length of the coil is shorter than its diameter. The 2D axi-symmetric model describes the coil as N Coil parallel conductor loops and solves the associated EM-field. Therefore, the coil's helical structure is not fully considered. The 3D model fully accounts for the coil's helical shape and solves the associated EM fields. In addition, the assumption  of a homogeneous electron density profile generally leads to an underestimation of the RF current transmitted by the coil in comparison to the holistic profile assumed by both the 3D model and the 2D axi-symmetric model [48].
In order to understand the difference in the results of the three models, as a first step, the amount of power deposited in several loss channels is studied. At first, the losses in the discharge vessel will be split into two categories: First, the losses inside the plasma volume (mainly excitation and ionization) and, second, losses by absorption of ions and electrons on the inner walls of the discharge vessel and the screen grid. The results are shown in Fig. 11. All three models show that the losses by impact of electrons and ions on the thruster's walls exceed those by processes inside the plasma volume at low input mass flows. Increasing the flow, a crossing point occurs, where losses in the volume start to become dominant. However, while for both the 3D and 2D model this crossing point is observed at a mass flow of 84 μg s −1 , in the case of the 0D model this point is already reached at a mass flow of 68 μg s −1 . Overall, the 0D model seems to overestimate the losses in the plasma volume. The 2D model estimated the smallest losses over the entire range of the simulation, for losses in the plasma volume as well as losses on the thruster's walls. To gain a better understanding of this behavior, one has to consider as well the plasma parameters and must look into the losses inside the plasma volume in more detail.
The calculated plasma parameters are depicted in Fig. 12. Raising the mass flow leads to an increase of both neutral and electron density and to a decrease of the electron temperature. The neutral density is in the order of a few 10 19 m −3 , while the electron density is smaller by about two orders of magnitude. The electron temperature of a few electron volts is typical for inductively coupled low temperature discharges [12]. The neutral density shows a very good agreement between all models over the entire range of mass flows studied. The electron temperature estimated by the 0D model is significantly lower compared with the 2D axi-symmetric and the 3D model; which are both in very good agreement. We guess that this discrepancy is caused by the difference in the assumed electron density profiles, as discussed in detail by Reeh et al. in a recent publication [48]. The electron density in the radial center of the discharge is very similar for the 3D model and the 0D model. The mean electron density estimated by the 0D model is significantly larger. For a mass flow of 75 μg s −1 , the mean electron density determined from the electron density profile calculated by the 3D model is about 82 % of the value in the center of  Plasma parameters calculated by the three models: a) Neutral gas density n 0 , and bulk electron density n + in the center of the discharge as function of the injected mass flow and b) electron energy k B T e /e. While the electron density determined by the 0D model is assumed to be radially homogeneous, for both the 2D and 3D model, the electron density profile is maximal in the center of the discharge, as shown by the insert the discharge. The reason for the slight discrepancy of the electron density between the 2D axi-symmetric and the 3D model is not known, yet. A possible reason is the use of a different cross section database. The radial profile of the electron density determined by the 2D axial-symmetric and the 3D model is identical, as shown in the insert of Fig. 12 a).
If one considers the spatially averaged electron density, the 0D model estimates a larger electron density compared to both other models. This can partially explain the larger estimated losses in the plasma volume, since all losses in the discharge volume are proportional to the product of electron density and neutral density. Figure 13 shows the power that is used for the ionization of xenon to achieve the designated screen current of 30 mA and the power which is lost due to excitation of the neutral atoms. The power required for ionization is virtually independent of the used mass flow. The power invested for ionization is similar for all models, although the 0D model estimated a slightly larger power of about 1 W − 1.5 W. As shown by all models, the power lost due to excitation rises with increasing mass flow. However, the absolute values determined by the models disagree by a significant margin.
At the largest mass flow operating the thruster, of about 94 μg s −1 , the difference is approximately 4 W between the 0D and both the 2D model and 3D model. At smaller flows, and therefore smaller neutral gas densities, the discrepancy becomes slightly less pronounced, but does not disappear. The difference can be explained by the different mean electron densities and electron temperatures. The small difference of the losses due to excitation by the 2D and 3D model can be explained solely by the difference in the estimated electron density, since the electron temperature and therefore, the rate coefficient for excitation, is identical. As a second series of simulations, the extracted beam current was varied while keeping the neutral gas density constant. In order to keep the neutral gas density constant while increasing the beam current, the injected mass flow is increased by the same amount as the equivalent additional mass flow leaving the discharge chamber through the apertures of the screen grid. For singly charged xenon ions, a current of 1 mA is equivalent to a mass flow of 1.36 μg s −1 . For the smallest beam current of 10 mA a mass flow 24.4 μg s −1 was chosen. The advantage of keeping the neutral gas density constant is that the electron temperature is constant as well.
This yields a rather simple situation, facilitating the discussion correlation of beam current I b and the ion transmission factor β i . Since the electron temperature and neutral gas density are essentially fixed, the beam current is solely dependent on β i and the electron density at the sheath edge. The simulations were performed for currents ranging from 10 mA to 32.5 mA, the results are given in Fig. 14.
For both significant larger and smaller beam currents than in the simulated range, the beam current leaving the thruster's grid system is no longer equal to the current transmitted through the screen grid, because large grid currents occur as ions impinge on acceleration or deceleration grids due to a detuning of the ion optics. The detuning occurs because the plasma meniscus depends on the plasma parameters. Such operational points have to be avoided while operating those thrusters and are, therefore, of less interest. As depicted in Fig. 14, an increasing beam current leads to an increase of the required electron density, while the ion transmission factor decreases. The averaged transmission coefficients calculated by all global models and their respective ion optics simulation show the same trend. The ion transmission coefficient β i decreases with an increasing beam current. However, especially for beam currents of 25 mA and above, the values derived by the three models deviate significantly. The difference between the 0D model and the other two models can be explained by different plasma parameters arising due to the model assumptions at the same operational points, both electron temperature and the mean electron density. The electron temperatures are virtually constant for each model. The values of 6.5 eV for the 2D and 3D model are in good agreement, while the electron temperature calculated by the 0D model is considerable lower by about 0.8 eV.
In contrast, the electron density calculated by the 0D model and the 3D model appear to be in very good agreement in Fig. 14 b), whilst that derived with the 2D model is somewhat lower. However, one has to keep in mind that the electron density in the radial center is plotted, and not the mean value, which is significantly smaller for the 3D model.
The reasons for the discrepancy of the ion transmission coefficients β i of the 2D axisymmetric and the 3D model are not well understood yet. While the calculated ion transmission coefficient is larger for the 2D axi-symmetric-model, the electron density calculated by the 3D model is larger by essentially the same factor at each point of operation. Therefore, the product of n + and β i is the same for both models, yielding the same beam current, since the electron temperature is virtually identical, as previously shown in the analysis of the performance mapping.
Computational time An important feature of global models is the rather short computational time required for modeling a performance curve. Using the 0D model, the performance curve shown can be calculated in less than half an hour including the simulation of the ion optics, which requires the majority of the computational time. If one calculates the performance mapping using a constant ion transmission coefficient, the required time is just 15 s. However, to obtain the required ion transmission coefficient at the desired current leaving the discharge, a simulation of a few minutes is required first. The simulation time was achieved with a office computer with 16 GiB RAM and a Intel © i5 hexacore CPU of the ninth generation. Despite the use of a CPU with six cores, most calculations are done using one core only. The 2D axi-symmetric-model and the 3D model require several hours to one day to simulate the ion optics as base for the simulation of a performance mapping. An exhaustive analysis of the ion optics over a wide range of ion densities, electron temperatures and, therefore, beam currents can require simulation times from several days to a few weeks. The remaining calculating time for the performance-mapping is comparably short, for the shown results about one hour. The longer simulation time compared to the 0D model is mainly caused by the required time for simulating the electromagnetic fields. The determination of the neutral gas transmission factor using SPARTA requires only few tens of minutes and must be performed only once for a given grid system. The hardware used for running the simulations of the 2D axisymmetric model contains a fairly up to date Intel © i7 quad core processor with 16 GiB RAM. The four cores are not used in parallel for all sub models, partially only one core is used. The 3D model simulation was executed on a workstation with two 8 core/16 thread CPUs (Intel © Xeon E5-2650) with 256 GiB RAM. Regarding the ions optics simulation the integration of the trajectories with the gather process as well as solving the electrical field is simulated in parallel on 32 threads while the remaining tasks, such as the scattering, are implemented on one thread. All further calculations of the remaining sub models are done using one core.

Conclusion
All global models presented provide at least a satisfactory prediction of a performance curve of the RIT-4. In particular, the 3D model can predict the necessary RFG power to operate the RIT-4 very accurately. However, in difference to the 2D axi-symmetric model and the 3D model, the 0D model cannot predict the coils current with sufficient accuracy. We guess, that this is caused by both the assumption of a long coil and the electron density profile used. Further studies concerning this issue have to be conducted. The calculated plasma parameters and losses show the qualitatively expected trends in all models. The neutral gas density is determined in the same fashion by all models. The electron temperature shows a very good agreement between two of the three models. The deviation of the third model, the 0D model, confirms the influence of the electron density profile on the determination of the electron temperature predicted by Reeh et al. [48]. However, there are significant differences in the determined electron densities, the cause of which is only partially known. The aforementioned deviations of the plasma parameters are mirrored in the different distribution of power losses such as the absorption of ions and electrons by the discharge vessel walls or due to excitation of electronic states of the xenon atom. In a second comparison, the beam current was varied. The results confirmed the differences of the plasma parameters. It was shown, that the discrepancy of the calculated electron density between the 2D axi-symmetric and the 3D model can be directly related to the different determined ion transmission coefficients. However, the reason causing the differences is not clear yet.
Future activities should, among other things, aim at a more accurate comparison between simulation and experiment. Of particular importance is a measurement of the radial electron density distribution; since in the context of this publication the influence on the different power dissipation channels in the discharge vessel was clearly shown. Furthermore, a further runtime optimization of the ion-optics simulations appears to be essential since in each of the presented models they accounted for the largest fraction of the total computational time. All models are already able to deliver useful information for