Similarity scaling‐application and limits for high‐efficiency‐multistage‐plasma‐thruster particle‐in‐cell modelling

To suit a wide variety of space mission profiles, different designs of ion thrusters were developed, such as the High‐Efficiency‐Multistage‐Plasma thrusters (HEMP‐T). In the past, the optimization of ion thrusters was a difficult and time‐consuming process and evolved experimentally. Because the construction of new designs is expensive, cheaper methods for optimization were sought‐after. Computer‐based simulations are a cheap and useful method towards predictive modelling. The physics in HEMP‐T requires a kinetic model. The Particle‐in‐Cell (PIC) method delivers self‐consistent solutions for the plasmas of ion thrusters, but it is limited by the high amount of computing time required to study a specific system. Therefore, it is not suited to explore a wide operational and design space. An approach to decrease computing time is self‐similarity scaling schemes, which can be derived from the kinetic equations. One specific self‐similarity scheme is investigated quantitatively in this work for selected HEMP‐Ts, using PIC simulations. The possible application of the scaling is explained and the limits of this approach are derived.


F I G U R E 1 Schematic setup of the HEMP-T-DP1
use them for full predictive modelling, since all relevant time and space scales have to be resolved, namely, the electron plasma frequency and Debye length, respectively. This leads to large domain sizes and hence high computational demands. These practical run time limits are a problem and code optimization is necessary. Analysing the Boltzmann equation, a self-similarity scaling scheme can be derived, [5] which can then be used to decrease the amount of computing time needed for PIC. The aim of this work is to evaluate the similarity scaling scheme for PIC simulations of HEMP-Ts and to study limits and deviations caused by the scaling, quantitatively. The specific system used for this is the HEMP-T DP1 thruster, [6] which is described in the next section. Afterwards, the similarity scaling scheme is derived and discussed. With a miniaturized DP1 system (the TDP1) similar to the μHEMP-T, [7] the validity of the analytically derived similarity scaling scheme is shown. This verification is accomplished through various TDP1 simulations with different scaling factors and analysis of the deviations of key plasma properties. Then, the results are investigated from PIC simulations of the DP1 thruster with different similarity scaling factors. Since the DP1 is larger than the TDP1, higher similarity scaling factors have to be applied due to run time restrictions. These illustrate the limits of the self-similarity scaling approach. Finally, a summary and an outlook are given.

THE HIGH-EFFICIENCY-MULTISTAGE-PLASMA-THRUSTER
The High-Efficiency-Multistage-Plasma-Thruster (HEMP-T) is an ion thruster developed by THALES Deutschland GmbH. [1] A magnetic multi-cusp arrangement surrounds the thruster channel of the HEMP-T with subsequent axially magnetized permanent magnet rings with opposite polarization. The discharge channel is coated with a dielectric with a high sputtering threshold, such as boron nitride. Anode and neutral gas inlet are placed at the upstream end of the channel. Usually xenon is used as operating gas, due its high mass and low reactivity. Outside of the channel an electron source, representing the cathode, serves as a neutralizer for the emitted ion beam and an electron source for the plasma discharge in the channel. The high magnetic field strength of up to ∼0.5 T leads to a magnetization of the electrons, resulting in Larmor radii much smaller than the channel diameter. Ions are not magnetized because of their higher mass. The cusp-like structure of the magnetic field resembles a magnetic mirror towards the dielectric wall in some channel regions. It reflects electrons impinging on the dielectric and thus traps them inside the thruster channel between the cusps. This trapping of electrons in the discharge channel results in high ionization rates. Another advantage of the cusp structure is the low plasma wall contact, which is limited to the cusps, resulting in very low sputtering yield rates. This leads to a long lifetime and a high efficiency of the thruster. The magnetic field enhances axial transport of the electrons in the channel which results in a flat potential structure. At the channel exit, the potential drops to vacuum level, accelerating the emitted ions and hence generating the thrust. This work considers the thruster design HEMP-T-DP1 (short DP1), developed by Koch et al. [6] Figure 1 shows a sketch of the DP1. The discharge channel is 64.0 mm long and has a radius of 15.1 mm. The magnetic circuit consists of three magnet rings with different lengths, resulting in a three-cusp magnetic field configuration. In the channel, the cusps are located at axial positions of about 15 and 30 mm, respectively. The last cusp is located at the channel exit. At the bottom of the discharge channel (left side in Figure 1), the anode is located with an applied voltage in a range of 100-2,000 V. The neutral gas inlet is placed at the centre of the anode. Xenon is used as the propellant gas with injection fluxes between 10-30 sccm.
To better understand the scaling of the physics of the low temperature plasma in the HEMP-T, the fundamental kinetic equations are used to derive a scaling scheme in the next section.

