Modeling air-breathing electric propulsion is complex since it involves a vast amount of physical processes, mainly relevant to the fields of rarefied gas dynamics, chemistry, and plasma physics. The modeling approach and results in terms of requirement on thruster performance and intake sensitivity we presented hereafter are applicable to air-breathing, double-staged electrostatic devices. In this propulsion concept, some of the atmospheric particles impinging on the spacecraft at orbital velocity are collected by the intake and thermalized at wall temperature. By colliding with the intake and thruster walls, the particle kinetic energy is passively converted into internal energy. The thruster’s first stage is an ionization stage offering a relatively large open inlet area and volume to collect and ionize the incoming flow. It uses electrical power available from solar arrays to ionize the compressed flow by impact with electrons emitted from a pipe-fed cathode mounted internally in the thruster assembly. Inside the ionization stage, a set of electromagnets generate a magnetic field able to significantly reduce electron mobility. A proper magnetic field topology and intensity enhance ionization and confine the ions from the chamber walls. By the combined effect of the magnetic field and voltage applied to the electrodes, the ionization stage increases flow compression by pumping and channeling the ionized particles into the second stage (the acceleration stage), providing atmospheric propellant conditions suitable for the specific accelerator technology. The proposed physical description consists of a system of 0D ODE equations. Each equation represents particle continuity of the s neutral species or s + singly ionized species in the ionization stage or acceleration stage control volumes, which are, respectively, labeled with the subscript i and a. Each equation solves for the s or s + species number density n in the ionization or acceleration stages. In the two thruster control volumes, particle continuity is coupled with electron energy conservation to compute the time-dependent electron temperatures T _e,i and T _e,a .

We make use of 3D Monte Carlo methods to describe neutral and ion dynamics and compute the set of parameters defined in Table 1, including particle transmission (also known as Clausing factors) and mean residence time in the two thruster control volumes together with mean axial velocity of the particles leaving the propulsion system domain. The neutral Monte Carlo, discussed in Section 3, is based on the Cercignani–Lord–Lampis wall collision model. The ion Monte Carlo, not relevant to the results provided in this study, consists of a 3D Boris pusher in static 2D EM fields computed by solving ampere and electron drift-diffusion PDEs in a FreeFemm++ [13] environment. More details are given as follows:

• for a given axisymmetric thruster design, we at first computed the applied magnetic field B = ∇ × A from the magnetic vector potential A = A _θ (r, z):

In line with the 0D approach, we approximated ∇ ⋅Γ _i /n ∼ κ, where κ = κ (t), a time-varying function to be matched consistently with the 0D ODE system solution in each thruster control volume. The 2D potential solution seems unimpacted by the value of κ, and a single iteration is usually sufficient to ensure consistency.

Now, consider the case in which the thruster is in the OFF state. In the ionization and acceleration stage control volumes, we express neutral s-particle continuity as follows: d n s , i d t = α s , i t k Γ s , ∞ V i − n s , i τ s , i , (9) d n s , a d t = V i V a α s , i n s , i τ s , i − n s , a τ s , a , (10) where Γ _s,∞ = n _s,∞ u _s,∞ A _itk is the particle flux incident to the intake inlet area A _itk , and V _i /V _a is the ionization to acceleration stage volume ratio. The steady state particle fluxes through the propulsion system control volumes are visualized in Figure 4A. For a propulsion system in the OFF condition, thruster transmission α _s,thr is equal to the product of the ionization stage and acceleration stage transmissions, α _s,i and α _s,a , respectively.

At steady-state, the number density in the ionization stage is increased by the compression ratio r, as follows: r = n s , i n s , ∞ = u s , ∞ A i t k V i α s , i t k τ s , i . (11)

