Simulations of xenon beam ions emitted from electric thrusters in Earth’s magnetosphere

This study sought to characterize the behavior of exhausted electric thruster xenon ions in the near-Earth magnetospheric environment as functions of various trajectory and particle attributes, neglecting effects of electric fields, plasma waves and particle collisions. This was done via simulation using the AeroTracer program, a software tool which computes ion trajectories within the magnetosphere by applying an adaptive step-size Runge-Kutta technique to the fully relativistic Lorentz equation. Over 3,800 independent simulations were performed, with variables including release position, release energy and direction, ion charge, and orbital phase. Initial release altitude was a major driver in determining whether the ion eventually fell to Earth (“Below Minimum Altitude” or BMA), remained trapped by the simulation’s end (“Maximum Number of Steps” or MNS), or traveled beyond the magnetosphere (“Lost to Space” or LTS). Ions expelled at the highest altitudes investigated - 60,000 km and above - almost invariably were lost to space. Like altitude, increasing inclination and energy were important factors that reduced trapping, affecting the outcome probabilities. Higher charge state produced strong improvement of trapping capability. Effects of orbital phase, day of year and solar cycle phase were also apparent. A transition region was found in the 20,000 km to 60,000 km altitude range, within which the sensitivity of outcomes to parameter variation increased. The ordered sequence MNS> BMA> LTS was found to be consistent with decreasing confinement capability, and it was manifested consistently as parameters were varied.


Introduction
Electric thrusters are utilized for a substantial fraction of space vehicle propulsion needs, including near-Earth orbit transfer, on-orbit station keeping, repositioning and end-oflife disposal, as well as deep space operations.With a propellant usage efficiency greater than the chemical propellant alternative, electric propulsion devices typically rely on electricity and magnetism to energize and accelerate propellants, creating thrust.Usage has increased over the past several decades [1,2], with explosive growth occurring recently [3,4].
Electric thrusters can operate with a variety of propellants (such as noble gases or liquid metals), depending on their design and mission considerations [5][6][7].While Xe + has been a standard propellant for years, Kr + became a major player with the advent of the Starlink megaconstellation, and more chemically interesting species such as iodine are contenders as well.Tens of thousands of Starlink satellites are planned, each relying on Kr + Hall thruster electric propulsion, and population of the constellation has been progressing.OneWeb has populated a constellation of >500 satellites using Xe + propulsion.Smaller numbers of satellites using more powerful Xe + ion thrusters, and relatively low power iodine, FEEP and electrospray devices are in service as well.As a result of this and other programs, vehicles with electric thrusters on board already dominate the satellite population larger than cubesats and will increasingly do so in the coming years.
Use of Hg + devices, despite the fact that this massive ion provides improved efficiency relative to Xe + and Kr + , and was once a standard propellant, has been greatly constrained by concerns about its toxicity.Electrospray thrusters emit ion species that can be much more massive and chemically complex than Xe + , itself a heavy ion.
Mission design influences the range of initial position coordinates and velocity vectors for exhaust particles (most commonly Xe + and Kr + ) as they exit the thruster.Spacecraft Earth orbits range from circular and very low, such as JAXA's Tsubame at 170 km, to extremely high elliptical like NASA's IBEX with its 300,000 km apogee [8,9].Electric thrusters often fire for long continuous periods, so that ions and a relatively small number of neutrals are released along a slowly changing orbital track.Each position of release on the orbital track will have a unique set of local magnetic field, latitude, longitude, altitude and Earth-Sun orientation parameters.In addition, ion thrusters typically have charge state and exit velocity distributions, significant ion beam divergence, and expel neutral propellant species as a minor constituent of the exhaust.Only ionized exhaust is considered in the present study.Similarly, while the effect of ion beam divergence on trajectories is certainly of interest, it was considered out of scope for the present study.
The study was intended to provide insight into the fate of ions exhausted by orbiting spacecraft throughout the magnetosphere, particularly those with trajectories terminating on or very near Earth's surface.This can be of concern for the more toxic propellants.In addition, better understanding of the dynamics of anthropogenic particle emissions and the influence of various release parameters was targeted.