SCALING LAWS IN PLASMA PHYSICS
A plasma is described by the distribution function for each particle species s in space, velocity, and time f s (r s , v s , t). The time evolution for one species is given by the Boltzmann equation with the force F s acting on a particle with the mass m s and the collision term In electromagnetic systems, F s is the Lorentz force. If the particle density is small, the collision term can be neglected and the Boltzmann equation simplifies to the Vlasov equation. [8] The dynamics of a collisionless plasma is then fully described by the coupled system of the Vlasov and Maxwell equations.
Two systems are called similar if all physical quantities of both systems obey a scaling law which relates them through a similarity scaling factor , e.g. the length x = ⋅x. The tilde indicates the quantities of the scaled system in all following descriptions. The scaling laws for the physical quantities are derived from the invariants of the Maxwell-Vlasov system for collisional plasmas, which were derived by Lacina et al. [9] The five invariants have to be conserved by the scaling scheme. Here, E is the electric field, B is the magnetic field, t is the time, m is the particle mass, q is the particle charge, v is the particle velocity, and n is the particle number density for each plasma species s. C 1 and C 2 are the invariants of trajectories in the electric and magnetic field. C 3 and C 4 describe the invariants of the self-induced electric and magnetic fields generated by internal currents. C 5 describes the similarity of a nonstationary process. Of special interest in a magnetized low-temperature plasma is the invariant C 2 , also called the Hall parameter which relates the gyro motion r g to the system dimensions L. Additionally, the scaling of the trajectories in electric fields C 1 and the evolution of trajectories in time C 5 have to be conserved. In HEMP-Ts magnetic fields induced by the plasma currents are small compared to the externally applied magnetic fields. Therefore, an electrostatic approximation is valid for the simulation of ion thrusters. As a result, the invariants C 3 and C 4 can be neglected. For ion thrusters, collisions have to be considered and therefore it is important to look for further invariants of the system. Applying binary collisions to the Boltzmann equation one can derive a sixth invariant, [9] where is the mean free path. The expression from Equation (6) represents the Knudsen number Kn, relating the mean free path to the system dimension x. This results in a total of four different scaling constants (C 1 , C 2 , C 5 , C 6 ), which have to remain invariant in a similarity scaling scheme for the low-temperature plasmas in ion thrusters. To derive the scaling approach most valid for PIC simulations, the PIC model is explained first in the next section.

