Neutral flow interaction with a magnetic dipole plasma. I. Theory and scaling

The interaction between a high-speed neutral gas flow and magnetized plasma is a problem of fundamental interest to space plasma physics,^1–4 plasma propulsion,^5,6 atmospheric reentry,^7,8 and plasma aerocapture concepts.^9,10 Theoretical understanding of this problem has primarily focused on the critical ionization velocity (CIV) phenomenon.^11–13 CIV theory hypothesizes that rapid ionization of a neutral flow will occur when the kinetic energy of the flow relative to a background plasma exceeds the ionization energy of the neutral particles. The corresponding velocity is commonly referred to as the critical ionization velocity. Evidence of this effect has been found in a number of laboratory and space experiments. Laboratory experiments have focused mainly on studying the interaction of a high-speed plasma flow with stationary neutral gas.¹⁴ Homopolar plasma devices, which possess a radial electric (E) field and an axial magnetic (B) field, found the E × B velocity of the ions to be limited to the critical ionization velocity following rapid ionization of the neutral gas.^11,15 Experiments with plasma guns observed that high-speed plasmas rapidly decelerate to the critical ionization velocity upon merging with a neutral gas.^16–19 Similar experimental support for CIV has come from studying the interaction of high-speed neutral gas flows in the Earth's ionosphere using rockets that possess shaped-charge explosives²⁰ and controlled gas releases from orbiting spacecraft.^21–23 

CIV theory generally considers relative motion between a neutral gas and isotropic plasma of infinite extent. As incoming neutral particles are ionized via electron impact ionization, the resulting fast-moving ions transfer a fraction of their kinetic energy to the electron population. This increased electron energy results in a corresponding ionization rate increase—a feedback process that rapidly ionizes the neutral flow. Theoretical research^3,12,24 has focused on understanding energy transfer between ions and electrons in the collisionless regime relevant to space and laboratory CIV experiments. Less emphasis has been placed on the influence of plasma anisotropies, finite length scales, and collisional energy transfer. These effects can be expected to impact CIV through a combination of finite Larmor radius effects, spatial variations in the ionization rate, and the exchange and diffusion of mass and energy. Furthermore, they are critical to understanding momentum transfer between the neutral flow and plasma and determining specific requirements of CIV on the electron energy confinement time.

In this paper, we address these effects in a theoretical model describing the interaction between a high-speed neutral flow and dipole plasma of finite extent. In a finite system, the magnetic field geometry is critical because the location at which a neutral is ionized determines the nature of its resulting trajectory. Therefore, we focus on understanding the spatial distribution of these trajectories to account for finite Larmor radius effects and the capture of mass and energy from trapped particles. This analysis enables derivation of general scaling laws for the transfer of mass and energy from the neutral flow to the plasma and for momentum coupling to the magnetic field. These laws are used to describe the processes governing distinct regimes of plasma–flow interaction and the critical ionization-like transitions between them. In Paper II of this series,²⁵ we derive a higher-fidelity global model that combines the scaling laws found in this current paper with physical models for diffusion, charge exchange (CX), ionization, and ion–electron energy transfer.

The comprehensive approach presented here is unique in addressing critical ionization in a finite, inhomogeneous magnetic topology. It is important to note that we primarily consider collisional interaction between species, such as charge exchange and electron-impact ionization, as a departure from many prior CIV models. These effects are most readily invoked for low-temperature (⁠$T<102$ eV), high density (⁠$n>1016$ m⁻³) plasmas. We assume a relatively high magnetic field strength (B > 100 G) so that the timescale for thermal equilibration is generally faster than that for particle diffusion from the system. Therefore, ions that become trapped in regions of high magnetic field strength provide a source of mass and energy for the plasma. Furthermore, we assume that magnetic pressure dominates plasma pressure (low-beta), so the magnetic field is static and unaffected by the plasma.

To formulate this model, we consider an axisymmetric magnetized plasma whose axis of symmetry is oriented parallel to an incident neutral flow (Fig. 1). This alignment allows the use of 2D cylindrical coordinates (r, z), where symmetry dictates that solutions are uniform in θ. The density and velocity of the neutral flow are denoted as $n∞$ and $u∞$⁠, respectively. We assume that the thermal energy of the neutral flow is negligible compared to its kinetic energy, $ε∞≈msnn∞u∞2/2$⁠, where $ε∞$ is the freestream energy density and m_sn represents the neutral particle mass.

For simplicity, we consider the field topology of a magnetic dipole formed by a closed loop of electric current. The magnetic field can be described by the flux function²⁶ 

with

Here, $ψ̂=ψ/(B0rc2)$ is the normalized flux function, B₀ is the magnetic field strength at the magnet center, and r_c is the magnet radius. K and E are complete elliptic integrals of the first and second kind, respectively. Cylindrical coordinates are normalized by the magnet radius: $r̂=r/rc$ and $ẑ=z/rc$⁠. The normalized magnetic field vector, $B̂=B/B0$⁠, can be found from

where $êθ$ is the unit vector in the θ direction.

A contour plot of $ψ̂$ is shown in Fig. 1 for normalized cylindrical coordinates. Here, lines of constant $ψ̂$ represent the intersection of magnetic flux surfaces with the $r̂$–$ẑ$ plane. We note that the magnetic topology presented in this coordinate system is independent of both r_c and B₀. This property will be leveraged in Secs. II C and II D when performing surface and volume integrals in $(r̂,ẑ)$ coordinates. Also note that this formulation assumes an infinitesimal current loop. We therefore neglect particle loss to the magnet surface, which is a fair approximation provided that diffusion scales inversely with B (or B²). This allows us to exploit the theoretical literature on plasma dipole equilibria for simplification of the pressure profile, the validity of which is discussed in Sec. II C.