Now, consider the case in which the thruster is switched ON, Figure 4B. In the ionization and acceleration stage control volumes, we express ion and neutral particle continuity as follows: d n s , i d t = α s , i t k Γ s , ∞ V i − n s , i τ s , i + 1 − α s + , i n s + , i τ s + , i + ω ̇ s , i , (12) d n s + , i d t = − n s + , i τ s + , i + ω ̇ s + , i , (13) where the four terms in the right-hand side of Eq. (12) represent the particle flux collected from the intake, flow convection, ion flow neutralized at walls (we assume all unchanneled ions to be neutralized at ionization stage walls), and a source/sink term ω s ̇ , respectively. The latter is related to ionization and dissociation phenomena involving energetic electrons in the thruster discharge, as follows: ω ̇ s , i = ∑ j ν s , j k j T e , i n j , i n e , i , (14) where k _j (T _e,i ) is the electron temperature-dependent reaction rate of the jth reaction; the neutral or ionized species number density n _{j, i} and electron number density n e , i = ∑ n s + , i (assuming quasi-neutrality) are the reactants of the jth reaction; and ν _s,j is the stoichiometry coefficient defining the amount of species s created or destroyed in the jth reaction. Electron impact reaction rates strongly depend on electron temperature, which, in each control volume, is computed by applying electron energy conservation (see Section 1 of Supplementary Material). Reference cross sections for electron impact ionization, dissociation, and dissociative ionization for N, N₂, O, and O₂ can be found in the studies mentioned in references [14–17]. Given the flow density in the altitude range of Figure 3 and typical discharge electron temperatures in the 10–100 eV range, we estimated chemical reactions at walls and among neutrals, ion-electron recombination, and negative ionization to be of seconday importance compared to electron-impact reactions. Nonetheless, these processes could be included as additional source/sink terms in the continuity equations. An equivalent formulation is used for the acceleration stage control volume, as follows: d n s , a d t = V i V a α s , i n s , i τ s , i − n s , a τ s , a + 1 − α s + , a n s + , a τ s + , a + ω ̇ s , a , (15) d n s + , a d t = V i V a α s + , i n s + , i τ s + , i − n s + , a τ s + , a + ω ̇ s + , a . (16)

Similarly to neutral particles, the thruster ion transmissions α _s+,i and α _s+,a and mean ion residence time τ _s+,i and τ _s+,a describe the ionized particle dynamics, and they are all figures of merit for the thruster operation. In particular, the ion transmission α _s+,i quantify the capability of the ionization stage to channel the ions toward the acceleration stage while avoiding their neutralization at walls. Particle confinement in the ionization stage, described by the mean residence time τ _s+,i , quantifies the capability of the thruster to increase plasma density and ionize more effectively.

Now, considering momentum conservation over the entire propulsion system control volume, the thrust T and the net thrust T _N are defined as follows: T = V a ∑ s + M s α s + , a n s + , a u s + , a τ s + , a + V a ∑ s M s α s , a n s , a u s , a τ s , a , (17) T N = V a ∑ s + M s α s + , a n s + , a u s + , a τ s + , a + V a ∑ s M s α s , a n s , a u s , a τ s , a − ∑ s M s u s , i t k 1 − α s , i t k Γ s , ∞ + V i 1 − α s , i n s , i τ s , i − ∑ s M s Γ s , ∞ u ∞ , (18) where the four terms in the right-hand side of Eq. (18), visualized in Figure 5, respectively, represents the momentum of the accelerated ions, the momentum of the neutrals effusing from the acceleration channel, the momentum of the particles not transmitted through the thruster and effusing back through the intake, and the inlet particle momentum. Consistent with Figure 5, it is now convenient to define three propulsion system performance indicators: 1) the collection efficiency η _c , as follows: η c = V a ∑ s + M s α s + , a n s + , a τ s + , a + ∑ s M s α s , a n s , a τ s , a / ∑ s M s n s , ∞ A i t k u ∞ , (19) 2) the effective exhaust velocity u _e , as follows: u e = V a ∑ s + M s α s + , a n s + , a u s + , a τ s + , a + ∑ s M s α s , a n s , a u s , a τ s , a / ρ ∞ u ∞ A i t k η c , (20) and 3) the effective inlet velocity u _∞,e or, equivalently, the drag coefficient C D i t k of the intake internal surface, as follows: u ∞ , e = C D i t k 2 u ∞ = u ∞ + 1 ρ ∞ u ∞ A i t k ∑ s M s u s , i t k 1 − α s , i t k Γ s , ∞ + V i 1 − α s , i n s , i τ s , i . (21)

The thrust and net thrust expressions simplify as follows: T = ρ ∞ u ∞ A i t k η c u e , (22) T N = ρ ∞ u ∞ A i t k η c u e − u ∞ , e . (23)