THE PIC METHOD
The PIC method is used for simulating low-temperature plasmas with non-Maxwellian distribution functions. [10,11] Due to the axial symmetry and the negligible influence of self-induced magnetic fields in the DP1, the used PIC-MCC code is an axisymmetric, electrostatic 2D3v PIC code operating in cylindrical coordinates (r,z). All three dimensions of the velocity space (v r , v z , v t ) are resolved to preserve energy and momentum in the collision algorithms. In regions where the electric field is perpendicular to the magnetic field, like at the exit of the HEMP-T, an ExB drift occurs in plasmas. This drift of the electrons leads to a Hall current in the azimuthal direction, resulting in density and electric field gradients which then lead to instabilities. This greatly increases the cross-field mobility of the electrons at the exit. For Hall-Effect-Thrusters, the influence of the anomalous transport has been thoroughly studied. [12] In the used 2D simulation the azimuthal coordinate is only resolved in the velocity space, therefore neglecting the possible density and electric field instabilities. This anomalous transport is implemented as a random walk in velocity space acting as a rotation of the velocity vector of the electrons. It directly scales with the local electric and magnetic field, using a Bohm-scaling. The implementation of the anomalous transport in the present simulation was derived by Kalentev et al. [13] using 3D simulations.
In the simulations presented here, neutral xenon particles Xe, electrons e, and singly Xe + and doubly charged xenon ions Xe 2+ are included. The density of doubly charged ions is only a small fraction (<10%) of the singly charged ions. Because the fraction of higher charge states less than 1% combined to the total number of ions, they are neglected. [14] The Coulomb collisions are handled as described by Takizuka and Abe. [15] Also included in this simulation are direct single and double e-Xe impact ionization, single e-Xe + impact ionization, integral elastic Xe + -Xe collisions (including charge exchange and momentum transfer), and integral elastic and inelastic e-Xe collisions, [13] using a binary collision model [16,17] which utilizes experimentally measured collision cross section. [18] Further information on the PIC model used here can be found in previous publications. [19,20]

Similarity Scaling in PIC
PIC has to resolve the Debye length, which is several orders of magnitude smaller than the system dimensions, leading to large domain sizes in terms of the total number of cells required to cover the domain. Large domains result in high computing times for solving the Poisson equation and a high number of simulated particles, which again increases the work load for particle motion and collisions. To reduce the computational effort of PIC simulations, the following self-similarity scaling scheme is introduced, based on the plasma invariants discussed Section 3. Keeping the plasma density and hence the spatial resolution of PIC constant, which are determined by the Debye length, it is useful to reduce the system size by the similarity scaling factor , This results in a scaling of the old domain length to x =x with > 1 and thus in a reduction of the number of cells in the simulation with n , where n is the number of spatial dimensions of the PIC model. With the scaling constant of time evolution C 5 from Equation (4) it follows that the time is scaled as t =t, to keep the particle velocity constant. The scaling of the domain size in combination with the resolution of the Debye length leads to a scaling of the cell volume relative to the system size with the scaling factor ∼ n . To keep the number density constant, the number of super-particles in each cell has to be scaled, too. The scaling of the particle numbers is then applied to the super-particle factor, increasing the amount of particles represented by one super-particle and keeping the total amount of particles to be tracked low.
In HEMP-Ts, the motion of particles is determined by the magnetic and electric fields. The scaling of both fields follows from the derived plasma constants of Equation (2) and the scaling of the system length Due to the electric field being scaled with −1 , the electric potential Φ remains nonscaled. To preserve ionization and other collisions, which are vital for thruster operation, the plasma invariant for binary collisions C 6 from equation Equation (6) has to be considered. The mean free path is then scaled as =̃through scaling the collision cross sections TA B L E 1 The most important quantities in the self-similarity scaling scheme. A tilde indicates the quantities of the down scaled system. is the similarity scaling factor