New ions with charge q_i and mass m_i are generated through charge exchange and electron impact ionization as the neutral flow encounters the dipole plasma. The net behavior of the plasma and neutral flow depends on the trajectories of newly created ions in the vicinity of the magnetic dipole. Understanding the impact of non-uniform magnetic fields on the resulting trajectories will inform not only momentum transfer to the magnet, but also the definition of a control volume for energy and mass accounting. Here, we examine ion trajectories as a function of the location at which they are created to understand their behavior in the presence of the dipole magnetic field.

The equation of motion for an ion in a magnetic field without collisions can be expressed in the dimensionless form as

where

represents a characteristic ion Larmor radius normalized by the magnet radius. Here, we have introduced the dimensionless variables: velocity vector $û=u/u∞$⁠, position vector $r̂=r/rc$⁠, and time $t̂=tu∞/rc$⁠. From Eqs. (4) and (5), the ion trajectory as a function of time, $r̂(t̂)$⁠, is uniquely determined by ρ_L and the initial ion location and velocity, $r̂(0)$ and $û(0)$⁠, respectively. Charge exchange (CX) and electron-impact ionization are the two mechanisms by which streaming neutrals become ions. In both cases, the newly formed ion will initially maintain the neutral flow velocity in the axial direction, $û(0)=êz$⁠, where $êz$ is the unit vector in the axial direction. Note that while the ion is created by a collision, we consider its resulting trajectory neglecting any further collisions. This is a fair approximation in the limit of an ion Hall parameter $Ωi>1$ applied only for the purpose of determining whether a given ion is trapped by the field.

We examine the influence of the non-uniform magnetic field by considering the trajectory of ions formed at various locations within the $r̂$–$ẑ$ plane. Equations (4) and (5) are propagated many times longer than the transit time of an ion across the magnet (⁠$t̂max=200$⁠) using the LSODA differential equation solver.²⁷ A few sample trajectories are depicted in Fig. 2. We observe that ions originating upstream of the magnet are either deflected or reflected by an amount that depends largely on their initial radial location, $r̂(0)$⁠. Ions formed away from the magnet (i.e., $r̂(0)>1$⁠) deflect by an amount that increases as $r̂(0)$ decreases [Figs. 2(a)–2(c)]. Ions formed on a collision course with the magnet [i.e., $r̂(0)∼1$] tend to reflect because they experience stronger magnetic fields throughout their trajectory [Fig. 2(d)]. The behavior of ions created near the axis of the magnet (⁠$ẑ∼0$⁠) is distinct from those created prior to reaching the magnet. Notably, a region of space exists where ions enter closed (or trapped) trajectories about the magnet [Figs. 2(e) and 2(f)]. A second trapped-orbit region exists at slightly larger initial radii whereby the ion trajectories exhibit orbits similar to the well-known banana orbit in tokamak plasmas²⁸ [Fig. 2(g)]. Finally, ions formed beyond this second trapped region exhibit deflected orbits [Fig. 2(h)].

The nature of the different ion orbits is best demonstrated using a heat map of the final axial velocity of each ion, $ûz,f$⁠, as a function of the location at which the ion was formed. Here, we consider ions formed on a uniform grid with spacing $Δr̂=Δẑ=0.1$ and calculate $ûz,f$ as the average of $ûz(t̂)$ for $t̂∈[100,200]$⁠. This averaging reduces the noise associated with trapped orbits for which $ûz$ is oscillatory with time. Note that $ûz,f=1$ for undeflected ions, $ûz,f=0$ for trapped ions, and $ûz,f=−1$ for fully reflected ions. The heat map shown in Fig. 2(i) clearly shows the spatial distribution of reflected, deflected, and trapped ions. Notably, a region of pass-through orbits exists on-axis because ions formed here do not encounter significant radial magnetic fields. We will ultimately use maps of $ûz,f$ to calculate the force of the reacting neutral flow on the magnet by taking the summation of the axial momentum change per ion multiplied by the rate at which ions are formed.

Ions with trapped orbits are particularly important because they generally reside near the magnet for timescales long compared with the collisional timescale, allowing them to reach thermal equilibrium with the plasma. As a consequence, trapped ions are a source of plasma mass and energy. Using distributions of $ûz,f$ for different ρ_L, we observe a characteristic boundary inside which all ions are trapped. This boundary is defined by a specific magnetic flux surface,

The resulting magnetic flux contour, depicted in Fig. 2, marks the transition between trapped ions and reflected/deflected ions. In other words, any ion formed from the flow in a location at which $ψ̂>ψ̂*$ will become trapped by the magnetic field. In 3D space, the $ψ̂*$-contour forms a toroidal control volume that we will use to examine mass and energy transfer between the neutral flow and plasma. We conservatively neglect the trapped-ion volume formed by the banana orbit region [Fig. 2(g)] due to the complexity associated with modeling this irregular shape. This region may be addressed by future higher-order modeling.

The results of the single particle model provide us with a crucial link between the high-speed neutral gas flow and confined plasma. As neutral particles interact with plasma ions and electrons, new ions can be formed through ionization or charge exchange collisions. Ions formed within the boundary defined by Eq. (7) provide a source of mass and energy to the plasma. Here, we use continuity to derive an equation for the neutral particle density spatial distribution as a function of key dimensionless parameters. The neutral density distribution is then used along with our expression for the trapped-particle boundary to determine the rate of mass and energy transfer from the neutral flow to the plasma.

We assume the ion and electron densities are equal everywhere (n_i = n_e) to preserve charge neutrality. In the presence of charge exchange reactions and electron impact ionization, the steady-state continuity equation then takes the following form:

Here, $Rtot=Rion+Rcx$ is the sum of the ionization (⁠$Rion$⁠) and charge exchange (⁠$Rcx$⁠) reaction rates. Note that we have introduced the subscript sn to distinguish between incident neutral particles (viz. stream neutrals) and neutral particles that are formed via charge exchange reactions (viz. secondary neutrals). Because the cross section for collisions between stream neutrals and secondary neutrals is much smaller than ionization and charge exchange cross sections for the stream energies of interest (>1 eV), the secondary neutrals do not play a significant role in mass and energy transfer from the neutral flow to the plasma. However, as we will discover in Paper II,²⁵ secondary neutrals play a key role in the global balance of mass and energy within the plasma.