As η _c defines the probability of an incident particle leaving the propulsion system though the acceleration channel, η _c is proportional to the probability α _itk of an incident particle being originally accepted into the ionization stage. Thruster transmission α _thr = η _c /α _itk and effective exhaust velocity u _e (or, equivalently, specific impulse I _sp = u _e /g ₀) are key indicators for the thruster performance, and the primary function of the ionization stage is to significantly increase α _thr as the thruster is switched ON.

At the satellite level, solar array and platform structure external surface produces a drag, as follows: D p l t = C D p l t 2 ρ ∞ u ∞ 2 A f , (24) where A _f is the platform frontal area, and C D p l t is the platform drag coefficient. Note that in free molecular flows, the drag coefficient can be significantly affected by the spacecraft’s lateral surfaces geometry. Full drag compensation is achieved when the following condition is met: T N > D p l t → u e > u ∞ η c C D i t k 2 + A f A i t k C D p l t 2 . (25)

In Eq. (23), we see that the condition to provide a positive net thrust T _N > 0 is as follows: u e > U e = u ∞ , e α i t k α t h r ∼ u ∞ α i t k α t h r . (26)

Depending on the achievable α _itk and α _thr , exhaust velocities on the order of 10 km/s to 100 km/s are required, Figure 6A. This range of exhaust velocity is commonly achieved by electrostatic devices, which are now a mature technology, counting on hundreds of successful flight applications. For electrostatic accelerators, the effective exhaust velocity is usually expressed as follows: u e = η T 2 e ϕ / M ∞ , (27) where ϕ is the applied acceleration voltage and η _T is the thrust efficiency. η _T accounts for several thrust detrimental effects, such as incomplete propellant ionization, ineffective ion acceleration, unstable thruster–cathode coupling, and plume divergence [29]. Figure 6B shows the achievable exhaust velocity vs. acceleration voltage and thrust efficiency for a reference atmospheric propellant mass of 21 amu, representative of a 50/50 N₂/O composition. For thrust efficiencies higher than 50%, which is a typical value for state-of-the-art electrostatic devices, the criterion u _e /U _e > 1 can easily be met. Figure 6C shows the behavior of u _e /U _e vs. acceleration voltage and overall system performance α _itk α _thr η _T . For 25% intake transmission (representative of the intake test case selected in Section 4.1) and considering 0.25, 0.5, and 0.75 as reference values for α _thr η _T , only acceleration voltages higher than 2000, 450, and 250 V can provide a positive net thrust, Figure 6D. Even if achieving sufficiently high thruster performance with atmospheric propellant is not obvious, acceleration voltage levels in the 300–3000 V range are commonly used by electric engines operating with xenon.

Now, considering any target air-breathing EP application scenario, a mission requirement can be formulated as follows: ∑ i = 1 I O N ∫ t i , 0 t i , 1 T N , O N t d t + ∑ j = 1 I O F F ∫ t j , 0 t j , 1 T N , O F F t d t ≥ ∑ i = 1 I O N ∫ t i , 0 t i , 1 D p l t t d t + ∑ j = 1 I O F F ∫ t j , 0 t j , 1 D p l t t d t , (28) where I _ON and I _OFF are the total number of thruster ON/OFF cycles and t _i,0, t _i,1, t _j,0, and t _j,1 are the start and end time of the ith/jth ON/OFF cycle. It is reasonable to assume that Ram-EP propulsion is constrained to operate when solar power is available. In this regard, SSO dawn–dusk orbits are a preferred choice as they offer minimum eclipse duration and an almost constant solar flux direction, which allows for reduced drag by orbiting with solar arrays aligned with the incoming flow. By considering time-averaged quantities, Eq. (28) can be simplified as follows: t O N η c u e − t O N C D i t k , O N 2 u ∞ − t O F F C D i t k , O F F 2 u ∞ > t O N + t O F F A f A i t k C D p l t 2 u ∞ , (29) where t _OFF and t _ON are the total OFF and ON time accumulated throughout the mission. Based on Eqs (19) and (21), we can see that the intake drag coefficient is related to the collection efficiency η _c and mean of the axial velocity distribution of the reflected particles u _z,itk as follows: C D i t k , O N ∼ 2 + ( 2 − 2 η c ) u s , i t k / u ∞ and C D i t k , O F F ∼ 2 + 2 u s , i t k / u ∞ , where a thruster transmission α _thr ∼ 0, thus a collection efficiency η _c ∼ 0 is considered when the propulsion is in the OFF state. The requirement of minimum thruster exhaust velocity is as follows: u e > u ∞ α i t k α t h r C D i t k , O N 2 + A f A i t k C D p l t 2 + t O F F t O N C D i t k , O F F 2 + A f A i t k C D p l t 2 . (30)