Quantity Scaling law
Length scale Number density n = ñ Super particle flux per unit area Γ sp =Γ sp for each collision process represented in the model, with the averaged velocity ⟨v⟩ and taking into account the other scaling laws. It is important to mention, that the numeric implementation of the anomalous transport is represented by a Bohm-like transport. The Bohm transport scales with the magnetic field, which scales with the self-similarity scaling factor. The anomalous transport is implemented according to the self-similarity scaling. A general overview of the resulting scaling scheme is given in Table 1.
One should note, that the presented scaling scheme is exact for all plasma volume processes, including particle trajectories, binary collisions, and even fluxes. But because the system size shrinks while the plasma density remains constant, the ratio of the Debye length Db to the size of the scaled system changes The increased ratio leads to an overestimation of charge separation, which occurs on the scale of the Debye length. Especially, in regions with smaller Debye lengths, such as the plasma sheath or the exit region and plume of ion thrusters, the scaling overestimates the charge separation. When extrapolating the results from the down-scaled system to the real one, this leads to regions where quasi-neutrality is seemingly violated.
The plasma boundary sheath with a typical dimension of a few Debye lengths increases relative to system size L. This leads to a higher sheath to plasma volume ratio for higher scaling factors. The increased charge separation leads to longer decay lengths of the potential, increasing the influence of the potential boundary conditions for larger similarity scaling factors. As long as the sheath is small compared to the plasma volume the similarity scaling is an accurate model for the discharge, but has an upper bound of applicable similarity scaling factors. For HEMP-Ts this influence manifests mainly in the plume and at the thruster exit. Additionally, if three-body processes become important the scaling lacks accuracy, because only binary collisions are covered in this self-similarity scaling Scheme. [9] This is usually not important for ion thrusters, but can be dominant in fusion edge plasmas. [21] The applicability of a scaling model was already shown for Hall thrusters by Taccogna et al. [5] The scaling presented here is equivalent to the scaling of the vacuum permittivity 0 [22] or a scaling of the plasma density n e , as was shown by Lacina. [9] Multiple groups applied such a scaling scheme already in the past without discussing their limits, for example, to simulate spoke frequencies in a Penning discharge [23,24] or to gain physical insights in Hall effect thrusters. [22,25,26] However, a quantitative study of the influence of the scaling factor on the results is missing. The aim of this work is to present a quantitative analysis of the influence of the self-similarity scaling scheme for the example of the HEMP-T system, including channel and plume solutions.

RESULTS
To study the influence of the scaling factor on the system, the following strategy is pursued. At first, the PIC simulation setup of the DP1 thruster is described to introduce the considered thruster design. A verification study of the similarity scaling is then carried out for a smaller system, the Tiny-DP1 (TDP1). The simulation results from the TDP1 are additionally used to describe the general characteristics of HEMP-Ts. Afterwards, the larger DP1 is used to study quantitative limits of the similarity scheme.

Setup of the simulation
To cover the channel and the near-exit region of the DP1 the simulation domain measures 35 × 100 mm (r,z). With an expected charged particle density of 10 13 cm −3 and an electron temperature of 10 eV, the spatial resolution Δr = 3.7 × 10 −4 cm and temporal resolution Δt = 2.8 × 10 −13 s are chosen. The target particle density is represented by six super-particles at the symmetry axis and increases linearly with radial distance. The anode potential is set to 500 V and the metal and the top domain boundary are set to ground potential, while the right domain boundary is set to an electric field of zero in axial direction to simulate a vacuum transition downstream. At the anode the neutral gas inlet is set to a mass flow rate of 15 sccm. An electron source is located at the right domain boundary, replacing the electron current from the neutralizer of the real system. The injected electron current is an input parameter for the PIC simulations.
In ion thruster physics, the operation mode is defined by the anode current, beam current, and thrust. Higher electron source currents lead to higher plasma densities and thereby increase anode and beam current. Therefore, the strength of the electron source current is the parameter to generate different operational modes in thruster simulations, assuming a constant anode potential and neutral gas mass flow rate. For the simulations of the DP1, an electron current of 360 mA is applied. The simulation domain is set up similar to Figure 1. The similarity scaling scheme introduced in Section 3 is applied to scale down the size of the simulated system and therefore decrease computing time. With increasing scaling factor, the number of cells in the domain decreases according to Table 2.