Two simplifying assumptions enable us to reduce Eq. (8) into a form amenable to analytical modeling. First, we assume that $usn=u∞êz$ everywhere, which is equivalent to saying that the stream neutral flow is directed along $êz$ with a Mach number $Msn≫1$⁠. Second, following Kesner et al.²⁹ we assume the dipole is driven to a stationary state and any fluctuations are damped out on timescales faster than the characteristic evolution time of the system. This second assumption allows us to approximate the electron density as $ne/ne,r=(ψ̂/ψ̂r)α$⁠. Here, $ne,r$ is the electron density along an arbitrary reference surface of constant scalar flux $ψ̂r$⁠. The parameter α is a constant that represents the steepness of the density profile. In our subsequent analysis, we take α = 4 to be consistent with theoretical predictions for plasma confined by a strong magnetic dipole.^29,30 Note that the density profile considered here neglects the effects of anisotropy and the interaction between the plasma and magnet surfaces.^30,31 Diffusion across the $ψ̂*$ surface is expected to far exceed any possible cross field loss. Despite their limitations, these simplifications enable us to examine the general scaling of mass and energy transfer over a wide parameter space. The inclusion of higher order effects is required for more accurate predictions within a narrower parameter space; however, this is outside the scope of the present analysis and is left for future work.

With the above simplifications, the solution to Eq. (8) is

where $n̂sn=nsn/n∞$ and

Here, I_n is a function of $r̂$ and $ẑ$ that depends only on the specific magnetic field topology. According to Eq. (9), the stream neutral density distribution forms a wake whose size depends on the dimensionless scalar quantity

We note that $ζtot$ scales with the ratio of the characteristic ion formation rate and the characteristic rate at which stream neutrals transit the magnet.

The stream neutral density distribution is shown for four different values of $ζtot$ in Fig. 3. It is clear that the wake formed by the interaction of the neutral flow with the dipole plasma increases in size as $ζtot$ increases. This result is due to an increased probability that a stream ion undergoes an ionization or charge exchange reaction with the plasma prior to transiting past the magnet. Also shown in Fig. 3 are magnetic flux contours that define the trapped-ion volume for three different values of ρ_L, as determined from Eq. (7). Notably, for each ρ_L, there exists a value of $ζtot$ above which the trapped-ion volume is shadowed by the wake (i.e., $n̂sn≈0$ for $ψ>ψ*$⁠). This shadowing effect will ultimately play an important role in the transfer of stream neutral flow mass and energy to the plasma.

The total capture rate of stream neutral particles within the trapped-ion control volume (⁠$Ṅcap$⁠) can be determined by integrating the volumetric stream neutral reaction rate over the volume $ψ>ψ*$⁠. This can be written symbolically as

Recognizing that each stream neutral brings with it $msnu∞2/2$ worth of kinetic energy, the capture rate of stream neutral energy or captured power (⁠$Pcap$⁠), is given by

It is important to note that the mass captured by the plasma is divided among the ion and secondary neutral populations by an amount that depends on the ratio of $Rion/Rcx$⁠. This division will be examined in more depth in Sec. III.

The equations for particle and energy capture can be simplified considerably by introducing the following normalization: $Ṅ̂≡Ṅ/(n∞u∞rc2)$ and $P̂≡2P/(msnn∞u∞3rc2)$⁠. Here, $Ṅ$ and P are normalized by the flux of neutral particles and kinetic energy on an area equal to $rc2$⁠, respectively. The dimensionless form of Eqs. (12) and (13) then becomes

where

For a given magnetic field topology, the integral I_sn is a function of only two dimensionless parameters: the bounds of the integral are governed by ρ_L via Eq. (7) and the spatial distribution of $n̂sn$ is governed by $ζtot$ via Eq. (9). This result suggests that variations in the mass and energy capture rates with the properties of the neutral flow (e.g., density, particle mass, and velocity) or magnetized plasma (e.g., magnet strength, size, plasma density, and temperature) are similar with respect to ρ_L and $ζtot$⁠.

The dimensionless particle (and power) capture rate is shown as a function of ρ_L and $ζtot$ in Fig. 4(a). For fixed ρ_L, $Ṅ̂cap$ increases linearly with $ζtot$ at small values of $ζtot$⁠. This linear regime results from a weak $ζtot$-dependence within the integral I_sn, which occurs when $n̂sn≈1$ for $ψ̂>ψ̂*$⁠. Physically speaking, the stream neutral wake does not shadow the ion-trapping volume at a low $ζtot$ (un-shadowed $ψ*$⁠). As $ζtot$ increases, a maximum capture rate eventually occurs. Just beyond this maximum value, further increases in $ζtot$ yield slight reductions in $Ṅ̂cap$ (partially shadowed $ψ*$⁠). Eventually, a regime is reached wherein $Ṅ̂cap$ decreases rapidly with increasing $ζtot$ (fully shadowed $ψ*$⁠).

The solid lines in Fig. 4(a) represent a simple analytical approximation for $Ṅ̂cap$⁠, described by Eqs. (A1)–(A6) in the  Appendix. Equations (A1)–(A6) are important because they enable the mass and energy capture rates to be quickly calculated without needing to solve the full integral in Eq. (15). The value of $ζtot$ that corresponds to the maximum capture rate can be estimated from this approximate form of $Ṅ̂cap$⁠, which to first order is $ζtot(m)≈0.3/ρL0.6$ for $ρL≪1$⁠. We find that full shadowing of $ψ*$ occurs for $ζtot>10ζtot(m)$⁠. Beyond this limit, the vast majority of neutral flow particles ionize and deflect prior to reaching the ion-trapping volume. In Paper II,²⁵ we incorporate the analytical approximation for $Ṅ̂cap$ within a self-consistent global model to examine the time-dependent neutral flow–plasma interaction.

