IFM Nano Thruster performance studied by experiments and numerical simulations

Field emission electric propulsion thrusters are characterised by their low thrust range, which makes them ideal for the precise control of spacecraft. Decisive for precise control are the properties of the thruster ion beam, which includes the beam divergence angle and the thrust vector. The analysis of these properties is also necessary in order to be able to estimate the interactions of the beam with components of the spacecraft. Due to such interactions, solar panels or electrical instruments on board the spacecraft could be damaged by sputtering effects. The spatial ion current density and energy distribution of a test crown emitter beam, with different specifications compared to the IFM Nano thruster, were examined experimentally with a diagnostic system, including Faraday cups and a retarding potential analyser. In addition to the analysis of the beam profile of an emitting crown, a single emitting needle was analysed. Based on these experimental analyses, an ion trajectory simulation model was developed to determine the theoretical ion current density distribution. This model includes the properties of a liquid metal ion source, where the ion trajectories start from their point of origin, the so-called Taylor cone jet cap. The benchmark of the model shows that the thrust vector and divergence angle correspond to the experimental results and shows the identical calculations for different thruster parameters, like emission current and electrode voltage. The simulation allows for the optimisation of existing and novel thruster geometries in terms of performance, reliability and longevity.


Introduction
Liquid metal ion source (LMIS) technology takes advantage of a Taylor cone, which is formed at a sharp needle tip, because Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. of its strong electric fields and liquid metal surface tension. The LMIS finds application, for example, in material preparation, microscopy and ion beam lithography, as well as in space, such as in spacecraft potential control or propulsion systems [1]. The IFM (Indium FEEP Multiemitter) Nano Thruster, developed at FOTEC and commercialised by ENPULSION, is based on this ion thruster technology, where liquid indium is used as propellant [2]. The thruster emitter consists of a porous tungsten crown composed of 28 needles, wetted with liquid indium. The other main components are the ring extractor and the housing structure, as represented in figure 1.  Research on LMIS technology started in the early 1960s when Taylor [3] analysed the shape of a liquid cone. He calculated that the cone has a half angle of 49.3 • at the static equilibrium between the electric field forces and the liquid metal surface tension. The cone structure is explained following the schematic drawing in figure 2. The Taylor cone half angle θ T is material independent [3] and also independent of the tip radius r tip [4]. Gomer developed one of the first theories on LMISs [5]. He showed that the space charge at the Taylor cone apex has a great influence on the cone shape even at low ion emission currents. Kingham and Swanson developed a theory to describe the Taylor cone shape including liquid flow and space charge effects [6], which was improved by Ljepojevic in 1992 [7]. In 1994, Praprotnik et al [8] took experimental measurements and Kingham and Ljepojevic made theoretical calculations on the Taylor cone half angle and the jet length depending on the emission current. Praprotnik observed the dynamic behaviour of an LMIS for the first time in situ at a tungsten tip coated with liquid indium in a high-voltage electron microscope during emission. He found that there is a linear dependency of the Taylor cone half angle θ T on the ion emission current I em due to the change of the field strength: with θ T0 = 51.1 • and dθ dI em = 0.298 For emission currents I em >10 µA, an extension of the Taylor cone jet occurs. The jet length increases under the influence of the electrostatic field to provide an equilibrium between the electric field forces and the liquid metal surface tension [9], [10]. Mair and Forbes [11] analysed a linear dependence of the jet length to the emission current: where γ is the surface tension, e/m the charge-to-mass ratio of the ions and V em the emitter voltage. Mair also noted that the increasing jet length with increasing emission current is part of the reason for the microdroplet emission [9]. Another consequence of increasing emission current is the increasing probability that multi Taylor cones are formed at one needle tip [8].
In 1989, Hornsey undertook an analysis of the energy distribution of a gallium LMIS [12]. He observed that an increasing emission current leads to energy broadening. One reason for the broadening effect is the longitudinal Taylor cone jet oscillation in the GHz range. These oscillation amplitudes are small (≈10% of the jet length) and do not influence the emitted current. The energy spectrum results from ions emitted at different phases of the Taylor cone jet oscillation. Hornsey also described the energy distribution as almost Gaussian and he assumed that the Taylor cone shape has a direct influence on the energy spread of the emitted ions. In addition, the energy broadening occurs because of Coulomb interactions of the ions in the beam, which is known as the Boersch effect [13]. Mair [14] measured an increase of the energy broadening (5-35 eV) and a decrease of the energy shift (∼5 eV) with increasing emission current (1-50 µA).
Beam diagnostic measurements of a capillary single-needle indium FEEP thruster, the predecessor of the IFM Nano Thruster, were made in the early 2000s. These measurements were carried out with wire probes and a Langmuir probe [15,16]. Tajmar et al observed that the shape of the beam is similar to a cosine distribution at lower thrust values and closer to a Gaussian distribution at higher thrust values. Based on these measurements, Vasiljevich made numerical simulations of ion trajectories and calculated the beam divergence. Here the ions started at 50 µm distance from the needle tip [17].
The analysis of the thruster beam shape and behaviour is of considerable interest for increasing the thrust efficiency and avoiding interactions with spacecraft components. In this paper, a detailed analysis of a test crown emitter ion beam  with simulation and experimental methods will be presented. Therein, the properties of an entire emitting crown and a single emitting needle will be analysed. The beam properties of experimental and simulation results will be compared. The shape of the ion current density profile, including the divergence angle, and the ion energy distribution, including energy shifting and broadening, will be analysed.