5.2
Verification of the self-similarity scaling scheme with the TDP1 To verify the similarity scaling scheme, a new thruster design, in the following referred to as the TDP1, is introduced to be able to study small similarity scaling factors. It is geometrically similar to the DP1, but the geometric dimensions are reduced. The channel length of the TDP1 is 21.33 mm with a radius of 5.033 mm. Its magnetic field topology is similar to the magnetic field of the DP1. For the same plasma parameters as the DP1, the total number of cells in the simulation is reduced. The simulation with a low scaling factor represents the reference case from which the similarity scaling factor can be incrementally increased, which is used to verify the similarity scaling. The plasma solution characteristics should not vary strongly across the different scaling factors for the verification. In the Figures A1 and A2, the resulting electron and ion density distributions are shown for different scaling factors relative to the reference case in the range between ref = 1… 1.6. For discussion of the HEMP-T characteristics, the simulation of the reference case of the TDP1 thruster with ref = 1 is considered here. The results of the PIC simulation of the TDP1 thruster model show the same general behaviour as other HEMP-Ts, namely, a flat potential in the acceleration channel, a potential drop at the exit and the trapping of the electrons between the magnetic cusps. The electric potential, as shown in Figure A3 is flat in the thruster channel due to the high electron mobility, in particular close to the axis.
Typical for HEMP-Ts is the cusp structure leading to low plasma wall contact, which is limited to the cusp region. The electrons in the thruster are trapped between the cusps, which act as magnetic mirrors, enhancing the ionization of the neutral xenon gas, which can be seen in the electron density distribution in Figure A1. In addition, the ion energies at the F I G U R E 2 Behaviour of global values of the HEMP-T, namely the maximum electrical field at the exit (blue), the global ionization (red), the anode current (yellow) and the thrust (black), with rising scaling factor. The values are normalized by the values from ref = 1 cusps are so low that they are below the sputter threshold and produce very low sputtering rates at the channel wall. [2] The magnetization of plume electrons at the exit cusp, leads to an an additional plasma in front of the exit. Due to the magnetic field structure the electrons follow the magnetic field lines towards the axis and into the discharge channel. In the density gap between the plasma in the exit cusp and the acceleration channel the potential drops and the electrons are accelerated by the electric field upstream into the discharge channel. The ion density distribution in Figure A2 is similar to the electrons in the acceleration channel. At the exit, they are accelerated by the electric field and are then ejected into the plume to generate the thrust. All these results are typical for HEMP-Ts. [2,27] Because the TDP1 shows these typical characteristics of a HEMP-T, the aim is now to verify the self-similarity scaling scheme. Because the scheme is exact for volume processes and wall fluxes, only little impact on the plasma distribution in quasi-neutral regions of the discharge is expected, i.e. in the acceleration channel.
The evolution of the system is now shown for increasing similarity scaling factors in relation to the reference case. A very important constraint of the self-similarity scheme results from the additional influence of charge separation producing electric fields. This is determined by the Debye length, which changes relative to the system size and hence causes deviations from the reference case. For low similarity scaling factors the ratio of Debye length to system size does not increase much, leading to a conservation of the plasma solution by the similarity scaling. This is shown by comparing some of the bulk parameters at different scaling factors, shown in Figure 2. Afterwards, the maximum of the electric field is determined at the exit. For the global ionization, the ionization collision rates are integrated over the whole channel volume. The anode current is obtained by integrating the axial particle flux crossing the left domain boundary and converting them into an electric current by multiplying with the correspondent particle charge and divide by the time step. The thrust is calculated by summation of the axial momentum of the emitted ions. For better statistics the data is averaged over 10 6 time steps of a converged simulation, which is equal to 1.68 × 10 −5 s in the scaled system. The averaging time covers all relevant transport processes in the system.
In Figure 2, it can be seen, that the maximum electric field at the thruster exit, which is responsible for the acceleration of emitted ions, stays nearly constant over the considered range of scaling factors. The same behaviour is observed for the global ionization, the anode current and the thrust of the TDP1. These selected parameters are representative for the global thruster solution and show that the similarity scaling conserves the solution for small scaling factors.
To further study the similarity scaling scheme, the time averaged axial electron density profile in the channel at r = 4 mm is shown in Figure 3 for different scaling factors. It shows the increase of density at the cusp, the electron density of trapped electrons between the cusps and the density drop at the exit of the thruster. All density profiles have a similar shape, with variations of the total values of up to 30%, which can be explained by the slight variations of the global ionization, see Figure 2.
The channel solution in front of the potential drop is similar for the applied scaling factors as long as the influence of the charge separation in the plasma sheath remains small. This limit is violated for large scaling factors when the radial sheath approaches the radial diameter of the channel, similar to systems such as the HEMP-T. [7] In small systems with an already high ratio of Debye length to system size this limit is reached with already comparably small scaling factors. Therefore, the similarity scaling factor is not the correct parameter to decide if realistic solutions are possible. Another parameter to be considered is the ratio of Debye length to scaled system sizẽ It is important to note that the limit, where one considers the scaling to produce too strong deviations, depends on the question to be addressed. In qualitative studies deviations of up to 50% can be sometimes acceptable, whereas for specific design studies even deviations of 10%, e.g. of thrust, can be too large. In general, deviations of several 10% can be still considered acceptable, given other uncertainties in designs. [28] The scaling is limited by simulations, where the sheath covers the whole radius of the thruster channel, which leads to unexpected and physically not correct behaviour. Some of these limits of the similarity scaling can be found in the work from Brandt et al., [7] where a down-scaled HEMP-T is discussed. Now, the influence of larger scaling factors on a large system, but in a similar range of̃is studied.