A force is produced on the dipole magnet when ionized stream neutral particles either deflect or are captured by the magnetic field. In this case, momentum is transferred to the magnet via the interaction between the magnet current and the current density that results from summing over individual ion orbits. Understanding the scaling of this force is of fundamental interest to space plasma physics.^32–34 For plasma aerocapture concepts,^9,10 this force represents the drag generated when the magnetized plasma encounters a planetary atmosphere at high velocities—a quantity critical to assessing the feasibility of this technology. Here, we use the particle trajectory maps (Fig. 2) and stream neutral density distribution (Fig. 3) to derive a scaling law for the drag force generated by a high-speed neutral flow interacting with a dipole plasma.

We define the drag force, F_D, as the axial force imparted by the neutral flow to the magnet. The force density can be written as a product of the volumetric reaction rate of stream neutrals (⁠$Rtotnsnne$⁠) and the change in momentum of the newly formed ion (⁠$msnΔuz$⁠), where $Δuz=uz,f−u∞$ is a function of r and z that represents the total change in axial velocity of an ion formed at location (r, z). The volume integral of the force density over all space yields the drag force

We define the normalized drag force as

Here, $F̂D$ effectively represents the drag of the system relative to the drag on a solid disk of radius r_c, and can be thought of as a characteristic drag coefficient that includes the effects of magnetic deflection. With this definition, Eq. (16) can be written as

where

We thus arrive at a second integral function that depends on only two dimensionless parameters for a given field topology: $n̂sn$ is governed by $ζtot$ via Eq. (9) and $Δûz$ depends on ρ_L as shown in Fig. 2.

The force density distribution is shown in Fig. 5 for a wide range of ρ_L and $ζtot$⁠. Here, the force density is defined as $f̂D=(2/π)ζtotn̂snΔûz(ψ̂/ψ̂r)α$⁠. Values of $f̂D$ are calculated on an adaptive grid that focuses gridpoints in regions of high-reactivity. This technique helps to increase both the resolution of the resulting spatial distribution and numerical precision of the integrated drag force. The general features of the force density distribution can be traced back to the three variables upon which it relies: the dark region near and behind the magnet results from wake formation in the function $n̂sn$⁠; the streaks in front of the magnet, the dark band on the magnet axis, and the lobes on the magnet periphery are all features associated with the ion orbit mapping (Fig. 2) contained in the function $Δûz$⁠; and the region over which the previous two effects are important is localized near a particular flux surface via the plasma density distribution, $n̂e∼ψ̂α$⁠. We observe that the radial extent of $f̂D$ generally increases with decreasing ρ_L. For a given ρ_L, the radial extent of $f̂D$ increases with $ζtot$ before eventually reaching a saturation point at which it ceases to expand. Finally, we note that the value of $ζtot$ corresponding to the saturation point increases as ρ_L decreases.

Numerical calculations for $F̂D$ are shown as a function of $ζtot$ and ρ_L in Fig. 4(b). As the $f̂D$-distributions in Fig. 5 suggest, $F̂D$ generally increases as ρ_L decreases. Furthermore, $F̂D$ increases with $ζtot$ for fixed ρ_L up to a saturation point, with saturation occurring at higher values of $ζtot$ as ρ_L decreases. Even with the saturation effect, drag forces up to three orders of magnitude greater than the purely aerodynamic drag are observed. The solid lines in Fig. 4(b) represent a simple analytical approximation to $F̂D$⁠, described by Eqs. (A7)–(A10) in the  Appendix. The scaling of the saturated drag value is revealed by taking the limit $ζtot→∞$ of the approximate form of $F̂D$ [Eq. (A7)], yielding $F̂D(s)≈1.2/ρL$ to first order. Interestingly, this is roughly 70% of the drag that would be generated if the neutral flow were to be replaced by a fully ionized flow with equivalent velocity and density. We find that the reactivity at which saturation occurs is roughly 1–5 times higher than $ζtot(m)$⁠. Therefore, it is possible to approach drag saturation while maintaining significant mass and energy transfer to the plasma. This point has important consequences for plasma aerocapture applications because it suggests the possibility of optimizing in situ flow utilization without sacrificing drag performance. Such applications are discussed in depth in Paper II.²⁵ 

Our analysis thus far has taken the plasma density (⁠$ne,r$⁠) and temperature (through $Rion$⁠) as independent variables. However, these parameters will inevitably evolve as mass and energy is transferred from the neutral flow to the plasma. Before advancing toward greater complexity with a time-resolved model, we seek physical insight into the steady-state behavior of the system by combining the scaling laws derived in Sec. II with a simple model for plasma mass and energy balance.

We consider the conservation of mass and energy inside a control volume defined by the ion-confining magnetic flux surface, $ψ*$ (Fig. 1). Upon reacting with the plasma, flowing (stream) neutral particles provide a source of mass and energy to the control volume, as described by Eq. (14). This mass and energy is distributed among ions, electrons, and post charge-exchange (secondary) neutrals prior to leaving the control volume. Mass conservation for these species can be written in the dimensionless form as

where source terms are on the left side of each equation and loss terms are on the right side. Similarly, energy conservation for ions, electrons, and secondary neutrals can be expressed as

Here, $N̂j$ and $T̂j$ are the total number of particles and temperature of species j, respectively; $τ̂j$ represents the timescale for particle diffusion from the control volume; $τ̂ie$ represents the thermal equilibration time between ions and electrons; $Ṅ̂2n→i$ is the ionization rate of secondary neutrals averaged over the control volume; and $ε̂ion$ represents the effective ionization energy, which also accounts for additional losses due to electron scattering and electron impact excitation. The fraction of stream neutrals captured via ionizing and charge exchange collisions is given by $ξion=Rion/Rtot$ and $ξcx=Rcx/Rtot$⁠, respectively. Finally, we consider the injection of external power into the plasma, where ξ_i and ξ_e represent the fraction of injected power (⁠$P̂inj$⁠) deposited into the ions and electrons, respectively. We note that Eqs. (20)–(25) utilize the following normalizations: $Nĵ=Nj/(n∞rc3), T̂j=3Tj/(msnu∞2), τ̂j=τj/(rc/u∞)$⁠, and $ε̂ion=2εion/(msnu∞2)$⁠.