Experimental setup
The thruster used for this analysis is a research model from FOTEC. The structure and emitter of this model are similar to the IFM Nano Thruster, but differs with regards to the emitter specifications such as the needle shape. For analysing the test crown emitter beam, a beam diagnostic system was developed at FOTEC in 2018 [18] and [19]. The diagnostic system is located in FOTEC's largest vacuum facility LIFET4, which has a length of 3 m and a diameter of 2.2 m. The diagnostic system consists of a remotely controlled semi-circle rotating arm equipped with 23 Faraday cups (FCs) and a retarding potential analyser (RPA). The arrangement of the Faraday cups and the resulting mesh are shown in figure 3 and figure 4. The cups are arranged at different intervals, 4 • , 8 • and 10 • in the peripheral direction. At the centre of the arm is a horizontal curved extension, equipped with an RPA at the ϕ = −10 • position. The diagnostic arm can be rotated stepwise from −80 • to +80 • in a distance of 95 cm around the thruster emitter centre point. All measurements were performed at a vacuum facility pressure of 2 × 10 −7 mbar.

Ion energy distribution
The RPA was used to analyse the ion energy distribution of the test crown emitter beam. It consists of four grid electrodes and a collector. The first grid is grounded and sets up a potential gradient against subsequent grids. Due to the negative potential (−30 V) of the second grid, electrons cannot enter from the outside and influence the collector signal. At the retarding electrode, a positive voltage V ret is swept from 0 V up to 10 kV to discriminate between ions of different energies. Only ions with an energy higher than eV ret can pass through the retarding electrode and reach the collector. The suppressor, a second negative electrode (−35 V), prevents secondary electrons, which are generated by ion bombardment, from reaching the collector. This electrode provides a higher resolution of the detected ion current. The last part of the RPA stack is the ion collector, which measures the current with respect to the retarding electrode voltage. A more detailed description can be found in [18]. Figure 5 represents an RPA measurement at the thruster emission current of 37 µA and emitter voltage of 6 kV with 1 of 28 needles firing at RPA position ϕ = 0 • . The measured ion current distribution during the voltage sweep at the retarding electrode is indicated in blue. The energy distribution f (eV ret ) represented in red results from the derivation of the ion current distribution [20]: where e is the elementary charge, I c is the total ion current at the collector and V ret is the bias voltage of the retarding electrode.
To be able to compute the derivative of the ion current density in figure 5, the Savitzky-Golay filter was applied with a window size of 20 and a polynomial order of 2. The ion energy distribution is characterised with a Gaussian fit, as Hornsey [12] described. This allows for the energy shift and the width of the distribution to be analysed. In figure 5 it can be seen that the energy is distributed around 6 kV, matching the applied thruster emitter voltage with a shift of +82 eV. This energy shift is expected and will be explained in section 4.1.

