A perturbation analysis of ionization oscillations, which cause low frequency oscillations of the discharge plasma, in Hall effect thrusters is presented including the electron energy equation in addition to heavy-species transport. Excitation and stabilization of such oscillations, often called the breathing mode, are discussed in terms of the growth rate obtained from the linear perturbation equations of the discharge plasma. The instability induced from the ionization occurs only when the perturbation in the electron energy is included while the neutral atom flow contributes to the damping of the oscillation. Effects of the electron energy loss mechanisms such as wall heat loss, inelastic collisions, and convective heat flux are discussed. It is shown that the ionization oscillations can be damped when the electron transport is reduced and the electron temperature increases so that the energy loss to the wall stabilizes the ionization instability.

## I. INTRODUCTION

Low frequency ionization oscillations in the range of 10–30 kHz in electric propulsion devices have been observed for over 40 yr in experiments^{1–3} and numerical simulations.^{4–6} Such plasma oscillations result in the oscillation of the discharge current, often referred to as the breathing mode, which is deleterious for the thruster efficiency^{7} and the electric circuit.^{8}

Recent studies on mode transition of Hall effect thrusters (HETs)^{7,9} showed that the discharge oscillation is damped when varying the plasma parameters. Although several theoretical frameworks have been developed,^{5,6,10} there has been no complete theoretical justification and explanation of the discharge and plasma oscillations.^{8} Yamamoto *et al.*^{10} have extended the predator-prey model proposed by Fife^{5} and allowed perturbation in the electron current. Barral *et al.* have developed a theoretical framework in which the growth rate of the discharge current is inherently assumed so that the criteria for growth and damping effects cannot be discussed.^{6,11} Peradzyński *et al.* have developed a simplified theory and shown implications that instability may be excited when including perturbation in the electron temperature.^{12}

In the present theoretical analysis, the stability and excitation criteria of ionization oscillation in HETs are discussed including the perturbation in the electron temperature. First, the predator-prey model is revisited in Sec. II. In Sec. III, the heavy species transport equations are derived and it is shown that ionization oscillations are unconditionally damped when only accounting for the heavy species transport. In Sec. IV, a complete perturbation theory is developed including the ion and neutral continuity equations, the ion momentum equation, and the electron energy equation. It is shown that an undamped solution appears when the perturbation in electron energy is taken into consideration.

## II. PREDATOR-PREY MODEL

The mechanism of ionization oscillations in HETs has been explained by insufficient neutral flow.^{4,5,13} Theory has been formulated using the transport of heavy species, i.e., neutral atoms and ions. In the HET community, the predator-prey model,^{5} of which the simplest form is also known as the Lotka-Volterra model, has been widely used and it has often been attempted to obtain a more correct scaling of the ionization length.^{14,15} In this section, we review the predator-prey model.

The model assumes that the plasma is contained in an ionization box whose length is *L* and the spatial variation is neglected, i.e., *∂*/*∂x* ∼ 1/*L*. The ionization rate coefficient *ξ _{ion}* and the velocities of ions and neutral atoms are constant in time. There is no ion flux entering the box and no neutral flux escaping the box. Using these assumptions, the continuity equations are written as

where *N* and *U* are the number density and mean velocity, and the subscripts *i* and *n* denote ions and neutral atoms, respectively. To study the linear perturbation, a quantity follows the form: $Q=Q0+Q\u2032\u2009exp(\u2212i\omega t)$, where *Q*_{0} and $Q\u2032$ are equilibrium and perturbation quantities, respectively. Using this first-order perturbation equation, an equation that exhibits a harmonic oscillation can be derived as

and the harmonic oscillator frequency is given by

Using the zeroth-order equilibrium condition, one obtains *ω* = (*U _{i}U_{n}*)

^{1∕2}/

*L*, where

*U*= (

_{i}*eV*/

_{D}*M*)

_{i}^{1∕2}is a function of the discharge voltage

*V*and

_{D}*M*is the ion mass. Although ionization oscillations depend on other parameters such as the mass flow, magnetic field topology, and channel wall materials, it is assumed that these are incorporated into the ionization length

_{i}*L*.

There are two assumptions that need to be reconsidered. First, as pointed out by Barral and Ahedo,^{13} the neutral inflow rate in the actual HET operation is fixed at the anode and the outflow varies, so Eq. (2) is inaccurate. The meaning of Eq. (2) is simply to satisfy the predator-prey model. Second, *L* is not well-defined and is often assumed to be on the order of the channel length. Therefore, a more complete model must be considered to study the mechanism of the ionization oscillation.