Equations (20)–(25) describe the volume-averaged, steady-state balance of mass and energy for a high-speed neutral flow interacting with a plasma. Before proceeding, we take a moment to briefly summarize the physics associated with each term within these equations. In Eqs. (20) and (21), the first and second source terms account for ion/electron creation from electron impact ionization of stream and secondary neutrals, respectively. The loss term in the ion and electron continuity equations represents diffusion from the control volume. Charge exchange collisions between stream neutrals and ions provides the only source of secondary neutrals, as described in Eq. (22). The two secondary neutral loss terms describe losses due to diffusion and electron impact ionization. In Eq. (23), the three ion energy source terms correspond to power injected into the ion population (e.g., via resonant wave heating), the kinetic power captured from the neutral flow as a result of both ionizing and charge exchange collisions, and thermal power addition from ionized secondary neutrals. The three ion energy loss terms represent the rate of ion energy lost to secondary neutrals during charge exchange collisions, ion thermal energy diffusion from the control volume, and the power transferred from the ions to electrons (e.g., via Coulomb collisions). The source terms in the electron energy equation, Eq. (24), represent power injected into the electron population and power transferred from the ions. Electron energy losses include diffusion of thermal energy and the rate of energy lost to inelastic collisions, including ionization of both secondary and stream neutrals. Finally, Eq. (25) shows that charge exchange collisions between the ions and stream neutrals provide a source of power for the secondary neutral population. This power is balanced by losses due to diffusion and electron impact ionization.

Finally, we note that the steady-state mass and energy balance model described above makes two key assumptions: (1) the thermal energy is much greater than the kinetic energy for ions, electrons, and secondary neutrals and (2) the temperature of each species is uniform. Two useful simplifications result from these assumptions. First, substitution of the plasma quasineutrality requirement (⁠$N̂e=N̂i$⁠) into Eqs. (20) and (21) yields $τ̂i=τ̂e$⁠. The diffusive loss rate of ions must therefore equal that of electrons to maintain zero net charge, which is equivalent to saying that charged particle transport is governed by ambipolar diffusion. Second, Eqs. (22) and (25) combine to give $T̂2n=T̂i$⁠. In other words, secondary neutrals remain in thermal equilibrium with the ions.

The overall power balance within the control volume can be obtained by summing over the power balance equations for each species, Eqs. (23)–(25). Noting that $ξion+ξcx=1$ and $ξe+ξi=1$⁠, the total steady-state power balance is described by

Here, the sources of power are the kinetic energy capture rate from the neutral flow and injected power. The first term on the right hand side represents the loss rate of electron thermal and frozen flow energy associated with charged particle diffusion. The second loss term is equivalent to the sum of thermal energy lost via secondary neutral and ion diffusion. In this analysis, we will refer to the left hand and right hand side of Eq. (26) as $P̂in$ and $P̂out$⁠, respectively.

An important transition can occur if the maximum rate of energy capture is much greater than the injected power or $P̂cap(ρL,ζtot(m))≫P̂inj$⁠. According to Fig. 4(a), captured power scales linearly with plasma density, $P̂cap≈N̂i/τ̂cap$⁠, for $ζtot≪ζtot(m)$⁠. Here, the characteristic particle capture time

is a dimensionless quantity that compares the particle capture timescale to the neutral particle transit time. If the dominant contribution to $τ̂e$ is independent of density (e.g., Bohm diffusion), the first power loss term in Eq. (26) also scales linearly with $N̂i$⁠. Equation (26) therefore suggests that a critical value of $τ̂e$ exists, defined by

The importance of $τ̂e*$ to steady-state power balance is demonstrated in Fig. 6, which shows the curves of $P̂out$ and $P̂in$ as a function of $N̂i$ for different values of $τ̂e$⁠. Here, steady-state equilibrium occurs when $P̂out=P̂in$⁠, satisfying Eq. (26). For $τ̂e<τ̂e*$⁠, the power lost exceeds the power captured irrespective of $N̂i$⁠. As a consequence, $P̂inj$ dominates the equilibrium power balance. For $τ̂e≈τ̂e*$⁠, the linear portion of the $P̂out$-curve approximately equals that of the $P̂in$-curve, indicating that small changes in $τ̂e$ yield large changes in the equilibrium state. Finally, for $τ̂e>τ̂e*$⁠, the equilibrium density increases to the point where the non-linear portion of $P̂cap$ balances $P̂out$⁠. In other words, increases in the charged particle diffusion time require large corresponding increases in plasma density to maintain power balance via the wake shadowing effect described in Sec. II C. The cumulation of these effects is a mode transition from a regime where power balance is dominated by injected power to a regime where it is dominated by power captured from the neutral flow.

Similar to the overall power balance in Eq. (26), a requirement on $τ̂e*$ also exists to satisfy the overall mass balance. The addition of Eqs. (20) and (22) yields

This equation is similar in form to Eq. (26) with the exception that there is an additional loss term on the right hand side due to secondary neutral diffusion. Conservation of mass therefore requires

This equation states that, for $N̂2n/τ̂2n>0$⁠, the timescale for charged particle diffusion must exceed the characteristic particle capture timescale if the plasma is to be sustained solely by mass and energy from the neutral flow.