Ion current density distribution
The FCs were used to map the ion current density distribution of the test crown emitter beam. During the ion current density measurements, the arm was moved from −80 • to +80 • in steps of 1 • . This step size and arrangement of the FCs result in the mesh shown in figure 4. Here, r is the constant distance of 95 cm of the FCs, ϕ is the diagnostic arm position and θ is the FC position. Each cell is filled with the density measured in its centre.
The FC repeller electrode was biased negatively (−15 V) to prevent thermionic electrons or slow charge exchange ions from entering the cup and to recapture secondary electrons that are released by ion bombardment. With a frequency of 1.5 kHz, the FC collector current was measured successively 3 . In order to reduce the noise of the measurement signal, it was measured 1000 times for a single-needle emitter and 10 times for the emitting crown. This results in a measurement time for the complete scan of 1 h for the single needle and 15 min for the crown emitter. During this period, the profile remained stable and no thermal drift was noticeable, which was verified by multiple measurements.
The data was collected on the virtual hemisphere, which is visualised in figure 6 as a projection onto a 2D plane. The figure shows an example of the ion current density distribution at a thruster emission current of 37 µA and an emitter voltage of 6 kV with 1 of 28 needles firing. The density profile shows a plateau with sharp edges. Within the region bounded by the dashed circle, 95% of the total ion current is measured, which corresponds to the beam divergence angle. This is a common  definition of the beam divergence angle used to compare electric propulsion thrusters [21]. The cross in the centre of the circle indicates the thrust vector.

Ion trajectory simulation model
For calculating the ion beam trajectories, a simulation model was set up with COMSOL Multiphysics. The system consisting of the test crown emitter geometry and detector sphere is located inside a vacuum cylinder as shown in figure 7. The ion trajectory colour represents the particles' kinetic energy, where an emitter voltage of 2 kV and an extractor voltage of −8.8 kV are applied. The particles are accelerated up to 7 kV and subsequently decelerated back to their initial potential of 2 kV. A half sphere with a 95 cm radius and with the emitter needles in the centre represents the detector where the properties of the simulated ions are recorded.
The simulation model setup was carried out based on lowemission currents to reduce the complexity of the dynamics occurring at the tip of an LMIS. For lower emission currents (<10 µA), the following can be assumed: • no Taylor cone jet extension exists [8]; • no microdroplet emission occurs [22] and [23]; • no current fluctuations and no multiple Taylor cones on the needle tip appear [24].
An ion trajectory analysis was carried out including an iteration procedure between calculating the electric field and the resulting space charge due to the ion trajectories. The main part of the simulation model is the inlet of the ions, which represents the radius of the Taylor cone jet cap r jet , as indicated in orange in figure 8. Kingham and Swanson described the shape of the Taylor cone jet cap as hemispherical [6], which is why the jet cap in the simulation model has a curvature radius of r cap . The physical diameter of the jet is around 5 nm, but the effective ion optical diameter is around 50 nm because of Coulomb interactions between the emitted ions [25]. Typically, the test crown emitter needle tips have a radius r tip between 1 and 10 µm [26]. The tip of the needle used for the measurements presented in this paper was scanned with an electron microscope (figure 8) and a radius of 2 µm was estimated. This shape was taken over to the simulation geometry of the needle and is indicated in red. As described in the introduction, the Taylor cone half angle is independent of the needle tip radius, which means that the beam divergence angle is not affected. According to [8], the Taylor cone half angle and the Taylor cone jet length are linear depending on the emission current (equations (1) and (2)). This dependency was simplified only by changing the Taylor cone half angle linearly, which is why the following pre-factors in equation (4) differ from equation (1). With this, the following function emerged by calibration for the simulation model: The calibration of this equation was based on two experimental measuring points: point 1, I em = 9 µA, V em = 2 kV, and point 2, I em = 37 µA, V em = 2 kV, which will be explained in section 4.2.1. All simulation results are based on this equation (4). For calculating the ion trajectories, the COMSOL particle tracing module was used. Here a Lagrangian description of a problem by solving ordinary differential equations using Newton's law of motion is applied [27], whereby Newton's law of motion is solved with the indium ions' particle mass m ion and forces acting on the particles ⃗ F. The forces are composed of the external electric field ⃗ F = q ⃗ E, which is computed with the finite element method, and composed of the space charge density by solving the Poisson equation.
This calculation is carried out by the particle-field interaction, which is a mesh-dependent iteration process. The emitter current is given that results in particles (super-particles); this represents a number of charged particles per unit time. The space charge density ρ(⃗ r) is assigned to the mesh elements the super-particle passed. This means there are no direct interactions between the ions, but the ion trajectories are influenced by the electric field of the particles of previous iterations. A quadratic Lagrangian shape function is used to define the variation of the electric field in each element. The calculation alternates between electric field and particle trajectories with a certain number of iterations, including the previous solutions.
The indium ions are initialised with a density distribution accordingly to the spherical surface at the Taylor cone jet: where r jet = r cap · cos (θ T ) is the radius of the inlet surface. A pre-factor of 1.3 was chosen to prevent zero particles from starting at the edge of the jet surface. The number of 1 was added to shift the parabolic distribution, resulting in the maximum density being reached in the centre of the surface. To match the experiment and simulation, equation (5) was used for all simulation results. The ions start with zero kinetic energy and with potential energy corresponding to the emitter voltage V em . For this reason, a high space charge density is formed in the area of the inlet. Therefore, the mesh elements in the area of the jet have a supporting role. At the inlet surface, a mesh consisting of triangles with a side length of 0.5-10 nm is used as represented in figure 9. Around the inlet, a sphere composite of tetrahedrons is located (tip mesh). The next larger mesh (beam mesh) has the shape of a cone because the ions spread in this area. Compared to previous simulation models (e.g. [17]), the space charge density at the point of origin is included.   Table 1 shows the experimental, geometric and numerical parameters used in the simulation model. The listed numerical parameters were converged when setting up the simulation model in order to achieve the best possible agreement with the experiment.
Timesteps are defined with a given tolerance by the timedependent solver. With this, the solver can adjust the size of the steps according to the change in the field strength.