Platform drag can be reduced by relaxing solar array area requirement, for example, by reducing the thruster power consumption or by increasing solar cells and power distribution efficiencies. In general, platform performance may be quantified as available power P _av per unit platform drag area C D p l t A f . Since the thruster discharge power is proportional to squared exhaust velocity P D ∝ u e 2 , we can see that relaxing the requirement of Eq. (30) is always beneficial.

We here define as “optimal” the intake designs that minimize the minimum required exhaust velocity, given the constraint on intake wake volume (Figure 2) for low-drag integration of main platform subsystems. The requirement on minimum exhaust velocity expressed in Eq. (30) can be relaxed by decreasing the ratio t _OFF /t _ON , which is driven by mission/platform requirements and is not influenced by the intake design, and by increasing the collection efficiency η _c and the intake to frontal area ratio A _itk /A _f . A proper intake and platform aerodynamic design reducing both intake and platform drag coefficients is also beneficial. In free-molecular flows, the drag coefficient C _D is affected by several variables, such as flow properties, wall material, temperature, and geometry [30], and can not be estimated without a detailed assessment of flow properties and platform geometry/thermal behavior. However, the impact on drag caused by the dependence of C _D on lateral surface shape and wall/flow temperature is usually of secondary importance when compared to the impact of the frontal area A _f . In addition, in the absence of detailed information on platform design, the impact of intake geometry on A _f can only be estimated by simple proportionality relations, implying that the intake design must be optimized iteratively throughout the whole platform design process. As such, we here consider as “weakly optimal” the intake designs maximizing the achievable η _c A _itk /A _f (which is the ratio of the propellant collecting area and the platform frontal area), given the platform constraint on intake wake volume. The resulting weakly optimal intake designs, discussed in Section 4, can then be used as an initial guess for further optimization during the advancement of air-breathing platform design and consolidation. In this study, we limit our analysis to a simplified intake geometry consisting of a truncated conical shape having an inlet radius R ₁, an outlet radius R ₂, and a length L ₁. This is the simplest shape to consider for connecting the intake to the thruster, allowing payload and other spacecraft systems housing in the free volume behind the intake wake. We described any conical intake geometry by means of two dimensionless parameters, defined as the intake aspect ratio d ₁ = L ₁/R ₁ and the intake to ionization stage inlet area ratio d 2 = ( R 1 / R 2 ) 2 . We supposed that the platform payload and main subsystem are not directly exposed to the incoming flow, and that solar arrays are the major drag source. For conical intake geometries, we defined the volume utilization ν as follows: ν = d 1 d 2 1.5 2 3 − 1 3 d 2 − 1 3 d 2 , (31) which is simply the ratio of intake wake volume and cubed thruster inlet radius R 2 3 . The solar array area S _S/A is proportional to the platform power demand, which in turn is proportional to the power P required to sustain the plasma discharge, as follows: S S / A ∝ P ∝ ρ ∞ u ∞ A i t k η c u e 2 . (32) Considering now quasi-continous thruster operation t _OFF /t _ON ∼ 0 and a platform frontal area comparable to the intake inlet area A _f ∼ A _itk , we have u _e ∝ u _∞ /η _c , and Eq. (32) can be reformulated as follows: S S / A ∝ ρ ∞ A i t k u ∞ 3 α i t k α t h r . (33) Assuming a minimum value of neutral density for which the thruster can efficiently operate, a higher intake compression allows the spacecraft to orbit at lower atmospheric densities ρ _∞ ∝ 1/r, where r is the intake compression ratio as defined in Eq. (11). Moreover, the transmission α _thr is not influenced by the intake design as it is driven by thruster performance only. We thus estimated the influence of the intake design on solar array area requirement as follows: S S / A ∝ A i t k r α i t k . (34) Considering solar arrays of length L _S/A and thickness h _S/A to be aligned with the atmospheric flow, as it would be in a SSO dawn-dusk mission scenario, we have the following: A f ∝ h S / A L S / A S S / A . (35) As h _S/A is not influenced by the intake design, and assuming the solar array length L _S/A to be comparable to intake length L i t k ∝ d 1 d 2 , we evaluated the impact of the intake design on the platform frontal area as follows: A f ∝ A i t k r α i t k d 1 d 2 , (36) and compute the to-be-maximized intake collecting area to platform frontal area as follows: η c A i t k A f ∝ α i t k 2 r d 1 d 2 , (37) where we used η _c ∝ α _itk , since the thruster transmission α _thr is not impacted by the intake design. In the following section, we presented a simple Monte Carlo algorithm aimed at computing the performance of conical intake geometries. Weakly optimal intake geometries, here defined as the ones maximizing α i t k 2 r d 1 d 2 given the constraint on volume utilization ν, are discussed in Section 4.1.