Quantitative self-similarity scaling study for the DP1
The larger DP1 system is now used to study the influence of the applied similarity scaling on the PIC results. The size of the DP1 leads to run-time restrictions for the PIC calculations. Therefore, high similarity scaling factors have to be applied to stay within reasonable computing times of weeks. Full PIC simulations of the DP1 for scaling factors in the range from = 60 to = 100 were performed, which results in different sized domains as listed in Table 2. For completeness, the density distributions and potential solutions can be found in the Appendix A. Since the self-similarity scheme is exact for volume processes and wall fluxes, only little impact on the plasma distribution in quasi-neutral regions of the discharge is expected, especially in the acceleration channel in front of the potential drop.
In the plasma density distribution in Figure B1 can be seen, that most of the channel distribution between anode and potential drop is unchanged by the applied scaling Scheme. As expected, a difference appears close to the exit region, where the plasma contracts further towards the inner channel with a higher scaling factor. This behaviour is generated due to the different decay length of the grounded domain boundary conditions in addition to the overestimated charge separation because of the increased Debye length relative to the system size with increased similarity scaling factor. With higher scaling factors the grounded potential from the metal at the thruster exit extends further, since the relative distance in terms of Debye lengths between the boundary and the plasma is smaller. This pushes the potential drop further into the channel, resulting in a decreased plasma bulk volume.
In Figure 4 a comparison of the electron density distribution in the channel is shown for two different scaling factors of = 60 and = 100. Most of the channel distribution is unchanged by the applied scaling scheme. As expected, a difference appears close to the exit region, where the plasma contracts further towards the inner channel with a higher scaling factor.
The decrease of the plasma volume in the exit region, see Figure B1, leads to reduced ionization in the plume and decreased ionization at the exit cusp with increasing similarity scaling factors. The dependence of the ionization current on the similarity scaling factor is shown in Figure 5. Less ionization leads to a lower plasma source for the channel plasma, which decreases the effective plasma volume and density. In front of the anode the electron density decreases with rising scaling factor. This is due to the already discussed decrease of plasma volume and density. Because of less electrons in the discharge, less electrons reach the anode region and lead to a lower electron density in the pre-anode region. In comparison to the main ionization areas in the DP1, namely, between the two cusps in the middle (z = 15-30 mm) and between the second cusp and the exit (z = 30-50), the changes in the anode region have only little impact on the overall plasma solution. Other plasma parameters in the thruster channel, such as thermal and excitation losses, remain nearly constant with the variation of the scaling factor. This demonstrates again the physical accuracy of the chosen self-similarity scaling scheme for volume process and wall flux dominated plasmas, where the charge separation on the Debye scale relative to the system size is not dominating the behaviour with increasing scaling factor the mean emission angle drops ( Figure B2). Now, the changes in performance parameters for different scaling factors are investigated. Thruster performance is mainly defined by the anode I a and beam current I b , the thrust T and the beam efficiency beam = I B ∕eṁ source at a given neutral gas mass flow rateṁ source . Especially, the ion angular current distribution is important and accessible experimentally using a Retarding-Potential-Analyser (RPA). [2] The mean emission angles of the ion angular current distributions for different scaling factors are shown in Figure 5. With increasing scaling factor the mean emission angle drops ( Figure B2). This leads to a higher thrust contribution of the ions due to a smaller ratio of radial to axial momentum because only the latter contributes to the thrust. To explain the drop of the mean emission angle, one has to take into account the potential solution.
In Figure 6, axial potential profiles at the axis and radial potential profiles at the exit are shown. The axial potential at the exit (z = 50 mm) is nearly at anode potential and drops downstream from the exit until it reaches vacuum potential. For reference see also Figure B3, where the full potential solutions are shown. With increasing scaling factor the plasma contracts, due to the increasing influence of the grounded metal at the thruster exit as discussed before. In addition, the increased charge separation creates lower electric fields at the exit, leading to reduced ionization processes in the exit region. This trend is visible also in the global ionization current, see Figure 5, which also decreases with rising scaling factor. The contracted plasma and the decreased potential in the exit region, as shown in Figure 6, draw the potential drop further into the thruster. This contracted potential profile decreases the radial potential drop ( Figure 6) outside of the thruster exit leading to lower radial electric fields. Therefore, the ions receive less radial momentum which leads to a decrease of the main emission angle of the exiting ions. The change of the ion angular beam current leads to the effect that the same electron source current does not represent the same operating state for different scaling factors because of the change in the plasma density distribution. In Table 3, the most important parameters defining the operation state calculated from the PIC simulations of the DP1 are listed for different scaling factors. In summary, the beam divergence is reduced with increasing scaling factors, increasing the contribution to the thrust of the ions due to lower emission angles. The decreasing plasma bulk volume and hence the decrease of the ion density volume compensates the decreasing emission angle, so that the thrust of the DP1 stays nearly constant for different scaling factors, except for the highest scaling factor = 100, where it drops sharply. The reduced ionization leads to a reduced beam current and therefore to a decreasing beam efficiency beam of the thruster. For practical applications, the limit for scaling factors is important, which is determined by the ratio of the Debye length to the scaled system size. As already mentioned, this range changes with the system size of the considered device, because higher scaling factors can be applied to bigger thrusters. To compare different system sizes, the following physical As discussed above, the biggest changes to the plasma properties are expected close to the exit, where the influence of the space charge should be the strongest. To identify this effect, the channel is divided into different characteristic volumes. In each of Figures 7-9, the upper figure represents axial cuts, representing the volume at the axis (r = 0-0.3 cm), in the centre of the discharge (r = 0.3-1.2 cm) and the sheath towards the dielectric (r = 1.2-1.4 cm). The figures at the bottom depict radially oriented slices, representing the anode area (z = 0-1 cm), the first cusp (z = 1-2 cm), the second cusp (z = 2-4.5 cm), the exit region (z = 4.5-6 cm) and the plume region (z = 6-10 cm).
From Figures 7-9 one can see that the scale of the relative deviation is different, showing that the axial cuts show in general lesser variations than the radial cuts. But with increasing ratios̃, the deviations increase for all considered properties. Especially the radial cuts which cover the exit region show strong deviations with increasing ratio of Debye length to the scaled system size.
When comparing the different plots, it is apparent that the rate of the increase of the deviation rises sharply, when the ratio of Debye length to scaled system size gets bigger than 0.01. For lower ratios the deviations from the simulation at a ratio of 0.006, are lower than 20-30%, which usually can be considered to be sufficient for design studies. [28] As a consequence for the DP1 and the TDP1, the scaled channel length should be at least 100 times larger than the Debye length, while the channel radius should be at least 20 times larger than the Debye length. This can be considered as a general guideline, at least for the considered thruster design.
The results in this section demonstrate that the similarity scaling factor influences the results of the PIC simulations mostly in the thruster exit region. One expects the simulation with a scaling factor of 1 to be the best representation In the present PIC studies, it shows that with increasing similarity factor the mean emission angle decreases, which in the contrary leads to higher emission angles at lower scaling factors. But in comparable HEMP-Ts, the experiment measures small emission angles ∼30 • . [2] This creates a mismatch of PIC and experiment, because PIC simulations with similarity scaling factors already tend to overestimate the emission angle of HEMP-Ts compared to the experiment. [29] As shown in Figure 5, the ionization rate at the exit changes with scaling factor. The position of the potential drop reacts to this and affects the ion dynamics being accelerated into the plume. This changes the emission angle. In general, the exit region and the plume, which are influenced by the potential and charge separation at the exit have to be handled with care. Here, a combination with other methods, such as hybrid codes, can overcome these problems. [30] It is expected that for lower similarity scaling factors the boundary conditions will have less impact on the plasma solution, as shown for the TDP1.