The magnetosphere
The magnetosphere shields life on Earth from solar and cosmic particle radiation, and prevents erosion of the atmosphere.When strongly perturbed by solar activity, its protective effect can be reduced and negative impacts on satellites and some other elements of commerce and industry may be experienced.The geomagnetic field diverts supersonic solar wind particles around the planet.The direction of solar wind changes at the bow shock boundary, with a fraction of particles crossing this boundary to populate the magnetosheath region.Some go further and cross the magnetopause boundary where outwardly-directed atmospheric pressure approximately equals solar wind pressure (Earth's magnetic field is weak here so magnetic pressure is low), and enter the magnetosphere -the region within which the dominant magnetic field is Earth's.Inside the magnetopause, the plasma is predominantly of terrestrial origin, and outside it is predominantly of solar wind origin.The magnetosphere is shaped by the solar wind interaction and changes dynamically -with the magnetopause boundary able to shift in excess of 100 km/s -always adjusting to balance solar wind pressure [10].
Once an ion is expelled from a near-Earth thruster, the trajectory is strongly influenced by Earth's magnetic field.Inside the magnetosphere, ions can be trapped for long periods of time, but they can also spiral into Earth's surface and be lost, or exit through the magnetopause and no longer be significantly influenced by Earth's field.As the geomagnetic field lines are "pushed" by the solar wind, they stretch out far behind Earth to create the magnetotail which contains a reservoir of plasma and energy [11].Ions created in the ionosphere and upflowing to higher altitude constitute a second major source of magnetospheric ions as shown in Fig. 1.The main ions here are H + , He + , and O + , with additional ion species contributing [12].A range of energies and latitudes, from polar to mid-latitude, are involved.
Due to solar wind influence, rotational axis offset ( ∼ 23.4 • ) with respect to the ecliptic plane's normal vector, time dependence of the magnetic pole direction, an imperfectly shaped oblate spheroid, and a variety of other factors, Earth's magnetic field direction and magnitude at the moment of ion release from an electric thruster is complex.It deviates substantially from a perfectly symmetric dipole and depends on solar activity, day of year, time of day, latitude, longitude, altitude, orbital phase and other orbital parameters [13,14].Magnetic field magnitude and direction are far more influenced by solar wind at high altitude than at Earth's surface.Intensity contours at the surface are shown in Fig. 2. For a perfect dipole, the intensity falls off as the inverse cube of radial distance, 1/r 3 , and because the geomagnetic field strength is so strongly dependent on distance, the spacecraft altitude is a critical determinant of exhaust ion trajectories and their ultimate fate.
Earth's magnetic field largely determines the trajectory for any charged particle within the magnetosphere.These particles may be introduced into the magnetospheric environment from multiple sources, including solar wind [16], nuclear explosions [17], electric thrusters [18], and the ionosphere [12].While the overwhelming majority of solar wind ions are energetic protons, the ionosphere supplies heavier and relatively slow species Fig. 1 Charged particle flow and density diagram for near-Earth space [12].Relevant regions of the magnetosphere and their particle energies are indicated including O + , N + , and NO + to the magnetosphere at the rate of about 10 26 ions per second, and this is now believed to be the dominant ion source for the magnetosphere [12].The heaviest of this group of ions, NO + , is about 30 amu.Kr + and Xe + are much heavier, about 84 and 131 amu, respectively.Transport of ions from the ionosphere is via the polar wind, as diagrammed in Fig. 1.
The small fraction of neutrals emitted by an ion thruster has been neglected in the present study, but we can expect that xenon atoms emitted at lower altitudes will become ionized [19,20].This is because the critical ionization velocity for xenon is only 4.2 km/s, which is exceeded for spacecraft over a wide range of orbital altitudes.In addition, the neutral xenon may be ionized via charge exchange interaction with ion species such as H + .The charge exchange process is an important mechanism causing decay of hydrogen ions (independent of xenon neutrals, which are historically of low abundance), an essential component of the magnetospheric ion population [21].They are sourced from both solar wind and ionospheric sources.Xe + may have significantly longer lifetime than H + since charge modification upon collision with atmospheric species such as O, H, OH, N, N2 and their ionized forms tends to be energetically unfavorable, allowing it to persist as an ion.Any potential interactions with atmospheric particles are not included in our simulations, however.
The rate of magnetospheric ion supply from the ionosphere could be matched by 50,000 electric thrusters emitting 0.5 milligram/s of Xe + , on average, and somewhat less for Kr + .A much smaller number of very powerful thrusters could do the same.30 kW Hall thrusters will typically exhaust on the order of 100 mg/s Xe + while operating, and a suitable satellite with these devices on board could potentially act as a space ferry or service platform.If the equivalent of 24 of these 30 kW thrusters operate continuously Fig. 2 Contours of total intensity for Earth's magnetic field at the surface, in nanotesla units [15] they would exhaust about 1E23 Xe + ions per second.The exhaust rate is then about 3 orders of magnitude below the estimated ionospheric supply rate, but due to its mass and persistence, Xe + may be very long lived with low average velocity, such that it may build up sufficient density to alter magnetospheric properties as effectively as the ionospheric O + .The spatial distribution associated with thruster and ionospheric sources would be expected to differ, as well.Due to the rapid growth of electric thrusters in space, in the future heavy ion population density in some portions of the magnetosphere may be determined primarily by space vehicle exhaust.The sheer numbers of ions to be released from electric propulsion devices, and the fact that they are unlike any that normally exist in the magnetosphere, potentially leads to new phenomena.
The magnetosphere is a dynamic environment [22], and its activity during geomagnetic substorms is hazardous in several ways.Induced current flows at ground level have been known to disrupt the electrical grid, causing power outages.Other potential effects include loss of communications, loss of radio and GPS navigation capability, increased drag on satellites and damage to satellite electronics.The ions and electrons present, whether they have been trapped for a long time or newly arrived from the ionosphere or solar wind, affect the way storm conditions develop and evolve.When the magnetosphere is disturbed by heightened solar wind activity, magnetic field levels are depressed at the planet's surface, and the ionospheric species already mentioned can significantly contribute to total energy content of the ring current and geomagnetic storm and substorm dynamics [23].Some of the energy increase stems from ion acceleration in elevated electric fields, potentially resulting in significant transport of ions from one location to another.Given the complexity of these phenomena and knowledge gaps regarding many of the specifics, the build up of high ion density in one region of the magnetosphere is potentially problematic for other regions.
Due to the asymmetric nature of the geomagnetic field, particle entrapment may be temporary or permanent.Details of the magnetic field such as strength, direction of field lines, weak/strong points and anomalies may result in escape points for a given set of ion parameters.Since in this study, we ignore the effect of other influences such as plasma waves, ion beam divergence, collisional losses, magnetospheric electric fields, etc., their fate is determined by parameters at the time of injection.Since collisional losses are neglected, the effect of the atmosphere is not included in these simulations.
The present study of ion trajectories investigates the interactions between ionsreleased with energies and altitudes consistent with xenon electric thruster operation on an orbiting spacecraft -and the geomagnetic field.The methodology is amenable to other electric propulsion propellants such as argon, krypton, iodine, hydrazine and its derivatives, cesium and ionic liquids, as well as those which may be of greater environmental concern such as mercury [24].In fact, such questions of environmental effects due to electric propulsion emissions have been raised since the 1970s with investigation of cesium [25].Even xenon propellant was considered in this context, specifically for the type of orbit raising and other mission parameters for WGS spacecraft [18].The present study also explores the fate of ions injected into the magnetosphere by electric thrusters, but examines the effect of many more parameters and addresses a different set of spacecraft missions.Behaviors identified in the present study may serve as a reference for future mission designs where the fate of the expelled ions is considered.
Xenon itself poses negligible risk upon collision with atmospheric species or Earth surface, but the methodology employed in our study is not exclusive to this propellant.The same approach could be applied to more hazardous propellants such as hydrazine, monomethylhydrazine, ionic liquids, mercury, cesium and others.Xenon was chosen due to its ubiquity in the electric propulsion world, but these same methods could be applied to whatever the propellant of choice might be for an interested party.In addition, the present results may be extrapolated to other species to obtain qualitative insights regarding trajectory outcomes.

Prior work
While previous investigations have sought to better understand ionic interactions within select regions of the magnetosphere, little research has been done into the behavior of recently expelled, ionized propellants in a realistic near Earth environment.Furthermore, little work has been done to characterize them throughout multiple regions of the magnetosphere.
Past work on understanding ionic interactions within the magnetosphere has generally focused on specific regions or topics, such as the magnetotail [26][27][28] and upstream of the bow shock [29], as well as vacuum and viscous environments [30].While some studies have focused on ion dynamics calculations [26] and computer simulations [30,31], others have investigated ion distributions via global magnetohydrodynamic and ion trajectory computations [27].Each of these studies, comprehensive in their own rights, tends to draw attention to an individual magnetospheric region.
The basic details related to motion of ions in the magnetosphere have been investigated for many years.Some examples include Kress' mechanism for the trapping of solar energetic particles (SEPs) in the inner magnetosphere [16], Damiano's proposed understanding of ion gyroradius effects in kinetic Alfvén waves along auroral field lines [32], and studies that go back much further [33].None of these studies, with the exception of [18], combine an accurate magnetic field model with ions emitted under representative electric thruster conditions, or explore the effect of orbital parameters, energy, inclination, day of year and so on.The exception mentioned was specific to the Wideband Global SATCOM (WGS) spacecraft and therefore pertained to a singular mission profile and addressed few of these parameters.Computational simulations of ion motion in Earth's magnetic field usually assume a perfect dipolar field, without the deviations caused by the oblate spheroid shape, South Atlantic Anomaly, strong influence of the solar wind, etc [34,35].
Perhaps one of the first studies of xenon beam coupling with ambient proton plasma was conducted by Wang, who identified three mechanisms by which these interactions take place: particle/particle Coulomb collisions, as well as "processes [that] operate in collisionless plasmas and are sensitive functions of... the angle between the solar wind velocity and the interplanetary magnetic field" [36].
Magnetospheric current flows, especially the build up and decay of ring current flows, has received considerable attention [21].Ions in the magnetosphere have been found to affect magnetic reconnection, [12] which mediates the transfer of mass and energy from the solar wind into the magnetosphere.Magnetic reconnection leads to magnetospheric storms and itself affects ion transport.Heavier ions from the ionosphere appear to significantly increase density and facilitate reconnection.Ions also create and interact with plasma instabilities, and their behavior in the magnetosphere is very complex with much remaining to be learned.
Past research on particle trapping within the Van Allen radiation belts provides additional context for trapping in all regions of interest.A stable inner zone is comprised of ultrarelativisitic electrons and energetic positive ions, which may exist in the belt upwards of years.The outer zone is comprised of MeV-range electrons, which may be trapped on the scale of days [37,38].A transient third ring has been observed to trap high-energy (>2 MeV) electrons upwards of four weeks.

Identification and utilization of coordinate systems
Ion trajectory simulations were performed using AeroTracer software, developed by Michael McNab at The Aerospace Corporation, while varying orbital and environmental parameters.The AeroTracer program requires the components of a state vector (three position coordinates, three velocity coordinates) to define a case.While conversion of particle properties (e.g.mass, energy) and epoch were relatively straightforward, coordinate transformations were necessary for conforming to AeroTracer's inertial coordinate system.Test cases were first visualized in terms of inclination, orbital phase and altitude.The coordinates were then transformed into the cartesian, geocentric equatorial inertial (GEI) coordinate system.Only circular orbits were considered, greatly simplifying the analysis.
Referring to Fig. 3, the variable denotes the angle from the X V axis in GEI coordi- nates to the satellite position at the moment of ion release.Universal time is always midnight (0:00 UTC) when ions are released.The orbital configuration in the figure is for In Fig. 3, +X V denotes the direction from Earth to Sun at the time of the Vernal Equi- nox.Projection of the +X 0 orbital plane axis onto the equatorial plane coincides with +X V , and the Y 0 and Y V axes are aligned with the intersection of the two planes.There- fore , the right ascension of ascending node relative to +X V , is 90 • .At the Winter Sol- stice shown in Fig. 4, the equatorial plane is tilted up from the solar direction and the ecliptic plane by about 23.4 • .
The −Y 0 / − Y V axis direction therefore points above the Sun.Expressions were derived for the satellite position and velocity within its circular orbit, in the GEI coordinate system.The x and y coordinates are projections of the satellite position onto Earth's Equatorial plane.The z coordinate is satellite altitude with respect to the Equatorial plane, i is the orbit inclination and is the polar coordinate for satellite position in its orbital plane.
Derived expressions for the x, y, and z coordinates as shown in Fig. 3, are Velocity components of satellite motion projected onto the Equatorial plane are obtained from time derivatives of the position coordinates in Eq. 1.
All ion trajectory calculations were performed with universal time set to 0 (i.e.midnight) so that Greenwich and 0 • longitude always faced away from the Sun.As a result, GEI coordinates associated with 0 • longitude varied with day of year.For the Winter Solstice it was aligned with +Y 0 / + Y V and at the Vernal Equinox with −X V and −X 0 .Earth has effectively undergone a +90 • axial rotation during its travel from Winter Solstice to (1) Vernal Equinox, which can be recognized in the Aerotracer input parameters by adding 90 • to .The effect of longitudinal variation for a fixed day of year is included in Aerotracer simulations by adding the longitudinal deltas to this adjusted parameter.is equal to 90 • at the Winter Solstice and 180 • at the Vernal Equinox, for 0 • longitude, and may be calculated in general from where and rel are in degrees and rel is the satellite longitude at the point of ion release.For Aerotracer calculations, longitude increments were typically 45 • , covering 8 equally spaced positions around the satellite orbit for a given inclination.The validity of these coordinate systems and transformations was confirmed by utilizing AGI's Systems Tool Kit (STK) to produce ephemeris data of satellites of identical trajectories to those modeled.

Ion trajectory simulation processes
In order to calculate and track the motion of charged ions (such as xenon), AeroTracer applies an adaptive step-size Runge-Kutta technique to the fully relativistic Lorentz equation.Two models were used to capture the Earth's internal and external magnetic field behaviors, respectively.The first -the International Geomagnetic Reference Field (IGRF) -approximates the Earth's internal geomagnetic field "with a truncated spherical harmonic series to represent the distorted dipole nature of the field" [39].The second -Tsyganenko (T96) model -takes satellite-produced magnetic field observations to utilize a best-fit parameterization for the external magnetic field [18,40].Notably, these are static and parameterized models that do not include waves and other magnetospheric activity produced by variations in solar wind, substorm conditions, or generated from the ion beam of the thruster itself.While these wave components are typically negligible, they may become more significant when integrated over paths about the Earth.This model also neglects the free energy carried from the ion beam that may excite plasma instabilities, which may lead to scattering of beam ions through wave-particle interactions outside of the thruster exit on the order of kilometers [36].
Through AeroTracer, iterative computational simulations were made to identify and record the final destinations and general behaviors for ions as functions of orbital parameters.A flowchart of the simulation process is given in Fig. 5.In an attempt to create as wide a breadth of geocentric missions as possible, multiple orbital parameters (e.g.inclination, altitude, true anomaly) were varied and simulated within STK.The times of these simulations were intended for four different days of the year: Vernal Equinox, Summer Solstice, Autumnal Equinox, and Winter Solstice.These days were chosen to separate the Earth's position relative to the Sun in distinct yet manageable sections for simulations.
After confirming the validity of the coordinate transformation into GEI coordinates (with the use of STK), tests were run over a variety of the aforementioned orbital parameters.Numerous aspects of a test's trajectory were recorded for analysis.Firstly, the particle's orbital trajectory was saved (as a position-time text file), followed by screenshots of the visualized trajectory (within the AeroTracer interface) as shown in Fig. 6, which included notes on the ion's terminating condition.
Three separate terminating conditions were utilized within the program to stop simulations: • Below Minimum Altitude (BMA): If the calculated altitude falls below the designated "minimum altitude" criterion at any point in an ion's simulated trajectory, the simulation stops and the outcome is classified as BMA.The minimum altitude was selected to be 0 km.• Lost To Space (LTS): At any point in an ion's simulated trajectory, if the ion's position crosses the magnetopause boundary, the simulation stops and the ion is considered outside the influence of the magnetosphere.• Maximum Number of Steps (MNS): If the simulation reaches the end of its designated maximum number of steps without satisfying either the BMA or LTS conditions, the ion is still undergoing magnetic mirroring (and is trapped) when the simulation ends.For the purposes of this study (and balancing computational cost with the number of total simulations to run), the maximum number of steps was normally set to two million.
Due to the adaptive step-size Runge-Kutta technique employed by AeroTracer, these two million steps can result in a wide range of "simulated" time, depending largely on the parameters used for any given test.A strong influence of altitude on simulated time is seen in Fig. 7, where groupings of time trend upward with altitude.Step size is related to "curvature" of the motion, which depends on magnetic field strength, energy, mass, charge state and so on, all of which relates to confinement capability.Since magnetic field strength varies with the inverse cube of distance (for a dipole), the effect of altitude dominates.Where confinement is robust (low altitude, low inclination, low energy, high charge state, low mass, high B field) the step size and elapsed time per step is small.
Lifetimes for MNS cases are left unspecified.It is possible that Xe + ions at some altitudes remain trapped for years, while lifetime at low altitude is shortened due to interactions with atmospheric species.At the very least, ion trajectories will be greatly altered by these interactions.Simulations run with identical input parameters always give the same result for a given set of run parameters.Ideally, simulations would be run for the equivalent of months or years to ensure that MNS outcomes are stable for the long term, in the absence of collisions.A few of the MNS cases were run much longer than the standard two million steps, to verify insensitivity to run time.
Everything in the figure is for 300 eV, Solar Minimum, and 0 inclination, encompassing all of the days and orbital phases (see Table 5).The Winter Solstice, Day 355, has a large proportion of the outliers in the plot.
Each terminating condition is represented pictorially (up to 10,000 computed points per case) in Fig. 8. Notably, the LTS simulation in Fig. 8 represents the same test case shown in Fig. 6 to serve as a comparison for the views between the AeroTracer interface and a MATLAB-generated position plot.
Despite the GSM coordinate system's use in visualizing the ion trajectories on Aero-Tracer, GEI was actually used for AeroTracer inputs, in-program calculations, and determinations of ion final destinations.Further, within the AeroTracer program, parameters used to set terminating conditions (e.g.minimum altitude, maximum number of steps) may be updated via the Options -Tracing Control panel.

Comment on ion velocities
The equation below describes satellite speed for a circular Earth orbit: Fig. 8 Examples for each terminating condition (viewed downward in the z-direction), with colors intended to reflect the terminating condition conventions throughout the study where v S is spacecraft velocity, µ E is Earth's standard gravitational parameter (or Gm E , where G is the gravitational constant and m E is mass of Earth), R E is the Earth's radius, and h is orbital altitude.
Because the geomagnetic field rotates with the planet, ion exhaust velocity must be considered with respect to the GEI coordinate system.A xenon ion's exhaust velocity with respect to the spacecraft is: where E is the exhaust energy and m Xe is the mass of the ejected xenon ion.
The equation below describes the ground speed, V G , of an Earth-orbiting satellite with respect to the Earth's rotation about its axis: where V E is Earth's rotational velocity about its axis, T SD is the period for a sidereal day, R S is the satellite's orbital radius, and is the geographic latitude of the point of interest on the Earth's surface.
For an observer standing on Earth, ion velocities for a geosynchronous satellite (which has a ground speed of V G = 0 km/s) are equivalent in magnitude regardless of whether the spacecraft is orbit lowering (OL) or orbit raising (OR).
Similarly, because the fixed geospherical coordinate system moves in tandem with the Earth's magnetic field, the field itself remains stationary in this coordinate system.Therefore, for the sake of simplicity, cases were first developed with respect to these fixed geospherical coordinates and later converted into geocentric inertial coordinates to be more closely compatible with the AeroTracer software.
Orbit lowering and raising maneuvers for a spacecraft are visualized in Fig. 9.During orbit raising maneuvers, ions are expelled in the opposite direction of the spacecraft's velocity vector, giving them higher absolute velocity and increasing the overall magnitude of spacecraft GEI velocity.Conversely, decreasing an orbit is done by expelling ions in the same direction as the spacecraft's velocity vector, giving lower absolute velocity and ultimately reducing orbital velocity and altitude.
Equations for orbit lowering and orbit raising ion velocities, respectively, are defined as follows: For the set of parameters related to this study, velocities are plotted in Fig. 10 as a function of altitude for both of the simulated thruster energies.
Velocities for the most common parameters tested throughout the study are displayed in Table 1. (4)

Results and discussion
The set of investigated parameters for the full set of simulations performed in the present study is given in Table 2.These apply to the moment of ion release.In Table 2, Alt is the release altitude of the ion, i is the inclination orbital inclination of the "satellite" releasing the ion, DOY is the day of the year, KE is the kinetic energy of the ion, and Direction refers to the ion velocity vector at time of release.It is of interest to determine the probability that specific exhaust ion characteristics (e.g.ion charge, solar cycle time) have on likelihood of an ion's return to Earth.The influence of specific parameters was therefore evaluated.Additionally, ion exhaust energies were commonly set to 300 and 1000 eV to represent generic Xe + exhaust energy for electric thrusters of the Hall and gridded ion engine types, respectively.For convenience, it was assumed the thrusters had zero divergence and were emitted parallel or anti-parallel to the spacecraft velocity vector.
The sum of terminating condition results as a function of altitude for the tests outlined in Table 2 is displayed in Fig. 11, where the sum of all terminating conditions is 2304.The altitude of an ion's release is a highly significant driver with respect to ion termination condition.This inverse relationship between altitude and propensity to return to Earth is sensible because stronger magnetic force results in more tightly controlled ion gyration about magnetic field lines and increased confinement capability.The Lorentz force of a moving charged particle in the absence of an electrical field is where q is ion charge, v is ion velocity, and B is magnetic field strength.

General cases
Seeking to capture ion behavior within the magnetosphere and compare results for a manageable set of cases, we will focus on the group of simulations indicated in Table 3.These will be referred to as general cases.
(8) The results from simulations laid out in Table 3 were segregated according to ion energy (directly related to ion velocity) and inclination, and displayed in Fig. 12.
It is apparent that -for both low and high inclination cases -higher ion velocity (energy) promotes ions being lost to space.This happens at the expense of MNS outcomes.Similarly, for ions of a given energy, higher orbital inclination strongly promotes loss to space (LTS).Here BMA outcomes increase together with LTS outcomes, entirely at the expense of MNS.
Next, results were separated according to longitude angles at time of release as shown in Fig. 13.In low inclination and low altitude scenarios, ions terminate at Earth's surface selectively at 315 • (-45 • ) longitude as evidenced by the spike of BMA cases there for both tested energies.This longitude coincides with the South Atlantic Anomaly (SAA), which appears as a local minimum of magnetic field intensity over the Atlantic Ocean (see Fig. 2).The effect of this pronounced local minimum is to cause the altitude of ions flying over to dip, which can result in ground impingement  and BMA loss for low altitude ions.One way the dip is manifested is by higher radiation levels on the ground, since energetic radiation belt particles are passing closer to the surface there [41].The inner Van Allen belt dips to an altitude of roughly 200 km.With a longitudinal span from +270 • (west) to +40 • (east) and latitudinal span from roughly -50 • to 0 • , the 315 • longitudinal tests at the equatorial plane fall within the upper region of the SAA [42].Based on the available evidence we conclude that the BMA anomaly at 315 • is indeed associated with the SAA.Further, tests were always performed at 0:00 UTC, with 0 • longitude at midnight, so that 0 • and ±45 • longitude are least perturbed by the solar wind.Earth's magnetic field is strongest at low altitudes, but higher energy increases the gyroradius and makes ground collision more likely.For the highest energy we see that 0 • and +45 • cases also produce the BMA result, in addition to -45 • , so there is symmetry about 0 • .The high energy apparently reduces the effect of magnetic field variations on the outcome, as might be anticipated.
To gauge the influence of inclination on ion behavior for the general cases, the lower energy (300 eV and 1000 eV) data from Fig. 12 were separated by inclination and plotted against altitude to produce Fig. 14a.The same data (300 + 1000 eV at 0-degrees inclination, all days included), were fitted using a logistic function Fig. 12 Outcome frequencies for each terminating condition within the general cases, separated by energy and inclination, where the total numbers of cases for each column sums to 384 Fig. 13 Results of BMA cases across all longitude values indicate anomalous outcomes at 315 • where x is altitude and A 1 , A 2 , x 0 , and p are fitting parameters.
The same group of simulations for 60-degrees inclination fits well to a similar function, but the transition begins at much lower altitude and spans a wider range.BMA outcomes for the same group of simulations and 60-degrees inclination have a rapid transition from high to low probability between 500 and 2000 km data points, and this behavior is well fit with the same logistic functional form as for LTS outcomes.Results are shown in Fig. 14b and Table 4.The transition occurs at even lower altitude for i=0, so that data set did not produce a suitable fit.Most of the LTS and BMA results could be fit with a simple exponential function in addition to the logistic function, but with reduced quality of fit.These results indicate the robust confinement that exists for many low inclination release cases, from low altitude to geosynchronous orbit and beyond.In contrast, high inclination cases result in loss much more frequently.
Further segregating the data according to energy, Figs. 15 and 16 are obtained.With both energies included in Fig. 14a, a stark difference is observed between inclinations 0 • and 60 • , as evidenced by the separation between MNS and BMA plot lines, as well as the altitude difference for the transition to predominately LTS outcomes.This same behavior is clearly seen in Figs. 15 and 16, demonstrating that energy is here less influential than inclination, since these figures are quite similar.Still, it is observed that the LTS transition comes at slightly lower altitude for 1000 eV compared to 300 eV, and (9)   there are minor differences at 500 km altitude as well.The influence of energy is especially difficult to see inside the inner radiation belt.
In sum, Figs. 15 and 16 suggest that high inclination, high energy, and high altitude each reduce trapping efficiency and enhance BMA and/or LTS outcomes.BMA outcomes are naturally favored when ions spiral with large radial amplitude, and high energy coincides with large radius of motion as well as the ability to penetrate further along converging magnetic field lines before mirroring occurs.A high altitude release implies that ions are subject to relatively weak magnetic field, with increased probability to reach the magnetopause boundary.The parameters that favor particle trapping over the widest altitude range are shown to primarily be low energy and low inclination.

Influence of energy and day of year
From this same data set of general cases, Day of Year (DoY) analysis was done to investigate its influence on the terminating condition.Results were tabulated at 20,000 km for singly charged ions during Solar Minimum, separating out the two energies and inclinations, all eight and longitude values and the firing directions.Together with 60,000 km and GEO, 20,000 km is considered "transitional", with elevated sensitivity of outcomes to parameter variation.
While the result is always MNS for zero inclination, the difference is again pronounced between 0 • and 60 • inclinations, as many LTS cases (no BMA) occur at 60 • inclination.At this inclination, LTS cases at 1000 eV outnumber those at 300 eV by a Fig. 17 Transition Region results for both lower and higher energy singly charged ions taken at 20,000 km altitude during Solar Minimum count of 31:23.We ascribe this to the more limited ability of magnetic field to control the motion of high energy ions.Release at highly inclined orbit includes release of ions at higher latitudes, depending on orbit phase, where spacecraft may be closer to open field lines and the magnetopause boundary.These factors imply that ions can escape the magnetosphere more easily.The LTS cases at 1000 eV can be largely interpreted as an expansion of LTS cases at 300 eV, as a result of reduced confinement capability at the higher energy.Complicating this interpretation, however, is the fact that several LTS cases at 300 eV became MNS outcomes at 1000 eV -presumably Fig. 18 Average of outcome frequency on altitude taken over all tested days of year.Outcomes were summed over , inclination, energy, and direction parameter results for the general test simulations Fig. 19 Deviation of outcome frequency for each day of year from the average.Outcomes were summed over , inclination, energy, and direction parameter results for the general test simulations because the larger gyroradius caused a collision with the ground before the particle could reach the magnetopause.
Figure 17 illustrates the great complexity and variation of outcomes that is involved, even after eliminating the main variable (altitude).Energy, inclination, orbital phase, day of year, and direction of release (in spacecraft direction of motion or opposite to it) are all influential, and there are instances where an outcome case is utterly isolated from any similar outcome.It results from the complexity of the planetary magnetic field, shifting Earth-Sun geometry and dependence of the particle interactions on their physical parameters.Because of the complexity involved the simulation is essential to obtain outcomes when confinement capability is somewhere between robust (normally get MNS) and weak (normally get LTS).
To compare influence of DOY, an averaged plot is shown in Fig. 18 along with each day's respective deviation from this average plotted in Fig. 19.Summing the test results for a given DOY, regardless of other parameters, removes any dependence on , energy, direction etc., potentially revealing the influence of rotational and magnetic axis pointing and right ascension of the ascending node.Figure 19 indicates that, for altitudes at GEO and below, MNS outcomes are more probable at the Winter Solstice (Day 355) relative to the other three days, with the exception of very low altitude on Day 264 where BMA outcome probability dips.The magnetic south pole recently migrated across the Canadian arctic, and is less than 5 • away from the rotational axis for the year 2020 assumed in these Aerotracer simulations, at about 180 • longitude.It is therefore displaced a few degrees further from Greenwich, than Earth's rotational axis.Since the simulations release ions when Greenwich is at midnight, the magnetic axis is more highly inclined to the ecliptic at the summer solstice as much as it is less inclined at the winter solstice.Therefore the reduced confinement probability at summer solstice, and the approximate mirror image of Day 172 and 355 deviations across the altitude range seems reasonable.The equinoctial days (81 and 264) are approximately inverted with respect to BMA outcomes at low altitude, despite similar magnetic and rotational axis orientation relative to the ecliptic.Presumably the BMA difference drives the MNS inversion, but we have not attempted a deep dive into these data to decipher the underlying cause of the BMA outcome.At higher altitudes days 81 and 264 do behave similarly.

Investigating cases of singly and doubly charged ions
Not all exhaust ions from an electric thruster are singly-charged.In the case of xenonbased electric thrusters, Xe 2+ flux is typically about an order of magnitude lower than Xe + [43].It was therefore of interest to investigate the effect of ion charge on trajectory behavior.The cases shown in Table 5 were run for Xe + and Xe 2+ .The number of BMA, MNS, and LTS outcomes are plotted for each ion charge state in Fig. 20.LTS probability decreased for the double charge case, particular at GEO altitude and above, enhancing MNS outcomes.At 2000 km, LTS probability is essentially zero and at 100,000 km it is virtually 100% for singly-charged ions.The data points in the middle typically have LTS probabilities between these extremes and are much more sensitive to the specific conditions.Three BMA cases for doubly charged ions produce a small unexpected blip at 60,000 km, but otherwise results are consistent with the magnetic Figures 21 and 22, respectively, segregate the full set of Xe 2+ and Xe + outcomes according to whether spacecraft orbit is being lowered or raised.The plots show that the lowering cases have higher probability for LTS outcomes at and above the transitional altitude.Orbit raising outcomes reflect better magnetic confinement, especially for doubly-charged cases.The difference of Xe + MNS results between orbit raising and lowering is greater, compared to Xe 2+ , stemming from more tenuous confinement of Xe + .The effect of elevated release velocity in the inertial frame (for orbit lowering) is enhanced when trapping is weak.At 500 km, the magnetic field is far stronger than at GEO, and little to no outcome difference is observed there between orbit raising and lowering circumstances, although the doubly charged ion is more tightly confined.We also note that if Xe 2+ energy were double the Xe + energy, as is approximately true for typical electric  An alternative way to compare orbit raising and lowering is to plot charge state termination results as a ratio, as shown in Fig. 23.This may be used to gain additional insights into release velocity direction and the impact it has.Here the drastic dip in BMA at 2,000 km is due to a single case.GEO altitude is of particular interest because it is toward the low end of the transition region, and should be sensitive to release velocity direction and other parameters.A large dip in the OR/OL ratio is observed for LTS.It recovers as altitude increases, reaching unity at 100,000 km.In contrast, not all Xe 2+ are LTS even at 100,000 km and the higher LTS probability for orbit lowering makes the OR/OL ratio less than 1.At GEO and below LTS does not occur for Xe 2+ .At 60,000 km altitude Xe + MNS outcomes are heavily biased toward orbit raising, so that the OR/OL ratio is maximized, whereas it takes higher altitude to reach peak ratio for double charge ions -as expected.

Investigating solar cycle extrema cases
Maximum solar activity during the Sun's 11 year solar cycle ("Solar Max") corresponds to greater distortion of Earth's magnetic field, and the effect of this on outcomes was investigated.Simulated tests are outlined in Table 6.
Comparing solar maximum and minimum results in Fig. 24 indicates that increased solar activity promotes a slightly higher LTS outcome and simultaneously reduces particle entrapment.This is most apparent in the ascent from GEO to a 60,000 km altitude, serving as a transitional region from few to many LTS conditions, and vice versa for MNS.Overall, solar cycle phase had less effect on results compared to charge state and energy.In addition we infer that the solar cycle parameter has little influence on the low altitude magnetic field and low altitude ion trajectories.However, it should be noted that geomagnetic substorm conditions could have far greater influence on outcomes than change in the average activity level during the solar cycle.
Splitting orbit raising and lowering as shown in Figs. 25 and 26 reveals that lowering cases account for essentially the entire difference between solar min and max outcomes.This difference is again concentrated overwhelmingly in the transition region where trapping probability begins to fall and there is elevated sensitivity to particle parameters.
Figure 27 displays the raising to lowering outcome ratios to accentuate the effect this parameter has on BMA and LTS results in the transitional region.It is not until higher altitudes (i.e.GEO and above) that notable discrepancies arise.A handful of orbit lowering LTS cases from GEO (three and four during solar minimum and maximum, Fig. 28 Tabulated terminating conditions for ions released with varying values of , ion velocity direction, and release energy.This considers singly charged particles released at 0 • inclination during Solar Minimum on Day 81 respectively) are responsible for the initial dip in the OR/OL ratio.Perhaps more interesting is the effect the solar cycle has on magnetic trapping at high altitude.Trapping probability is typically higher during orbit raising compared to orbit lowering, but the difference is greatly exaggerated during a period of peak solar activity.
Again, the distinguishing characteristic of note between orbit raising and lowering is presumably ion velocity in the GEI system, which is much larger for orbit lowering.During orbit raising the LTS probability is more limited, and the magnetic field can confine ions more effectively.This capability is reduced as ion velocity increases and magnetic field strength weakens.

Energy comparisons of 300 eV, 1000eV and 39390 eV
Effects of ion energy are presented in Figs.28 and 29, tabulated and plotted according to each direction for ion release, respectively.A 39390 eV energy was also tested as an extreme case to better elucidate the trend.No effect was observed at lower or middle altitudes.Differences do arise in the transition region (GEO -60,000 km), specifically for 0, 45, and 315 degree longitudes.Since zero longitude always corresponds to midnight, these longitudes are opposite the Sun and 45 degrees in both directions.As expected, particle trapping is reduced for orbit lowering (primarily around 0 • longitude), due to the additional energy in the inertial frame.Increased energy within orbit lowering and orbit raising does exert some influence, but the subtle difference are difficult to fit within the usual paradigm.One interesting case sequence involves GEO at 45 degrees, where outcome changes from MNS at the lowest energy to LTS at 1000 eV and finally to BMA at 39390 eV.The switch back to BMA at the high energy seems surprising, but may be related to the difficulty of reflecting such an energetic particle before it penetrates to ground level.In fact, the 39390 eV energy generates far more BMA outcomes than the other energies.All of these appear at either 500 km or GEO altitude, and are part of the progression away from MNS outcomes as energy increases.At 39390 eV, 2000 and 20,000 km altitudes stand out with respect to MNS outcomes, because confinement capability has been compressed from both ends of the altitude range.This is consistent with the fact that highly energetic protons are a dominant feature of the inner radiation belt and not the outer, whereas electrons with even higher energy are featured in the outer belt.The gyroradius is proportional to the product of mass and velocity, so that protons and electrons are much easier to confine than xenon ions of comparable energy.It explains why the former particles can be trapped in the belts, despite energy exceeding 10 million eV.
Data from Fig. 28 were also plotted against the Orbit Raising/Lowering ratio as shown in Fig. 29.As before the ratio mainly deviates from unity in the mid-altitudes and above.With the OR/OL ratio greater than unity for MNS outcomes at GEO, it seems clear that orbit raising continues to favor particle entrapment.Conversely, at GEO and above, the OR/OL ratio dips below unity for BMA and LTS conditions, because ions exhausted during orbit lowering at high altitude tend to have too much energy to remain trapped.

Conclusion
This study sought to characterize the fate of ions exhausted from a xenon electric thruster on an orbiting spacecraft in the near-Earth magnetospheric environment, depending on various trajectory and particle attributes.Effects of release position (altitude, inclination, phase and direction of the initial motion), energy, charge, day of year and solar cycle phase were considered.The scope was limited to circular orbits and steady state conditions, and did not include divergence, plasma waves, collisions and other effects that may affect ion trajectories and lifetime.Each simulation result produced one of three outcomes, with probability for each outcome related to its rate of occurrence.Initial release altitude was a major driver in determining whether the ion eventually fell to Earth, was trapped indefinitely, or was considered lost to space.Ions expelled at the highest altitudes investigated -60,000 km and above -almost invariably were lost to space.Increasing inclination and energy, like altitude, were important factors that reduced trapping probability, affecting the outcome probabilities.Higher charge state produced strong improvement of trapping capability.Effects of orbital phase, day of year and solar cycle phase were also apparent.A transition region was found in the 20,000 km to 60,000 km altitude range, within which the sensitivity of outcomes to parameter variation increased.The ordered sequence MNS> BMA > LTS was found to be consistent with decreasing confinement capability, and it was manifested consistently as parameters were varied.The altitude range of high confinement capability coincides roughly with the extent of the Van Allen belts, although exhaust ions from typical electric thrusters are confined more readily than the highly energetic protons found in the inner belt.The altitude range of confinement is compressed for more energetic xenon ions, from both ends, and tails off in the outer belt.The methodology can readily be extended to other ion propellants and initial conditions, but it is also possible to extrapolate present results.While we have not explained all of our results in detail, we believe they lay the ground work for future study.Future work could include divergence effects, filling in the parameter space, performing more comprehensive data fits, simulating propellants other than xenon, including particle collisions and ion scattering from wave-particle interactions, performing more detailed trajectory  both from an increase in LTS cases but also a small spike in BMA outcomes.These changes illustrate that magnetic confinement has been degraded due to the high energy.
To gauge the influence of inclination on ion behavior for the general cases at elevated energy, Fig. 32 breaks out these results for 39390 eV.
These high-energy results reveal stark differences in terminating conditions between inclinations.A lower-altitude shift of BMA and LTS cases is still observed for the more

Fig. 3
Fig.3Relationship between satellite inclined circular orbit and GEI coordinates at the Winter Solstice.The right ascension of ascending node, , was chosen to be 90•

Fig. 4
Fig. 4 Earth's tilt orientation and position relative to the Sun throughout the year

Fig. 5
Fig. 5 Flowchart of process for data collection over various mission designs

Fig. 6
Fig. 6 Example of the AeroTracer user layout (plotted in GSM coordinates), with labeled sections identifying key features of the program.The equatorial plane view shows that the ion's trajectory terminates upon reaching the magnetopause boundary

Fig. 7
Fig. 7 Illustration of the wide variance in simulated time for Section III. C. (Table 5) cases resulting in a "Maximum Number of Steps" terminating condition.Data are for 300eV, zero inclination and solar minimum

Fig. 9 Fig. 10
Fig. 9 Orbit altitude change maneuvers for an Earth-orbiting spacecraft, where the spacecraft's velocity vector points in the counterclockwise direction around the Earth

Fig. 11
Fig. 11 Number of terminating condition as a function of altitude across all values, where pink and purple shadings represent the inner and outer radiation belts, respectively

Fig. 14 a
Fig. 14 a Number of terminating conditions separated by inclinations of 0 • and 60 • .b LTS and BMA outcomes for same data as Fig. 14a, fit by a logistic function

Fig. 15 Fig. 16
Fig.15 Separation of inclinations at 300 eV for general cases

Fig. 20 Fig. 21
Fig.20 Results for double and single ion release (including both orbit raising and lowering cases)

Fig. 22
Fig. 22 Orbit raising results for double and single ion release

Fig. 23
Fig.23 Investigating the effect that ion charge and ion velocity direction have on an ion's terminating condition as a function of altitude

Fig. 24 Fig. 25 Fig. 26 Fig. 27
Fig.24 Results for solar maximum and solar minimum ion release (including both orbit raising and lowering cases)

Fig. 29
Fig. 29 Effects of release energy and ion velocity direction on an ion's terminating condition as a function of altitude, considering singly charged particles released at 0 • inclination during Solar Minimum on Day 81

Fig. 30
Fig. 30 Frequency of MNS terminating condition as a function of energy across all general cases, presented on a semi-logarithmic scale

Fig. 31
Fig. 31 Frequency of terminating condition as a function of altitude across all energies and values, where pink and purple shadings represent the inner and outer radiation belts, respectively

Fig. 33
Fig.33 Average of outcome frequency on altitude taken over all tested days of year.Outcomes were summed over , inclination, all energies, and direction parameter results for the general test simulations

Fig. 34
Fig.34 Deviation of outcome frequency for each day of year from the average.Outcomes were summed over , inclination, all energies, and direction parameter results for the general test simulations

Table 1
Tabulated geocentric inertial velocities for the altitudes, energies, and ion velocity directions that were simulated in the study

Table 2
Listing of parameters for a subset of test cases (referred to as general cases).Unless otherwise noted, all possible combinations of a given parameter were simulated and bold (blue) text indicates that partial completion for a given characteristic.In total, 3840 cases were tested for this study ⋆ Varied parameters include altitude, , and direction.†Varied parameters include altitude, , direction, and DOY.Varied parameters include altitude, , direction, and DOY

Table 3
Summary of parameters for general tests at the moment of ion release, where all 1536 possible combinations of a given parameter were simulated

Table 5
Summary of parameters for doubly charged ion tests at the moment of ion release, where all 768 possible combinations of a given parameter were simulated

Table 6
Summary of parameters for ion tests occurring at the extrema of the solar cycle, where all 768 possible combinations of a given parameter were simulated