## III. HEAVY SPECIES TRANSPORT

In order to construct the correct theory of ionization oscillations, we first define the geometry. Figure 1 shows a simplified HET schematic used in the present model. The discharge plasma is assumed to be confined inside the discharge channel. Although the one-dimensional flow in the axial direction is of interest, the radial plasma diffusion also plays an important role. One can take an approach similar to a finite-volume method in which the state variables are volume averaged and the fluxes at interfaces are modeled. Without employing the undefined ionization length, ion and neutral continuity equations are more correctly given by

where *N _{int}* is the number density of neutral atoms at the anode,

*U*

_{i}_{,}

_{w}= (

*eT*/

_{e}*M*)

_{i}^{1/2}is the ion acoustic speed,

*L*is the channel length, and

_{ch}*R*

_{Δ}is the channel width. Ion diffusion toward the channel walls is taken into account assuming that Bohm's condition is satisfied at the sheath edge near the channel walls. Any spatial variations are neglected inside the box due to the zero-dimensional assumption. Since there is no radial diffusion for the neutral atoms, their transport can be described only by accounting for the axial transport. Note that the ionization length

*L*can be described using geometric parameters by relating Eq. (5) and the ion continuity equation in the predator-prey model (Eq. (1))

Ionization length is now defined using the geometric parameters and plasma properties for the first time. Equation (7) shows that *L* decreases when *T _{e}* increases since the plasma diffusion becomes stronger, i.e.,

*U*

_{i}_{,}

_{w}increases for given

*U*,

_{i}*L*, and

_{ch}*R*

_{Δ}. Thus,

*L*is a function of

*T*and the geometry of the channel.

_{e}From Eqs. (5)–(7), the equilibrium densities of ions and neutral atoms are given by

where subscript 0 denotes the equilibrium condition. Both *N _{n}*

_{,0}and

*N*

_{i}_{,0}are a function of

*T*because

_{e}*ξ*=

_{ion}*ξ*(

_{ion}*T*). Equation (9) provides the first condition to have a steady-state plasma generated in the channel, i.e.,

_{e}*N*

_{i}_{,0}> 0 or

*N*>

_{int}*N*

_{n}_{,0}must be satisfied. From Eq. (8),

where the left hand side is the ionization frequency when ionizing all the inflow neutral gas and the right hand side is the characteristic frequency due to ion acceleration. For instance, ionization needs to be sufficiently large to sustain a steady-state plasma when the ion outflow is fast. Given all the other parameters, the minimum *T _{e}* (=

*T*

_{e}_{,}

_{min}) can be calculated when the two terms in Eq. (10) are equal. Increasing

*U*(via discharge voltage) or decreasing

_{i}*N*(via anode mass flow) increases

_{int}*T*

_{e}_{,}

_{min}. Thus, the condition of having a steady-state plasma becomes more severe.

Figure 2 shows *T _{e}*

_{,}

_{min}for different

*N*and

_{int}*V*when the expression for the ionization rate coefficient

_{D}*ξ*of ground-state xenon atoms proposed in Ref. 15 is used

_{ion}where *A* = −1.0 × 10^{−24}, *B* = 6.386 × 10^{−20}, and *C* = 12.13. The configuration of a Stationary Plasma Thruster (SPT)-100 ML thruster^{3} is assumed: *L _{ch}* = 2.5 cm and

*R*

_{Δ}= 2 cm. For instance, minimum electron temperature required to have a steady-state plasma is

*T*

_{e}_{,}

_{min}= 12.4 eV for

*V*= 300 V and

_{D}*N*= 1.6 × 10

_{int}^{19 }m

^{−3}. As shown in Fig. 2, a higher electron temperature is required to sustain the discharge plasma for higher discharge voltage and/or higher mass flow rate. In addition, the maximum electron temperature can vary due to the space charge limited sheath depending on the wall material.

^{3}A steady-state discharge plasma cannot be sustained when

*T*

_{e}_{,}

_{min}is larger than the maximum

*T*due to space charge saturation.

