Perturbation analysis of ionization oscillations in Hall effect thrusters

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⁷ and the electric circuit.⁸ 

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.⁸ Yamamoto et al.¹⁰ have extended the predator-prey model proposed by Fife⁵ 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.¹² 

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.

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,⁵ 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′ exp(−iωt)$⁠, where Q₀ and $Q′$ 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_iU_n)^1∕2/L, where U_i = (eV_D/M_i)^1∕2 is a function of the discharge voltage V_D and M_i 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 L.

There are two assumptions that need to be reconsidered. First, as pointed out by Barral and Ahedo,¹³ 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.

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_ch is the channel length, and 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_ch, and R_Δ. Thus, L is a function of T_e and the geometry of the channel.

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_e because ξ_ion = ξ_ion(T_e). Equation (9) provides the first condition to have a steady-state plasma generated in the channel, i.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_i (via discharge voltage) or decreasing N_int (via anode mass flow) increases T_e,min. Thus, the condition of having a steady-state plasma becomes more severe.

Figure 2 shows T_e,min for different N_int and V_D when the expression for the ionization rate coefficient ξ_ion of ground-state xenon atoms proposed in Ref. 15 is used

where A = −1.0 × 10⁻²⁴, B = 6.386 × 10⁻²⁰, and C = 12.13. The configuration of a Stationary Plasma Thruster (SPT)-100 ML thruster³ 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_D = 300 V and N_int = 1.6 × 10¹⁹m⁻³. 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.³ A steady-state discharge plasma cannot be sustained when T_e,min is larger than the maximum T_e due to space charge saturation.

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_i increases as the discharge voltage increases, so L/L_ch increases. Note that L/L_ch in Eq. (7) is not dependent on the neutral atom density. In Fig. 3, T_e,min to sustain a steady-state plasma is obtained assuming N_int = 1.6 × 10¹⁹m⁻³.

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

for $Q→=[Ni′,Nn′]T$⁠. The wave frequency can be calculated by solving the determinant of the matrix: $−ω2−iωNintNi,0ξion/(Nint−Nn,0)+Nn,0Ni,0ξion2=0$⁠. The frequency can be written as the real frequency and the growth rate as ω = ω_r + iγ, where ω_r and γ are always real numbers, assuming the perturbation to follow $exp(−iω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⁵rad/s for T_e,min ( = 12.4 eV) < T_e < 30 eV when V_D = 300 V. Although Eq. (13) becomes identical to the predator-prey model (Eq. (4)) when γ = 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.

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_e. In addition to the continuity equations for ions and neutral atoms, the perturbations in ion momentum and electron temperature are now taken into consideration.

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_eU_eE is the energy gain due to Joule heating, S_wall = N_eϵ_w(T_e) ν_w(T_e) is the energy loss due to the channel walls, and $Scoll=∑jNeNnξj(Te)ϵj$ denotes inelastic collisions, 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_e 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 × B drift is strong.⁹ 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⁶ is given by

where σ is the effective secondary electron emission (SEE) rate. The energy loss to the wall¹⁵ is written as

where $ϕw=−Te log[(1−σ)/(2π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_coll is further assumed as S_coll = χN_iN_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_e) for the range of T_e that is considered, e.g., 5 < T_e < 30 eV.

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→=[Ni′,Nn′,Ui′,Te′]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_e instead of the operational parameters such as B and V_D. The balance between the loss mechanisms and the Joule heating determines T_e,0 dynamically and non-locally in actual HETs.⁹ Also, κ = 1 is assumed. The determinant of Eq. (22) can be written as –ω⁴+ iω³F₁ + ω²F₂ + iωF₃ + F₄ = 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.

In addition to the values used in Sec. III, the SEE model proposed in Refs. 6 and 9 is used. $σ=max(σ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_e = 24.6 eV. The nominal conditions are N_int = 1.6 × 10¹⁹m⁻³ and V_D = 300 V. In order to investigate the stability of ionization oscillations in the SPT-100 thruster, U_e is varied from 0 to –U_i and T_e is varied from T_e,min = 12.4 eV to 25 eV. U_e is related to the current utilization efficiency η_c ≡ I_i/I_D, where I_i is the ion current and I_D is the total discharge current. Assuming only singly charge ions and quasineutrality,

From the literature, 0.72 < η_c < 0.86 for the NASA-173Mv2 thruster,¹⁶ 0.70 < η_c < 0.92 for the 6 kW laboratory Hall thruster at the University of Michigan,¹⁷ and 0.68 < η_c < 0.85 for the BHT-2000 thruster¹⁸ 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 η_c 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 U_e considered is from –U_i to 0, which corresponds to 0.5 < η_c < 1 from Eq. (24).

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_e particularly in large $|Ue/Ui|$⁠, or small η_c. However, for small $|Ue/Ui|$⁠, the increase in ω_r as increasing 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,⁹ which is compared with the experiments of the SPT-100 (Ref. 3) and H6 thrusters.⁷ 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_e approaches the space charge limited regime as the discharge oscillation stabilizes when increasing the magnetic field strength 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_e (or η_c) and 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.²³ This phenomenon can be related to the change in U_e from the plume into the channel.

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).

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⁷s⁻¹. 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_e ≈ 24.6 eV.

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_e is the diffusion coefficient and ν_eff is the effective collision frequency. In this study, Λ has contributions from the electron convective heat flux, wall heat flux, and inelastic collisions.

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.⁶ 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.

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²⁷ and can be applied for traveling plasma waves.

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 η_c) and T_e. 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.

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.