The mode transition described above has significant consequences on the fundamental interaction between a high-speed neutral flow and dipole-confined plasma. Taking the case shown in Fig. 6 as an example, a ninefold increase in $τ̂e$ yields a nearly four order of magnitude increase in both $N̂i$ and $P̂cap$⁠. We take a moment to examine the physical meaning of $τ̂e*$ to determine how the mode transition depends on energy transfer processes between the flow and plasma. Rearranging Eq. (28) yields $(1−T̂i)/τ̂cap=(T̂e+ε̂ion)/τ̂e*$⁠. The left hand side of this equation represents the net power deposited into the ions assuming neutral capture occurs only through charge exchange collisions (⁠$ξcx=1$⁠). This quantity can be thought of as the minimum amount of power the ions can absorb from flow interaction. The right hand side of the equation represents the net power lost from the electron population due to diffusion and inelastic collisions. If $τ̂e<τ̂e*$⁠, the electron energy loss per captured particle outweighs the ion energy gain, thus the plasma cannot be sustained without an additional power source. $τ̂e*$ therefore represents the minimum diffusion time required to balance power captured into the ion population with power lost from the electron population. This highlights an important consideration for the occurrence of critical ionization. Because flow energy deposits directly into the ion population, an energy channel must exist between the ions and electrons (e.g., Coulomb collisions) to sustain the electron temperature and balance frozen-flow losses. Thus, $τ̂e*$ also places bounds on the efficiency of this energy channel, as ion–electron energy transfer must be generally faster than the critical diffusion time.

The solution to Eq. (26) can be used to examine how the total input power and drag force on the dipole scale with ρ_L and $τ̂e/τ̂e*$ (Fig. 7). As anticipated, we observe the plasma to transition from a self-sustained (via $P̂inj$⁠) to a neutral-flow-sustained regime at $τ̂e/τ̂e*=1$⁠. Following the mode transition, the input power into the dipole plasma reaches a maximum that increases with ρ_L before eventually decreasing with $τ̂e/τ̂e*$ due to dipole shadowing. A similar trend exists for the drag force with the exception that $F̂D$ increases monotonically with $τ̂e/τ̂e*$ beyond the mode transition. This difference is due to the fact that the integral for $P̂cap$ is over $ψ>ψ*$⁠, while the integral for $F̂D$ is over all of space. Notably for $ρL=10−3$⁠, the plasma is observed to capture nearly four orders of magnitude more power than the injected power and upwards of three orders of magnitude more drag compared to the aerodynamic drag on a solid object with size equivalent to the magnet.

The preceding analysis considered $τ̂e$ as an independent variable to simplify the steady-state analysis and obtain intuitive physical understanding of the observed mode transition. In a real system, $τ̂e$ will depend on the properties of both the plasma and magnet. For example, increasing the magnetic field strength for a given flow will increase $τ̂e$ and decrease ρ_L. A critical magnetic field strength will therefore exist that corresponds to $τ̂e/τ̂e*=1$⁠, after which the captured power and drag will continue to increase as the systems transitions between curves of constant ρ_L.

A second mode transition occurs when $Ṅ̂2n→i≫N̂2n/τ̂2n$⁠. In this limit, nearly all of the secondary neutral particles are ionized prior to diffusing from the control volume, resulting in $ξcxṄ̂cap≈Ṅ̂2n→i$⁠. Insertion of this equation into the ion mass balance [Eq. (20)] yields

Simply stating, every neutral flow particle captured by the plasma ultimately diffuses from the control volume as an ion. We note that Eq. (31) is similar in nature to Eq. (26). Following the same logic used for the derivation of Eq. (28), a critical value of $τ̂e$ exists for this mode transition defined by

Physically, for $τ̂e<τ̂e**$⁠, the loss rate of electrons and ions via diffusion is larger than the capture rate from the neutral flow. When this occurs, the plasma is unsustainable without an additional source of charged particles. The requirement that $τ̂e>τ̂e**$ for plasma sustainment is equivalent to the Townsend criterion.³⁴ 

Elimination of secondary neutral diffusion has a significant impact on energy balance within the control volume. Substitution of Eq. (31) into Eq. (26) yields the following total energy balance equation:

Here, the sum of captured kinetic power and injected power balances the diffusion of thermal and frozen-flow power from the control volume. Considering again the limit $Ṅ̂cap(ρL,ζtot(m))≫P̂inj$⁠, energy balance can be written as

which states that the kinetic energy of each captured particle must equal the sum of the electron and ion thermal energies plus the effective ionization energy. To exist in this regime, the neutral flow must possess enough kinetic power to ionize and heat itself to the plasma temperature. The relationship between the second mode transition and Alfvén's critical ionization velocity phenomenon will be discussed in Sec. III E.

The fact that kinetic power from the flow deposits directly into the ion population places an additional requirement on $τ̂e$ with respect to the timescale for ion–electron energy transfer, $τ̂ie$⁠. The electron energy equation [Eq. (24)] can be rewritten as

where the simplification on the right hand side results from inserting the ion mass conservation equation with $N̂e=N̂i$ and $τ̂e=τ̂i$⁠. This equation states that the power lost from the electron population must be balanced by a combination of injected power and power transferred from the ions. However, the mode transitions are characterized by a transformation of the plasma into a state in which $P̂inj≪N̂i(T̂e+ε̂ion)/τ̂e$⁠. The electron energy equation therefore yields a critical value of $τ̂ie$⁠,

For $τ̂ie>τ̂ie*$⁠, the plasma cannot be sustained solely on power deposited by the neutral flow. As $τ̂ie$ decreases eventually $τ̂ie=τ̂ie*$⁠, at which point power transfer from the ions exactly balances power lost to diffusion and inelastic collisions. Once this mode transition occurs, Eq. (36) suggests that the temperature of the plasma species will adjust to maintain $τ̂ie=τ̂ie*$⁠.

The previous analysis suggests that the interaction between a high speed neutral flow and magnetized plasma exhibits three distinct physical regimes. In the first regime, the plasma is unsustainable without a source of mass and energy that is independent of the neutral flow. Due to its dependence on injected mass/power, we refer to this as the injection (I) regime. The first mode transition (Sec. III A) describes the transition of the plasma to a state where it is sustained entirely by the neutral flow. In this regime, secondary neutrals resulting from charge exchange are the dominant source for ionization. For this reason, we refer to it as the charge exchange (CX) regime. The CX regime is possible because diffusion of charge exchange neutrals enables stream neutral particle kinetic energy to be absorbed by the plasma without incurring the additional energy cost of ionization. As the stream neutral kinetic energy increases, eventually a point is reached where the kinetic energy per captured particle is large enough to overcome the ionization and thermal energy loss per diffused particle. The second mode transition (Sec. III B) occurs at this point, which is marked by an abrupt increase in stream neutral ionization. We refer to this regime as the critical ionization velocity (CIV) regime due to the well-known, eponymous velocity associated with the kinetic energy requirement of the mode transition.