_{e}Figure 3 shows the ionization length as a function of discharge voltage and electron temperature under the same geometric configuration as in Fig. 2. As the electron temperature increases, the ionization length decreases since the plasma diffusion is enhanced due to higher electron temperature via *U _{i}*

_{,}

_{w}in Eq. (7). On the other hand, the axial ion velocity

*U*increases as the discharge voltage increases, so

_{i}*L*/

*L*increases. Note that

_{ch}*L*/

*L*in Eq. (7) is not dependent on the neutral atom density. In Fig. 3,

_{ch}*T*

_{e}_{,}

_{min}to sustain a steady-state plasma is obtained assuming

*N*= 1.6 × 10

_{int}^{19 }m

^{−3}.

The first-order perturbation equations of Eqs. (5) and (6) can be described as a system of equations

for $Q\u2192=[Ni\u2032,Nn\u2032]T$. The wave frequency can be calculated by solving the determinant of the matrix: $\u2212\omega 2\u2212i\omega NintNi,0\xi ion/(Nint\u2212Nn,0)+Nn,0Ni,0\xi ion2=0$. The frequency can be written as the real frequency and the growth rate as *ω* = *ω _{r}* +

*iγ*, where

*ω*and

_{r}*γ*are always real numbers, assuming the perturbation to follow $exp(\u2212i\omega t)$. The solution of Eq. (12) is given by

Under the same condition as in Fig. 3, the real part of the solution is 0 < *ω _{r}* < 2 × 10

^{5 }rad/s for

*T*

_{e}_{,}

_{min}( = 12.4 eV) <

*T*< 30 eV when

_{e}*V*= 300 V. Although Eq. (13) becomes identical to the predator-prey model (Eq. (4)) when

_{D}*γ*= 0, Eq. (14) shows that the growth rate of the oscillation is always negative, i.e.,

*γ*< 0 since

*N*>

_{int}*N*

_{n}_{,0}to have a steady-state plasma. Therefore, only solving the two continuity equations shows that the oscillation is always damped.

It is also worth noting that ionization oscillations can be unstable if *T _{e}* <

*T*

_{e}_{,}

_{min}due to the inability to sustain the plasma discharge, which may occur when electron heating is insufficient. However, such oscillations are different from the breathing mode. The breathing mode can result in stable operation but is highly oscillatory since the ionization oscillations occur as an unstable solution from the linear perturbations around a steady-state equilibrium discharge plasma.

## IV. IONIZATION OSCILLATION THEORY

It is shown from Eqs. (13) and (14) that ionization oscillations always damp when there is no perturbation in the ionization rate coefficient *ξ _{ion}*, or equivalently the electron temperature

*T*. In addition to the continuity equations for ions and neutral atoms, the perturbations in ion momentum and electron temperature are now taken into consideration.

_{e}First, using the relation in Eq. (7), the ion momentum equation in the axial direction is given by

where *E* is the electric field. The pressure term is neglected since the ions are assumed to be cold, i.e., *T _{i}* ≈ 0. Next, the electron energy equation is taken into consideration

where *S _{joule}* = –

*N*is the energy gain due to Joule heating,

_{e}U_{e}E*S*=

_{wall}*N*

_{e}*ϵ*

_{w}(

*T*)

_{e}*ν*(

_{w}*T*) is the energy loss due to the channel walls, and $Scoll=\u2211jNeNn\xi j(Te)\u03f5j$ denotes inelastic collisions,

_{e}*j*, that contribute to energy loss. For simplicity, the contributions from electron pressure and heat conductivity are neglected. The convective energy flux at the exit is carried into the system, since the electron velocity

*U*is negative in HETs and it is assumed that the electron energy is sufficiently low at the anode so that there is no convective flux towards the anode. In the present model, it is also assumed that the electron mobility and the electric field are constant in time and the electron kinetic energy is negligible although the effect of the kinetic energy should be included in HETs where the

_{e}*E*×

*B*drift is strong.

^{9}Finally, a quasineutral assumption is used

*N*=

_{i}*N*.

_{e}In order to investigate the perturbation in the right hand side of Eq. (16), we introduce models for each term. The wall collision frequency^{6} is given by

where *σ* is the effective secondary electron emission (SEE) rate. The energy loss to the wall^{15} is written as

where $\varphi w=\u2212Te\u2009log[(1\u2212\sigma )/(2\pi me/Mi)1/2]$ is the sheath potential and *m _{e}* is the electron mass. Only singly charge ionization and excitation from the ground state atoms are considered.