CONCLUSIONS AND OUTLOOK
In this work, principles and limits of a self-similarity scaling scheme applied to PIC simulations of HEMP-Ts simulations are shown. The aim of this scheme is the reduction of the computing time of PIC simulations, at the cost of an increased ratio of Debye scale to system size. For the investigation of the influence of the scaling factor, PIC results for two different thruster designs and different scaling factors were considered. First, the scaling scheme was verified with the smaller thruster design of the TDP1 thruster, which allowed the application of low scaling factors. It was found that the scaling scheme delivers nearly identical results in the plasma acceleration channel, including fluxes and ionization currents. However, simulations of a larger thruster, the DP1, demonstrated the influence of the potential boundary conditions and the overestimation of the charge separation for higher similarity scaling factors, which are responsible for the reduction of the plasma bulk volume in the exit region of the discharge channel. The contracted plasma bulk leads to a decreased potential at the thruster exit and hence to a decreased emission angle, while generating comparable amounts of thrust. It was shown that the scaling factor has a non-negligible impact on the beam divergence and thrust. Therefore, if predictive simulations are required, the simulation results have to be benchmarked against unscaled simulations. It was found that the ratio of electron Debye length to scaled system size can be used to estimate practical limits for the application of the similarity scaling. To generate similar potential solutions the down-scaled channel length should be at least 100 times larger than the Debye length, while the down-scaled channel radius should be at least 20 times larger the Debye length. For a better scaling of particle densities and volume ionization processes even stricter ratios have to be chosen. Similarity scaling is a powerful tool to reduce computation time. Simulations with high scaling factors can be taken as a new starting point for further simulations with lower scaling factors, since the initializing phase of PIC simulations until convergence is reached is quite time and resource consuming otherwise.
To overcome the changes in the plume region, different solutions can be applied in the future. A hybrid solution, [30] where the electrons are treated as a collisionless fluid, could be a better representation since it is not influenced by charge separation on the Debye scale. However, the disadvantage of the hybrid method is the loss of the full kinetic information. Another method could be a more application-oriented approach, such as a Multi-Objective Design Optimization (MDO). [31,32] Here, experimental measurements can be combined with simulation characteristics for design optimization. As shown in this work and in ref., [5] the solution of the acceleration channel in front of the potential drop does not change for different scaling factors, and thus, these plasma properties can be used for MDO, which is based on a zero-dimensional power balance equation system. With this approach a large design space can be explored, where the lack of physical accuracy is minimized with the insight gained from PIC, which promises a fast and robust tool for thruster optimization in the future.