Comparison of experimental and simulation results
The beam properties of the test crown emitter were measured at a distance of 95 cm, both by experimental measurements and simulations. The energy and the ion current density distribution were analysed.

Ion energy distribution
The energy distributions during the experimental investigations were measured for one needle firing with an emission current of 37 µA at different emitter voltages, varying from 2 to 7 kV, as shown in figure 10. The higher the emitter voltage V em , the higher the number of collected ions. With higher  V em at constant emission current, the beam divergence angle decreases and the ion current density of the beam profile increases. This leads to an increased number of ions reaching the RPA collector at ϕ = 0 • , θ = 0 • . Furthermore, the higher the V em , the more widened the energy distribution and the greater it is shifted to higher energies. The measurement points are fitted with a Gaussian distribution. The detailed results are given in table 2. ∆x c indicates the shift of the energy distribution centre and ∆E the full width at half maximum (FWHM) of the Gaussian distributions. In the V em range of 2 to 7 kV, an energy shift variation of 62 eV and a broadening of 80 eV are observed. The centre of the Gaussian distribution x c in table 2 is about 1.3% higher than the applied thruster emitter voltage V em . One explanation for this energy shift and broadening is a potential dimple in the centre of the RPA retarding electrode grid apertures, as represented in figure 11. The colour scale shows the electric potential within the vacuum area, and the grid with a retarding potential of 7 kV is represented in white. The percentage distribution of the potential stays the same for the different retarding voltages applied. This creates a maximum potential dimple of 4.2% through which the ions with lower energy can pass. Another explanation of the energy shift was described by Knauer [28], where the Coulomb potential in front of the Taylor cone increases the energy of the ions.
When looking at the width of the energy distribution in table 2, increasing energy broadening ∆E with increasing emitter voltage V em is determined. Hornsey stated that the energy distribution results from ions emitted at different phases of the longitudinal Taylor cone jet oscillation [12]. Considering equation (2), the jet length increases with V em and thus also ∆E, which explains the energy broadening trend with V em .
The ions starting from the needle tip have, besides their main longitudinal components, transversal components of motion, which are very low compared to the former. Because of the high space charge density at the tip, energy redistribution takes place, which leads to an energy broadening effect of ≈5 eV for currents <2 µA [13]. For higher currents the broadening increases rapidly [25]. Table 3 shows the experimental Gaussian fit results for 6 kV emitter voltage and different emission currents. In the I em range from 37 to 148 µA, an energy shift of 16 eV and an FWHM variation of 40 eV are determined. Similar behaviour as measured by Mair [14] is observed. An increasing emission current leads to an increase in the energy broadening and a decrease in the energy shift. Additionally, the ratio of the emission current to FWHM also moves toward saturation.
With increasing current, the production of droplets increases. Droplets are produced very close to the Taylor cone [29]. If a droplet decays due to ion bombardment close to the tip, the energy distribution is broadened and there are more ions with lower energy.
In the simulation model, the ions have the same kinetic energy as the emitter voltage with a numerical distribution of 1%. The effects responsible for energy broadening and shifting are not fully covered with the developed simulation model.