*S*is further assumed as

_{coll}*S*=

_{coll}*χN*(

_{i}N_{n}ξ_{ion}*T*)

_{e}*ϵ*

_{ion}, where

*χ*is the ionization cost that includes the effect of excitation and

*ϵ*

_{ion}is the ionization energy loss. Since the expression for the ionization rate coefficient is complicated as can be seen from Eq. (11), the perturbation form of that is assumed to follow

*ξ*(

_{ion}*T*) ≈

_{e}*ξ*

_{ion}_{,0}(

*T*/

_{e}*T*

_{e}_{,0})

^{κ}for simplicity, where

*T*

_{e}_{,0}is the equilibrium electron temperature and

*ξ*

_{ion}_{,0}is the ionization rate coefficient for

*T*

_{e}_{,0}given in Eq. (11). This form has been chosen due to the monotonic increase in

*ξ*(

_{ion}*T*) for the range of

_{e}*T*that is considered, e.g., 5 <

_{e}*T*< 30 eV.

_{e}For simplicity, the perturbations of *L*, *χ*, and *σ* are neglected, i.e., frozen. The first-order terms of the quantities above are given by

Using the quasineutral assumption, the perturbation equations for $Q\u2192=[Ni\u2032,Nn\u2032,Ui\u2032,Te\u2032]T$ can be written as

where *N _{i}*

_{,0}> 0 to have a steady-state plasma, and

is the effective electron energy relaxation frequency that has contributions from the energy loss mechanisms such as wall loss, inelastic collisions, and the convective energy flux.

In the present study, *T _{e}*

_{,0}is taken as a parameter as well as

*U*instead of the operational parameters such as

_{e}*B*and

*V*. The balance between the loss mechanisms and the Joule heating determines

_{D}*T*

_{e}_{,0}dynamically and non-locally in actual HETs.

^{9}Also,

*κ*= 1 is assumed. The determinant of Eq. (22) can be written as –

*ω*

^{4 }+

*iω*

^{3}

*F*

_{1}+

*ω*

^{2}

*F*

_{2}+

*iωF*

_{3}+

*F*

_{4}= 0, where the coefficients

*F*(

_{k}*k*= 1, 2, 3, and 4) are functions of the plasma parameters. Since it is difficult to derive the solutions analytically, we analyze the stability criteria using the representative values of the flow in Sec. V.

## V. RESULTS

In addition to the values used in Sec. III, the SEE model proposed in Refs. 6 and 9 is used. $\sigma =max(\sigma max,Te/25)$, where *σ _{max}* = 0.986 is the space charge limited SEE coefficient. A space charge limited sheath is assumed to form at

*T*= 24.6 eV. The nominal conditions are

_{e}*N*= 1.6 × 10

_{int}^{19 }m

^{−3}and

*V*= 300 V. In order to investigate the stability of ionization oscillations in the SPT-100 thruster,

_{D}*U*is varied from 0 to –

_{e}*U*and

_{i}*T*is varied from

_{e}*T*

_{e}_{,}

_{min}= 12.4 eV to 25 eV.

*U*is related to the current utilization efficiency

_{e}*η*≡

_{c}*I*/

_{i}*I*, where

_{D}*I*is the ion current and

_{i}*I*is the total discharge current. Assuming only singly charge ions and quasineutrality,

_{D}From the literature, 0.72 < *η _{c}* < 0.86 for the NASA-173Mv2 thruster,

^{16}0.70 <

*η*< 0.92 for the 6 kW laboratory Hall thruster at the University of Michigan,

_{c}^{17}and 0.68 <

*η*< 0.85 for the BHT-2000 thruster

_{c}^{18}across a range of parameters. While these thrusters exhibit 50%–70% total efficiency, SPT-100 thrusters exhibit total efficiency of 35%–50%,

^{19–21}so it can be extrapolated that

*η*becomes worse than that for NASA-173Mv2 or BHT-2000 thrusters. In fact, it has been reported in Ref. 22 that the current utilization efficiency for a SPT-100 thruster is

_{c}*η*≈ 0.6 at optimal operation and can decrease down to

_{c}*η*≈ 0.4 at lower discharge voltage. In this study, the range for

_{c}*U*considered is from –

_{e}*U*to 0, which corresponds to 0.5 <

_{i}*η*< 1 from Eq. (24).