The particular regime in which a flow/plasma system will exist depends on the timescales for electron confinement, stream neutral particle capture, and ion–electron energy transfer. Equations (28), (30), (32), (34), and (36) provide equations that describe the timescale requirements of a particular regime as a function of species temperatures and the ratio of the effective ionization energy to the kinetic energy of flowing neutrals. Determination of the steady-state species temperatures for a given system requires the full solution to Eqs. (20)–(25) with appropriate physical models for the various diffusion, ion–electron energy transfer, and collisional reaction processes. However, it is possible to generate an approximate mapping of the different regimes by recognizing that the ionization energy represents a minimum energy requirement for plasma sustainment.

First, we consider the limit where ion–electron energy transfer occurs much faster than particle capture, or $τ̂ie≪τ̂cap$⁠. The minimum requirement on $τ̂e$ for the CX regime can be found by letting both $T̂i→0$ and $Tê→0$ in Eqs. (28) and (30), which gives

Taking the same limit for Eqs. (32) and (34) gives the following requirements on $ε̂ion$ and $τ̂e$ for the CIV regime:

Equations (37)–(39) are plotted in Fig. 8(a), which provides a map of the different regime requirements as a function of the quantities $τ̂e/τ̂cap$ and $ε̂ion$⁠. For $τ̂e/τ̂cap>1$⁠, Fig. 8(a) indicates that decreasing $ε̂ion$ (e.g., by increasing $u∞$⁠) for a system in the I regime can induce the I–CX mode transition. Further decreases in $ε̂ion$ eventually lead to the CX–CIV mode transition. For $τ̂e/τ̂cap<1$⁠, the system exists in the I regime irrespective of $ε̂ion$⁠. At a fixed value of $ε̂ion$⁠, increasing $τ̂e/τ̂cap$ (e.g., by increasing B₀) leads to either an I–CIV mode transition or an I–CX mode transition. Notably, increasing $τ̂e/τ̂cap$ cannot induce a CX–CIV mode transition unless it is accompanied by a corresponding decrease in $ε̂ion$⁠.

The regime boundaries are modified when the ion–electron energy transfer timescale becomes on the order of or greater than the neutral particle capture timescale. From Eq. (26), it is obvious that $T̂i<1$ for the CX regime. Equation (36) therefore provides the following requirement for the CX regime:

Similarly, power balance for the CIV regime [Eq. (34)] requires $T̂i−Tê=1−2T̂e−ε̂ion<1−ε̂ion$⁠. According to Eq. (36), the CIV regime therefore requires

Figure 8(b) illustrates how the electron energy requirements influence the regime boundaries for a fixed ratio of $τ̂ie/τ̂cap$ greater than unity. The main consequence is a shift in both the CX and CIV regime boundaries toward smaller $ε̂ion$⁠. This shift occurs because a smaller percentage of the ion energy is transferred to the electron population. As a consequence, more ion energy is required per captured neutral particle in order to sustain the plasma. At a large $τ̂e/τ̂cap$⁠, the CIV boundary approaches that of Fig. 8(a) or $ε̂ion=1$⁠. In other words, the CIV boundary does not depend on the ion–electron energy transfer timescale so long as the energy is transferred before it escapes the volume via diffusion.

The boundaries presented in Fig. 8 represent the minimum requirements on $τ̂e$ as a function of $τ̂cap, τ̂ie$⁠, and $ε̂ion$⁠. Using appropriate physical models for the electron diffusion coefficient, charge exchange cross section, ionization reaction rates, and effective ionization energy, Eqs. (37)–(41) provide an easy method to determine if a system can physically exist in a given regime. Identification of the precise mode transition boundary requires knowledge of $T̂i$ and $Tê$⁠, necessitating a full solution to the mass and energy balance equations, as presented in Paper II.²⁵ 

We end our analysis with a brief discussion of the previous results in the context of Alfvén's theory for the critical ionization velocity. CIV theory predicts that rapid ionization of the neutral gas will occur when the flow kinetic energy exceeds the ionization energy of the neutral gas. The velocity corresponding to this threshold or critical ionization velocity, is traditionally defined as¹ 

where U_iz is the ionization energy of a neutral particle with mass m_sn.

The actual energy cost of creating an ion is greater than the ionization energy due to the cumulative effects of electron scattering and electron impact excitation of neutral particles. Including these effects into an effective ionization energy, $εion$⁠, the following parameter can be defined:

Here, C_i is a function of T_e due to the dependence of electron impact ionization and excitation cross sections on electron energy. In general, $Ci>1$ and decreases with increasing T_e. The dimensionless effective ionization energy used to determine the flow regime boundaries in Sec. III D can therefore be written as

For a given neutral flow species, it is possible to decrease $ε̂ion$ by increasing either T_e or $u∞$⁠.

Referencing Fig. 8(a), for $τ̂e≥τ̂cap$ and $τ̂ie≪τ̂cap$⁠, an I–CX mode transition will occur when $u∞$ exceeds a threshold velocity. Inserting Eq. (44) into Eq. (37), we find that the threshold velocity for the CX regime is bounded by

Similarly, using Eq. (38) we find that the threshold velocity for the CX–CIV mode transition must satisfy