Ion current density distribution
In this section, the experimental and simulated beam profiles of one emitting needle and an emitter crown are compared. Therein, the behaviour for different thruster emission current and electrode potential ranges are discussed. For higher emission currents, microdroplet emission takes place, which is why this emission range was not covered by the simulation model. The following investigations were carried out for emitter voltages of 2-7 kV. An emission current of 1 mA was used for the crown. For the single-needle measurements, the emission current was scaled down ( 1 27 mA), which consequently does not correspond to the emission current of an entire crown.

Emission of a single needle.
The first beam profile analyses were carried out on a single emitting needle. This needle was located in the top left of the crown as shown in figure 12(a) left. The corresponding beam profile on the right side shows a shift to the top left according to the needle position. Here the experimentally measured ion current density distribution at an emission current of 9 µA is presented. It is characterised by a circular distribution with a sharp edge.
The corresponding beam profile resulting from the simulation model is represented in figure 12 (b). A circular distribution with a shift to the top left and sharp edges is also observed. When comparing the ion current density from the simulation model ( figure 12(b)) to the experimental results ( figure 12(a)), an overall higher density is determined. In the experiment, the resolution is affected by the geometry of the FCs and with this by sputtering effects occurring at the entrance aperture. Ions  hitting the negative graphite repeller electrode release electrons (electron emission yield ≈0. , of which 40% will be registered inside the cup and 60% will leave the cup 4 . These secondary electrons reduce the measured current, whereas in the simulation model no electrons are included. The profiles of the experiment and simulation show similar behaviour for different thruster parameters. This is shown when comparing figure 12 with figure 13. In both figures, the same emission current of 9 µA is applied, but different electrode potentials. The lower the emitter voltage, and with this, larger extractor voltage, the more the beam is widened. The divergence angle increases by a factor of 1.04 ± 0.06, while the extractor voltage drops by a factor of −(1.30 ± 0.27). This experimental trend is reproduced with the simulation model with almost perfect agreement.
In order to analyse the structure of the ion current density distribution more precisely, a cut through the beam profiles for different emitter voltages and an emission current of 37 µA is represented in figure 14, for both the experiment and simulation. Since with an emission current of 37 µA, the thrust vector is located at θ = 8 • , the cut through the beam profiles is at this angle. A steep increase occurs at ∼−50 • and a steep decrease at ∼ 20 • . The density within the plateaus increases towards the centre. The higher the extractor voltage V ex , the more the profile is expanded and is directed towards the nearest surface of the extractor electrode (−ϕ). This becomes distinct due to the increased distance between the profiles at the rise on the left side compared to the relegation on the right side. The density of the plateau increases towards the centre for both the experiment and simulation, and the amount of the rise depends on the emitter voltage, which confirms that the beam is more focused with a higher emitter voltage. While the shape of the simulated and measured profile shows excellent agreement, the absolute value differs by a factor of 1.4. This is probably due to secondary electrons as discussed above.
Compared to the measurements of Tajmar [15], no cosine distribution for one emitting needle is observed. Swanson [30] also measured that the angular ion current density distribution is rather flat near the 0 • axis. In Tajmar's measurements, a wire probe was used, which has no acceptance angle compared to the FCs used here, which have an acceptance angle of 46 • . Therefore, it is possible to resolve the sharp profile edges.
In order to investigate the behaviour of the beam profiles in more detail, the divergence angle α div was analysed for different thruster parameters. The divergence angle is the opening angle of the beam cone around the thrust vector, which contains 95% of the total current density. In figure 15, the experimental and simulated divergence angles are plotted against the emitter voltage for different emission currents. There is a close-to-linear dependency of the emission current on the divergence angle, as well as of the emitter voltage on the divergence angle. This is determined in both cases, for the simulation and experiment. This means the higher the emitter voltage (lower extractor voltage respectively) and the lower the emission current, the more focused the ion beam. The effect, α div ∝ V ex , occurs because with higher extraction voltage, the ion trajectories are pulled outwards towards the extractor. The second effect, α div ∝ I em , is explained with the Coulomb potential in front of the Taylor cone, which occurs by the beam itself. Due to the small emitting area, a high current density occurs in front of the Taylor cone. This leads to Coulomb interactions of the ions and with this to an increasing space charge. The Coulomb potential formed in front of the tip is influenced by the Taylor cone shape, and the latter by the emitted current.
The higher the emitted current, the higher the Coulomb potential, as shown in figure 16. Furthermore, the Taylor cone angle decreases, which has an influence on the shape of the Coulomb potential, which fans out. At the same time, this leads to a deflection of the ion trajectories, which causes a broadening and defocusing of the beam.
The divergence angle α div can be determined with the input parameters emission current I em and emitter voltage V em in the range below 37 µA and 7 kV. This angle is independent of the location of the Taylor cone on the needle tip.

Emission of an emitter crown.
In order to investigate the behaviour of the entire thruster ion beam, diagnostic measurements were carried out on an emitting crown. For reasons related to its manufacturing, 1 of the 28 needles emitted weakly; this needle was located on the left side of the crown, as marked in figure 17(a) left. The corresponding beam profile on the right side shows a slight shift to the right due to the weakly emitting needle. Here the experimentally measured ion current density distribution at an emission current of 1 mA is presented. In comparison to the single-needle profile, no sharp edges can be observed due to the overlapping of the individual circular profiles. The corresponding beam profile resulting from the simulation model is represented in figure 17(b). Here an overall higher current in the simulation than in the experiment was also measured, which can again be explained by the limitation of the FCs. The profiles show a similar shape with increasing ion current density towards the centre. The experimental beam divergence angle in figure 17 is 63 • compared to 49 • in the simulation. The simulation model is idealised, which means every Taylor cone has the same shape and is directed exactly vertically upwards. In addition, in the simulation model, the same current value of 1 27 mA is emitted from each needle tip. This cannot be influenced in the experiment, where only the regulation of the entire emission current is possible. This implies that some needles emit lower and higher currents.
In order to analyse the structure of the ion current distribution more precisely, a cut through the beam profiles at θ = 0 • for different emitter voltages is represented in figure 18 for both the experiment and simulation. The results can be compared with the results of one emitting needle in figure 14 (top), because of the same emitted current of 37 µA on average per needle. The lower the emitter voltage and the higher the extractor voltage, the more the beam profile flattens, which leads to a defocussing effect. This is the same trend as analysed for the single emitter. In comparison to other electric propulsion thrusters, the beam profile cannot be described with a Gaussian distribution. For higher emitter voltages it is a conical distribution with a distinctive tip in the profile centre. With decreasing voltage, this tip is getting flatter. The simulation results in figure 18 (bottom) show the same behaviour but an overall more focused beam profile than in the experiment.

Conclusion
Our detailed analysis of a test crown emitter beam provides an overview of the beam profile behaviour for a single emitting needle and an entire emitting crown in different emission regimes. Beam profile and ion energy distribution measurements were carried out. It should be noted that the emitter used in this work was manufactured to different specifications than the IFM Nano Thruster. The experimental energy distribution measurements have shown that with increasing emission current, an increase of the energy broadening occurs.
The COMSOL Multiphysics particle tracing module was used to simulate the ion beam trajectories. As an input parameter, the linear dependency of the Taylor cone angle on the emission current was included in the model. The model was calibrated using two experimental measurement points. Based on this, all other values were calculated without fitting. In the model, the space charge was calculated based on the slow ions in front of the nanometre-sized Taylor cone jet, so that the ions can start from their point of origin. The beam profile, including the divergence angle and density distribution of a single emitting needle, can be simulated in a low-emission current area. The experiment and simulation have shown a linear dependence of the divergence angle on the emission current and on the thruster electrode voltages. With increasing emission current and also with increasing extractor voltage, the divergence angle increases, which leads to a defocusing effect of the beam. Hereby the simulation model has shown that the inclusion of the space charge, and with this the Coulomb interaction, plays an important role in ion trajectory calculation. The model has access to parameters that are hard to determine experimentally, e.g. the space charge and the Coulomb potential, due to the ion beam itself. Due to the excellent benchmark of the model, it will be used in the near future to optimise the focusing of the thruster beam by adjusting the extraction geometry.