_{c}As shown in Fig. 4, there is an undamped oscillation solution, i.e., *γ* > 0, when including the perturbation in electron temperature, whereas the result without the electron energy perturbation has only damped solutions as shown in Eq. (14). Note that only including the ion momentum equation in addition to the two continuity equations results unconditionally in damping, similar to Eq. (14) in Sec. III. Therefore, the electron energy perturbation is required for the excitation of ionization oscillations. The corresponding real frequency in the unstable region is shown in Fig. 5. It can be seen that *ω _{r}* increases as

*T*

_{e}_{,0}increases but decreases after there is a maximum peak for a constant

*U*particularly in large $|Ue/Ui|$, or small

_{e}*η*. However, for small $|Ue/Ui|$, the increase in

_{c}*ω*as increasing

_{r}*T*

_{e}_{,0}is rather monotonic. With these theoretical predictions of the growth rates, we can discuss the cause of the discharge oscillation mode transition and the direct cause of discharge oscillations in HETs. For convenience, the two stable regions are labeled as regions I and II as shown in Figs. 4 and 5.

The present theory supports the observation in a recent numerical simulation of the mode transition,^{9} which is compared with the experiments of the SPT-100 (Ref. 3) and H6 thrusters.^{7} As shown in Figs. 4 and 5, increased wall losses via increased *T _{e}*

_{,0}result in damping of the ionization oscillations in region II. From the simulation results in Ref. 9, $|Ue|$ is reduced and

*T*approaches the space charge limited regime as the discharge oscillation stabilizes when increasing the magnetic field strength

_{e}*B*. Since $|Ue|$ is usually a decreasing function of

*B*, the increase in

*B*results in the decrease in $|Ue|$. The mode transition observed in Ref. 9 corresponds to the mode transition from the linearly unstable region to region II in Figs. 4 and 5. However, in order to directly compare experiments and the present theory, the internal electron temperature profile and the electron current need to be known in detail. The values of

*U*(or

_{e}*η*) and

_{c}*T*

_{e}_{,0}need to be examined in the experiments since the mode transition of the ionization oscillation is determined by the electron transport. For instance, a small change in the cathode mass flow may also change the mode transition.

^{23}This phenomenon can be related to the change in

*U*from the plume into the channel.

_{e}The cause of the ionization oscillations in HETs has often been explained by insufficient neutral atom flow using the predator-prey model,^{4,8,13} in which electron transport has been neglected. In the present study, it is shown that the ionization oscillations in HETs are strongly related to the electron energy perturbations, in particular, the electron energy relaxation frequency Λ as shown in Fig. 6. This framework can be used to estimate the stability region of ionization oscillations in HETs, although it is difficult to obtain a simple analytic formula due to the complexity of the equations as shown in Eq. (22).

## VI. DISCUSSIONS

### A. Electron energy relaxation frequency

As can be seen from Figs. 4–6, smaller Λ corresponds to larger growth rate. However, region I predicts damping of the ionization oscillations at small *T _{e}* and small $|Ue|$. This is possibly due to small ionization rate and the neutral atom flow. Note the region where

*T*<

_{e}*T*

_{e}_{,}

_{min}is unstable since the steady-state plasma is not satisfied. The other damped solution (region II) appears when the heat loss due to the near-wall sheath increases so that Λ increases as

*T*

_{e}_{,0}increases. The relaxation of the electron energy is too fast in comparison to the other frequencies in the system so that there will be no low frequency ionization oscillations. Most importantly, Fig. 6 shows that the stability condition for regime II corresponds to Λ ≈ 1 × 10

^{7 }s

^{−1}. This boundary is almost a straight line indicating that

*T*

_{e}_{,0}mainly plays an important role for the stability. The electron temperature for

*γ*= 0 is approximately 23 eV, which is slightly below the electron temperature for the space charge limited sheath,

*T*≈ 24.6 eV.

_{e}Similar concepts as the electron energy relaxation frequency in the present study have been used in the low temperature plasma community to discuss the effect of nonequilibrium electrons.^{24,25} However, the discharge plasma in the HETs and other low temperature plasmas can be significantly different. The discharge plasma in HETs is more complicated since there is electron flow, heat loss through the channel walls, and inelastic collisions, whereas the electrons mainly undergo diffusion and inelastic collisions in low temperature plasma systems and the energy relaxation length is often defined as *λ*_{ϵ} ∼ (*D _{e}*/

*ν*)