Without loss of generality, we computed intake performance for a pure N₂ flow. Performance sensitivity on flow composition is then assessed in Section 4.2. Figures 9A–C, respectively, show the impact of intake geometry on N₂ transmission α _itk , compression r, and volume utilization ν for the intake aspect ratio d ₁ and the area ratio d ₂ ranging between 2 and 20. As the intake transmission decreases monotonically with d ₂, the flow rate ∝ α _itk d ₂ accepted into the thruster increases less than linearly with the intake inlet area. For area ratios d ₂ < 3, the high degree of inlet flow collimation allows achieving a flow transmission greater than 50% also for large aspect ratios d ₁ > 10. Higher values of d ₁ increase the mean particle residence time, being more and more difficult for a thermalized particle in the ionization volume to be scattered back to the external environment. As discussed in Section 2.3, a weakly optimal intake design maximizes the achievable collection area to the frontal area ratio ∝ α i t k 2 r d 1 d 2 given the constraint on volume utilization ν. The front of maximum achievable α i t k 2 r d 1 d 2 for volume utilization constraint ranging from 50 to 450 and associated design solution are respectively visualized in Figures 9D,E as computed by means of a sorting algorithm from the dataset of Figures 9A,C. Intake performance in terms of transmission and compression associated to the optimal solutions is provided in Figure 9F. According to our definition of intake fitness, we see that the optimal intake aspect ratio d ₁ is a decreasing function with the area ratio d ₂, ranging from a value d ₁ ∼ 18 at d ₂ = 2 down to d ₁ ∼ 8 at d ₂ = 20. The resulting intake transmission and compression range from about 0.4 to 0.1 and 180 to 240, respectively. Depending on platform needs, an optimal trade-off exists between wake volume available for platform subsystems integration and performance achievable by the air-breathing propulsion system. For the purpose of the sensitivity study addressed hereafter, we considered the reference intake test case maximizing the product between ν and α i t k 2 r d 1 d 2 , that is, the d ₁ = 9 and d ₂ = 7.5 intake design depicted in Figure 7.

The effect of flow composition on intake performance is evaluated in Figure 10 for several N₂/O/O₂ mixtures, including pure N₂, pure O, and pure O₂ flows. Both transmission and compression increase by increasing the average propellant molecular mass, since a smaller molecular mass implies higher thermal velocities as well as a reduction of flow collimation and residence time. Comparing pure O with respect to pure O₂ flow, transmission decreases from about 24 to 19%, while the achievable compression decreases from about 240 down to 140. Even if O₂ is a minor constituent in VLEO, it is a species of interest for simulating very low orbits during on-ground testing of air-breathing systems, as storing oxygen in atomic form is not feasible and achieving a complete O₂ dissociation in any kind of VLEO flow simulator is complex. In the following sensitivity analysis, we used a 50/50 N₂/O composition as a reference, representative of a 200 km of altitude flight scenario during medium solar activity.