From Eq. (46), it is apparent that the CIV mode transition occurs when $u∞>ucr$⁠. The fact that excitation losses increase the threshold velocity for critical ionization beyond u_cr was examined theoretically by Möbius et al.¹² and is supported by numerous CIV laboratory experiments.¹⁴ Equation (45) predicts the existence of another threshold velocity that is less than the CIV threshold velocity and can also be less than u_cr if $τ̂e/τ̂cap>Ci$⁠. We will refer to this second threshold velocity as the critical charge exchange velocity (CXV) due to its forcing of the I–CX mode transition. The appearance of a second threshold velocity resulting from charge-exchange collisions is distinct from the theoretical results of Möbius et al.¹² and McNeil et al.,²⁴ who conclude that CX collisions tend to decrease the CIV threshold velocity. This distinction is due to the fact that their analysis neglects: (1) the loss of ion energy during a CX collision [first term on the RHS of Eq. (23)] and (2) the recycling of charge-exchange neutrals into plasma ions (i.e., terms related to $Ṅ̂2n→i$⁠). Notably, we have found no experimental evidence of a distinct CXV threshold in the open literature. One possible reason for this is that traditional CIV experiments are performed using open magnetic field configurations in which $τ̂e≈τ̂cap$⁠. According to Eqs. (45) and (46), the two threshold velocities become indistinguishable in this limit.

For $τ̂ie≥τ̂cap$⁠, the shift of the mode transition boundary toward lower $ε̂ion$ brings with it an increase in the CIV and CXV threshold velocities. Again restricting our analysis to $τ̂e≥τ̂cap$⁠, Eqs. (40) and (41) indicate that the CXV and CIV threshold velocities are bounded by

and

respectively. In both cases, the threshold velocity increases with the ratio $τ̂ie/τ̂e$ because of a corresponding decrease in the fraction of ion energy transferred to the electron population. This result agrees with Formisano et al.,³ who found that ion–electron energy transfer inefficiencies increased the CIV threshold velocity in cometary plasmas.

In this paper, we theoretically investigated the interaction between a high-speed neutral gas flow and dipole magnetized plasma. This problem represents a departure from classical CIV theories, which generally consider an isotropic plasma of infinite extent flowing against neutral gas within a uniform, open magnetic field configuration. Single particle trajectories were examined as a function of the location at which neutral particles undergo ionizing reactions, either via charge exchange with plasma ions or electron impact ionization. This trajectory analysis provided a method to incorporate non-uniform magnetic field effects within the fluid theory typically used to analyze CIV. Specifically, by identifying trapped and deflected ion orbits and integrating spatially, we were able to derive equations for the transfer of mass, momentum, and energy from the neutral gas flow to the dipole plasma. The resulting equations were incorporated into a simple model for the conservation of mass and energy in steady-state, from which we identified mode transitions indicative of CIV-like effects.

The main contributions of our analysis can be summarized as follows:

1. There exists a critical magnetic flux surface, $ψ*$⁠, within which newly formed ions enter into closed trajectories around the magnetic dipole. In general, ions formed outside of $ψ*$ either bypass or deflect from the dipole magnetic field. We found that $ψ*∼ρL$⁠, where ρ_L is the characteristic Larmor radius of newly formed ions normalized by the magnet radius. Furthermore, the volume bounded by a particular flux surface increases as ψ decreases. Therefore, the scaling of $ψ*$ with ρ_L describes how the spatial domain of trapped-ion orbits changes with the magnet radius, magnetic field strength, and neutral flow mass and velocity.

2. Equations for the transfer of mass, momentum, and energy from the neutral flow to the dipole plasma were derived as a function of two dimensionless quantities: ρ_L and $ζtot$⁠, as defined in Eqs. (6) and (11), respectively. Here, ρ_L dictates the spatial distribution of captured and deflected ion orbits (described above), while $ζtot$ determines the probability that a streaming neutral particle will undergo an ionizing reaction at a particular location. It was found that mass and energy transfer increase monotonically as ρ_L decreases. A critical value of $ζtot$ was discovered to maximize mass and energy transfer to the plasma for a given ρ_L. The appearance of a critical value is explained by the fact that in a highly reactive plasma (large $ζtot$⁠) a majority of neutral flow particles are ionized and deflected prior to reaching the ion-trapping volume (i.e., wake shadowing). We determined that momentum transfer increases monotonically with $ζtot$ and approaches an asymptotic value (i.e., drag saturation), which depends only on ρ_L and is roughly equal to 70% of the drag that would be produced if the neutral flow was instead fully ionized.

3. The dependence of the particle capture rate on plasma density drives a mode transition from a regime where injected mass and power are required to sustain the plasma (I regime) to a regime where the plasma can be sustained entirely by the neutral flow. The flow-sustained regime can be further divided into two separate modes. In the CX regime, diffusion of charge-exchange neutrals enables flow energy to sustain the plasma without overwhelming the electron population with ionization energy losses. In the CIV regime, flow energy is sufficiently large to ionize and heat the entire population of captured neutral particles. Transitions between these regimes were derived as a function of the electron confinement timescale, ion–electron energy transfer timescale, neutral flow particle capture timescale, and the ratio of the effective ionization energy to the kinetic energy of neutral flow particles. Threshold velocities were derived for each mode transition that were found to be similar to Alfvén's critical ionization velocity, but with additional dependencies on the timescales listed above.

The results presented here are general in the sense that we do not specify physical models for diffusion, charge exchange and ionization reactions, ionization energy losses, and ion–electron energy transfer. One limitation that results from this generality is that the time-evolution of the density and temperature of the different species cannot be determined self-consistently. This limitation is addressed in Paper II,²⁵ where we use the results of this current paper to derive a global model for the plasma/flow interaction that includes physical models for the various diffusion, reaction, and energy transfer processes for each species. This self-consistent model will ultimately provide a framework for us to examine the characteristics of, and transitions between, the different flow interaction regimes identified in this work.

We thank Dr. Robert Moses for the insightful conversations related to this research and its applications to plasma aerocapture. C. L. Kelly's effort was supported by the NASA Space Technology Research Fellowship No. 80NSSC18K1191.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Simple analytical approximations for the normalized particle capture rate [Eq. (14)] and normalized drag [Eq. (18)] are included here as a function of the dimensionless quantities $ζtot$ and ρ_L. The functions were obtained for α = 4 and $ψ̂r=1$⁠.

For $ρL≲1$ and $ζtot>0.1$⁠, the particle capture integral is well-approximated by the function