_{eff}^{1∕2}, where

*D*is the diffusion coefficient and

_{e}*ν*is the effective collision frequency. In this study, Λ has contributions from the electron convective heat flux, wall heat flux, and inelastic collisions.

_{eff}### B. Nonlinear saturation of ionization oscillations

The non-sinusoidal shape of the discharge current can be explained from the growth rate of the linearized oscillation wave. Some examples of the discharge current oscillations can be found in Refs. 3, 7, and 9. When the growth rate of the linear perturbation of the ionization oscillation is positive, the perturbation in all quantities will grow exponentially, and the perturbed values may reach some threshold value, say, *N _{i}*

_{,0}> 0 and

*N*>

_{int}*N*

_{n}_{,0}> 0. This results in so-called avalanche ionization.

^{6}After the avalanche ionization, the decay is also exponential. Therefore, the shape of the discharge current becomes an exponential rise and then exponential decay. In addition, as the avalanche ionization is a fast process, the plasma needs to wait for the depleted neutral atoms to fill in the ionization box, which is the slowest process in the system.

The growth rate of the linear perturbations in the present study may not exactly predict the behavior of discharge oscillations in HETs since the plasma parameters of the ionization oscillations can change significantly and the oscillations are likely to be a nonlinear phenomenon. However, the significance of this analysis is the observation of the unstable solutions from the linear perturbations in the discharge plasma of HETs. Linear instability grows and nonlinearity can play an important role, whereas a linearly stable system will damp the ionization oscillations.

### C. Comparison to experiments

Discharge oscillations also have been observed when further increasing the magnetic field after a stable discharge mode. In Ref. 26, it was discussed that a different strong ionization oscillation may be related to plasma diffusion towards the anode. The plasma diffusion or acceleration to the anode has not been taken into consideration in this model. Including this effect into the model is reserved for future work. In addition, the theory can be applied when changing the propellant. For instance, lighter noble gases such as argon and krypton are used in some thrusters. The mass of the heavy species can alter the ionization length in Eq. (7), the inelastic collision rate coefficients in Eq. (11), and the wall heat flux in Eqs. (17) and (18).

The framework of Eq. (22) theoretically can be applied to any plasma to investigate the ionization oscillation. In this study, it is assumed that the wall loss is significant in the SPT-type HETs because the ionization region is known to be inside the discharge channel. In addition, the formulations are based on the radial and axial transport assuming that there is no azimuthal component. The present zero dimensional framework may also be extended including the spatial variation in the azimuthal direction to investigate the spatial effects including azimuthally rotating spokes^{27} and can be applied for traveling plasma waves.

## VII. SUMMARY

A complete perturbation theory of the ionization instability in HETs is developed including the ion and neutral continuity, ion momentum, and electron energy equations. There are three observations in comparison with the previous predator-prey type models.

Without perturbations in electron energy, the ionization oscillations are always damped. The fixed neutral atom flow from the anode contributes to damping of the oscillation.

The ionization length

*L*, which is frequently used in the literature, is defined using the plasma and geometric parameters.Minimum electron temperature to sustain a steady-state discharge plasma is obtained from the continuity equations of ion and neutral atom.

It is further shown that the ionization oscillations in HETs are caused by the perturbation in the electron energy. This significantly advances the understanding of the breathing mode, for which only heavy species transport has been considered in the literature. Inclusion of the electron energy perturbation allows an undamped solution for the linear perturbations of the ionization oscillation wave. The parameters discussed in the present model are *U _{e}* (or

*η*) and

_{c}*T*. In addition, the present theory suggests the stabilization mechanism of ionization oscillations. Reduced electron transport and increased electron temperature yields a transition from an undamped oscillation mode to a stabilized mode. This theoretical observation supports recent numerical simulation results in Ref. 9. In order to investigate the ionization oscillations in experiments of the HETs, the electron transport properties such as the electron current and electron temperature must be measured.

_{e}## ACKNOWLEDGMENTS

The authors acknowledge the support for this work provided by the U.S. Department of Energy Office of Science, Fusion Energy Sciences Program, Grant No. DESC0001939, and the Air Force Research Laboratory through the Michigan/Air Force Center of Excellence in Electric Propulsion (MACEEP) Grant No. F9550-09-1-0695.