Not only flow composition but also its temperature largely varies depending on orbital parameters, time and solar activity. Figures 3E,F show a yearly flow temperature evolution at 180 and 250 km SSO dawn dusk orbits, where typical temperature values are between 700 and 1000 K. In principle, a higher inlet flow temperature has a detrimental impact on intake performance, as the increased flow thermal velocity decreases the inlet flow collimation. However, as illustrated in Figure 11, the actual impact on performance turns out to be minor and on the order of a few percentage points between the two considered temperature extremes. Admitting the possibility of a Ram-EP S/C to target not only circular but also elliptical orbits, and in general to be capable of performing different maneuvers, the impact of inlet flow velocity (in the 5 km/s to 15 km/s magnitude range) is also assessed in Figure 11. A higher inlet flow velocity increases significantly both intake transmission and compression, but, of course, it implies a higher drag and a more demanding requirement on the thruster discharge power consumption needed to provide a positive net thrust.

As discussed in Section 3, in our propulsion system the flow dynamics are dominated by particle-wall collisions, and, as indicated in Table 2, the reflected velocity distribution is affected by wall temperature. In general, intake and thruster wall temperatures may vary during the mission, depending on thruster discharge power (which, in turn, depends on the inlet flow density and exhaust velocity) and other S/C internal heat sources, together with direct Sun radiation, Earth’s albedo, IR radiation, and free-molecular heating from the atmospheric flow. Figures 12A–C show the performance sensitivity on intake and thruster wall temperature ranging between 300 and 800 K for a reference 50/50 N₂/O composition. A higher intake temperature slightly decreases the flow transmission, from about 22 to 18% at the two temperature extremes considered. Higher intake and thruster wall temperatures both reduce the mean particle residence time in the ionization stage. As a consequence, the resulting compression ratio is reduced from a maximum of 220 to a minimum of 120 among the considered wall and thruster temperature extremes. Thus, cooler walls improve intake performance.

Intake performance sensitivity on normal energy and tangential momentum accommodation coefficients α _n and α _t is visualized in Figures 12D–F. Consistently with most materials used in practical applications [31], we let α _n and α _t range from 0.1 to 1 and 0.8 to 1, respectively. A larger normal and a lower tangential accommodation improve intake performance. Intake transmission increases from about 15 to 29% as α _t is reduced from 1 to 0.8 and α _n is increased from 0.1 to 1. This behavior is consistent with the calculations for conical holes reported in [31]. The impact of wall accommodation add uncertainty in predicting the in-orbit performance of air-breathing systems, as α _n and α _t are difficult to assess with good precision and may change during the mission lifetime due to atmospheric oxygen contamination.

In previous analyses, we always assumed the intake axis to be perfectly aligned with the relative spacecraft-flow velocity vector. In general, this condition is not always met throughout the mission duration, depending on orbital dynamics, spacecraft attitude and external disturbances. Considering that any intake—flow velocity misalignment affects the output particle transmissions and distributions, a dedicated sensitivity study was performed to address this issue. Spacecraft’s rotations can directly be related to orbital velocity drifts along the S/C body-fixed reference frame. Consistently with Figure 7, we identify the roll x-axis as parallel to the S/C longitudinal direction and directed toward the intake inlet area, the S/C pitch y-axis as directed toward the zenith pointing S/C solar wing and the yaw z-axis to form a right-handed orthonormal frame. The inlet particle velocity distribution is approximated as f v , ∞ = π − 3 / 2 v t h , ∞ − 3 e x p − v x − u ∞ ⁡ cos ⁡ β ⁡ cos ⁡ γ 2 + v y − u ∞ ⁡ cos ⁡ β ⁡ sin ⁡ γ 2 + v z − u ∞ ⁡ sin ⁡ β 2 / v t h , ∞ 2 , (42) where β and γ are the S/C pitch and yaw angles with respect to the incoming flow direction or, equivalently, orbital velocity.

Impact of S/C attitude on intake performance is shown in Figure 13 for yaw and pitch angles ranging between −10 and 10°. Both transmission and compression are greatly affected by intake axis—flow direction misalignment. At 10° of yaw or pitch, there is an up to 90% performance reduction, coming down to a 20% reduction at 5° yaw/pitch. Intake sensitivity to S/C attitude can however be reduced by shortening the overall intake length, at the expanse of a reduced volume available behind the intake wake. As a comparison, at drag levels below 50 mN the GOCE S/C was capable of guaranteeing ±6° yaw/pitch by using magnetorquers only [7].