Three-dimensional coupling of electron cyclotron drift instability and ion–ion two stream instability

Low-temperature, partially magnetized plasmas are used for a wide range of applications, including material processing, space propulsion, and beam sources. In these plasmas, electrons are magnetized, whereas ions are unmagnetized. While crossed electric and magnetic fields are used to confine electrons, the electron drifts, such as the $E×B$ drift (where E is the applied electric field, and B is the applied magnetic field), cause plasma instabilities, leading to plasma turbulence and anomalous electron transport across magnetic fields.

The main kinetic instabilities that occur in partially magnetized plasmas include the electron cyclotron drift instability (ECDI)^1,2 and modified two-stream instability (MTSI)^3,4 driven by the electron $E×B$ drift. In laboratory cross field discharge sources, coherent Thomson scattering experiments have shown the existence of plasma waves driven by such kinetic instabilities in the direction of $E×B$ drift.^5,6 A number of kinetic simulations have recently revisited these instabilities and their effect on cross field electron transport.^7–11 In addition to the plasma waves in the $E×B$ direction, Tsikata et al.¹² have found experimental evidence of plasma waves that propagate in the direction of the applied electric field (E) resulting from an ion–ion two-stream instability (IITSI) driven by the relative drift velocity between singly and doubly charged ions in the cross field plasma discharge. Hara and Tsikata¹³ further observed using a particle-in-cell (PIC) simulation that the coupling of ECDI and IITSI in a two-dimensional (2D) configuration results in microturbulence, which enhances electron diffusion across the magnetic fields. Furthermore, the cross field electron transport is found to be affected by the magnitude and wavelength of the plasma waves.^14,15

ECDI occurs due to coupling of electron Bernstein waves and unmagnetized ions. The Bernstein mode in three-dimensional (3D), i.e., when spatial fluctuations are accounted for in the direction parallel and perpendicular to the magnetic field, results in the use of the plasma dispersion function assuming that the electrons follow a Maxwellian velocity distribution function (VDF), while the 1D and 2D assumptions reduce the plasma dispersion function to resonance-type modes that can be described using polynomials.¹⁶ Cavalier et al.¹⁷ developed a numerical method to obtain the solution of the 3D ECDI dispersion relation for singly charged ions. For such cases, there are two solutions that are complex conjugates, allowing for an iterative method to obtain the most unstable root. However, when multiply charged ion streams are considered, IITSI occurs in addition to ECDI, and the resulting modes are no longer a superposition of the two instabilities. As shown later, when ECDI and IITSI coexist, more than one unstable root exists, making the iterative method infeasible. Thus, a generalized numerical solver to calculate an arbitrary number of unstable roots of a kinetic dispersion relation is needed.

In this paper, we investigate the three-dimensional coupling of ECDI and IITSI in a partially magnetized plasma through the development of a numerical technique for finding complex roots of the coupled dispersion relation. A theoretical examination of ECDI and IITSI and the root-finding algorithm are presented in Sec. II. Section III discusses the dispersion relation for singly charged ions and magnetized electrons, validating the root-finding algorithm with published results. In Sec. IV, the dispersion relation for singly and doubly charged ions and magnetized electrons is presented. Here, the coupled ECDI and IITSI yield more than one unstable root, which the generalized numerical algorithm can obtain.

Plasma instabilities can be discussed in the context of a linear perturbation to the equilibrium distribution function of the form $exp [−i(ωt−k·r)]$ with $ω=ωr+iγ$⁠, where ω is the wave frequency, k is the wave vector, t is the time, r is the position vector, ω_r is the real frequency, and γ is the growth rate. In general, the electrostatic dispersion relation can be obtained from Poisson's equation as

where χ_e is the electron susceptibility, and χ_i is the ion susceptibility. The susceptibilities can be constructed from the Vlasov equation, assuming a collisionless plasma.

In this study, the plasma is assumed to be partially magnetized, i.e., the electrons are magnetized, and the ions are non-magnetized. It is further assumed that (i) the ions are cold, and (ii) the electrons follow a Maxwellian VDF for the equilibrium condition, from which the linear perturbation evolves. These assumptions follow the physics of cross field discharges, e.g., Hall effect thrusters (HETs). The ions are typically generated via electron-impact ionization from neutral gas, making the plasma non-equilibrium (⁠$Ti≪Te$⁠, where T_i and T_e are the ion and electron temperatures) in the region where $E×B$ drift is largest.^18,19 Additionally, while a few numerical investigations have studied the effects of non-Maxwellian electron distributions,^20,21 a Maxwellian VDF is often considered for electrons in various computational models^22,23 and linear instability studies.^17,24–26 The effects of finite ion temperature (e.g., ion Landau damping²⁷) and non-Maxwellian VDFs on the kinetic instabilities will be reserved for future work.

The assumptions used in this study are consistent with previous ECDI studies.^12,13,17 In this section, the 3D (and 2D) ECDI theory is first reviewed to discuss the fixed-point iteration method by Cavalier et al.^17,28 Then, the dispersion relation due to the coupling of ECDI and IITSI, i.e., when another ion stream with different velocity exists, is introduced. Finally, the root finding algorithm and results are shown.

Taking into account the electron dynamics in directions parallel and perpendicular to the magnetic field while assuming that the electrons are magnetized and the equilibrium distribution function is a Maxwellian VDF, the electron susceptibility can be written as

where $k=|k|=(k⊥2+k∥2)1/2, k⊥=(kx2+ky2)1/2$ is the wavenumber perpendicular to the magnetic field, $k∥=kz$ is the wave number parallel to the magnetic field, $λD=(ϵ0kBTe/e2n0)1/2$ is the Debye length, $Ud$ is the drift velocity for the electrons (later we consider the $E×B$ drift), $Ui+$ is the singly charged ion bulk velocity, $vth,s=(kBTs/ms)1/2$ is the thermal velocity of species $s=e,i$⁠, and $ωce=eB0/me$ is the electron cyclotron frequency. Here, m_e and m_i are the electron and ion masses, ϵ₀ is the permittivity of free space, k_B is the Boltzmann constant, e is the elementary charge, n₀ is the quasineutral plasma density, $b=k⊥2rL2, rL=vth,e/ωce$ is the electron Larmor radius at the thermal velocity, Z is the plasma dispersion function given by $Z(ξ)=π−1/2∫−∞∞ exp (−t2)(t−ξ)−1dt$⁠, and I_n is the modified Bessel function of the first kind and order n. It is to be noted that the assumption that the electrons follow a Maxwellian VDF for the equilibrium condition results in the well-known electron Bernstein mode,¹⁶ as shown in Eq. (2). For non-Maxwellian VDFs, the electron susceptibility needs to be reevaluated.²⁹ 

Assuming cold ions (⁠$ω/k≫vth,i$⁠), the ion susceptibility term can be rewritten as

where the ion plasma frequency is given by $ωpi=(e2n0/ϵ0mi)1/2$⁠.

Using Eqs. (1)–(3), the ECDI dispersion relation of cold, singly charged ions, and magnetized electrons in a three-dimensional (3D) configuration¹⁷ can be written as

As shown later in Sec. II A 3, this 3D ECDI dispersion relation has only two roots which are complex conjugates. Thus, a fixed-point iteration scheme works to solve Eq. (4).

Another particular instability that arises from the 3D ECDI dispersion relation is called the modified two-stream instability (MTSI), which exists for finite, non-zero $k∥λD$⁠. The dispersion relation for MTSI is derived by taking the fluid limit,³ i.e., $krL≪1$⁠, of Eq. (4) in which only the n = 0 term is retained in the sum of the electron susceptibility

The motion of electrons both parallel and perpendicular to magnetic field lines is taken into account. It is to be noted that this dispersion relation is called the MTSI as it resembles the two-stream instability, in which two species counter stream. The final term in Eq. (5) can be considered that electrons have an effective mass of $(k2/k∥2)me$⁠, which can be large for small $k∥λD$⁠. Unstable modes exist that are largely perpendicular to the applied magnetic field with $k∥/k∼(me/mi)1/2$ such that the effective electron mass is on the order of the ion mass.³ Because of this fact, these waves cause electron heating parallel to the magnetic field that is on the order of ion heating perpendicular to the magnetic field.^9,30

In the limit of $k∥λD→0$⁠, the dispersion relation can be considered to be two-dimensional (in the $k⊥$ direction). The argument of the plasma dispersion function, ξ, approaches infinity, allowing the use of the asymptotic expansion for $Z(ξ)$⁠. Thus, Eq. (4) reduces to

It can be seen from Eq. (6) that resonances exist for all higher harmonics $n=1,2,…$ if $(ω−k·Ud)2=(nωce)2$⁠. In this limit of wave propagation perpendicular to the magnetic field, Ducrocq et al.³¹ discusses the existence of instability for angular frequencies in the drift frame close to cyclotron resonances, i.e., $k·Ud≈nωce$ for positive integer n. These electron Bernstein waves³² dominate the behavior of ECDI in the 2D (⁠$k⊥$⁠) limit. Additionally, it can be seen that Eq. (6), i.e., 2D dispersion relation, can be written in a form using polynomials, whereas the 3D dispersion relation requires the solution of the plasma dispersion function as shown in Eq. (4).

In 3D, i.e., for $k∥λD≠0$⁠, the plasma wave along magnetic field lines plays an important role in how the instabilities grow. Smaller wavelength oscillations that propagate parallel to the magnetic field results in a finite $k∥vth,e$⁠, and the resonant structure existing for perpendicular propagation broadens, eventually becoming a broadband ion acoustic-like spectrum.¹³ Resonance broadening occurs as the summation of Bessel functions become smaller perturbations to the ion acoustic-like terms.⁹ However, it is important to emphasize the distinction between the ion acoustic-like instability in three-dimensional ECDI and the ion acoustic instability in unmagnetized plasmas. The kinetic dispersion relation for ion acoustic waves in an unmagnetized, Maxwellian plasma for $vth,i≪ωr/k≪vth,e$ is given below:

Equation (7) can be obtained³³ from Eq. (4) when $b→∞$⁠, i.e., $B→0$⁠.

To calculate the roots of the 3D ECDI dispersion relation shown in Eq. (4), Cavalier et al. developed a numerical algorithm in Ref. 17. This algorithm relies on the structure of the ECDI dielectric function in order to rewrite the dispersion relation in the following form:

where all terms with a tilde are normalized: k is normalized by the Debye length λ_D, ω is normalized by the ion plasma frequency ω_pi, velocities (⁠$Ui+$ and U_d) are normalized by the ion acoustic speed $cs=kBTe/mi$⁠, and $M̃=mi/me$ is the ratio of the ion mass to electron mass.¹⁷ The function $g(Ω,X,Y)$ in Eq. (8) is referred to as the Gordeev function³⁴ 

Writing the 3D ECDI dispersion relation in the form of Eq. (8) enables the use of a fixed-point iteration scheme, i.e., the unstable root of Eq. (8) is obtained iteratively until convergence is achieved. To initialize the iteration, $ω0=0$ is used to evaluate the Gordeev function.¹⁷ This scheme leads to two solutions that have different values for ω_r (off by $k̃·Ũi+$⁠) and complex conjugate values for γ, i.e., the magnitude of the imaginary part of the two roots is always identical. This allows one to find the most unstable root iteratively.

The discussion in Sec. II A assumes that all ions are cold, singly charged, and unmagnetized. However, if there is some fraction of cold ions that are doubly charged, i.e., the doubly charged ion temperature is also much smaller than the electron temperature, given by $α=2ni2+/ne$ where $ni2+$ is the number density of doubly charged ions that carry two positive charges, and n_e is the electron number density; then, the contribution of the doubly charged ions to the dispersion relation must be accounted for. With two ion species, the quasineutrality condition is $n0=ne=ni++2ni2+$⁠. For the sake of consistency, the ion plasma frequency, ω_pi, is still defined using n₀, the electron number density. Hence, in the presence of a mixture of singly and doubly charged ions, the 3D kinetic dispersion relation can be written as

where χ_e is the electron contribution to the dispersion relation, identical to Eq. (2), and $Ui2+$ is the velocity of the doubly charged ion stream. Equation (10) can be non-dimensionalized with the same quantities as Eq. (8) in the following way:

where $ω̃d$ is the non-dimensionalized form of $ωd=ω−k · Ui+=ωd,r+iγ$⁠, the angular frequency shifted into the frame moving at the singly charged ion bulk velocity, and $ΔŨi$ is the non-dimensionalized form of $ΔUi=Ui2+−Ui+$⁠, i.e., the difference in the singly and doubly charged bulk velocities. It can be seen that Eq. (11) reduces to Eq. (8) when α = 0, i.e., when doubly charged ions are absent.

While the dispersion relation in the case of singly charged ions and magnetized electrons in Eq. (8) has two solutions that are complex conjugates (i.e., resulting in one unstable root), there are four possible solutions to Eq. (11) leading to two unstable roots. Therefore, the fixed-point iteration algorithm by Cavalier et al.¹⁷ struggles to find the most unstable root in the presence of two (or more) unstable roots. As a result, a new approach must be taken for equations of the form given in Eq. (10).

To numerically evaluate kinetic dispersion relations in the form of Eq. (10), the plasma dispersion function must be calculated for complex inputs. First, the method by which the plasma dispersion function is numerically evaluated is verified. The plasma dispersion function can be reformulated as $Z(ξ)=iπ exp (−ξ2)erfc(−iξ)=iπw(ξ)$⁠, where $erfc$ is the complementary error function, and $w(ξ)$ is the Faddeeva function. The Faddeeva function is numerically evaluated using two different methods, depending on the modulus of the argument ξ. For sufficiently large $|ξ|$⁠, a continued fraction expansion, such as the one reported by Gautschi,³⁵ is used. For smaller $|ξ|$⁠, the method by Zaghloul and Ali³⁶ is employed, which evaluates the Faddeeva function along a carefully selected integration contour between 0 and $z=x+iy$⁠. To verify the accuracy of the computation of $Z(ξ)$⁠, a benchmark was performed against the work by Baalrud³⁷ for both real values of $ξ∈[0,10]$ and complex values of $ξ=x+iy$ for $x∈[0,3.5]$ and $y∈[−3,3]$⁠.

For the coupled ECDI and IITSI cases, the solution of Eq. (10) can result in two or more roots that have $γ>0$⁠, corresponding to linear instability growth. Therefore, to study the coupling of ECDI and IITSI, it is necessary to use an algorithm that is capable of enumerating all roots of the dielectric function for given values of λ_D, $Ud, vth,e$⁠, ω_ce, $Ui+, Ui2+$⁠, α, and k. This is achieved in the following manner: for a given set of values for the aforementioned parameters, the dispersion relation given by Eq. (10) is evaluated for a range of values for ω_r and γ. With this information, it is possible to analyze the $ϵr(ω)=0$ and $ϵi(ω)=0$ contours for fixed k, where ϵ_r and ϵ_i are the real and imaginary components of $ϵ(ω,k)$⁠. Roots of the dielectric function exist for values of ω_r and γ which simultaneously lie on the $ϵr=0$ and $ϵi=0$ contours.

In order to locate the intersection(s) of the two contours, a standard root-finding algorithm³⁸ is considered. First, for a given condition, ϵ_r and ϵ_i are calculated on discrete ω_r and γ. Next, the line segments of $ϵr=0$ and $ϵi=0$ are obtained. The line segments are iterated over pairwise to check for points of intersection. To eliminate the unnecessary overhead of solving for the point(s) of intersection of line segments which cannot possibly cross, the pairs of endpoints of each line segment are compared to determine if an intersection is possible. If so, a matrix inversion is performed to determine the location of the intersection along the line segments. The runtime of the algorithm is affected by how finely the contours are discretized and the structure of the contours themselves, which is a function of the electric and magnetic fields, the plasma density, electron temperature, ion bulk velocities, k, and α. Acceleration of the root finding algorithm is reserved for future investigation.

For the remainder of the paper, we consider an applied electric field $E=E0x̂$ and an applied magnetic field $B=B0ẑ$⁠. Hence, $kz=k∥$⁠. This results in an $E×B$ drift that occurs in the −y direction.

Figure 1 illustrates the contours of $ϵi=0$ (red) and $ϵr=0$ (black) for three different conditions for a given k: $k̃x$ = 0.5, $k̃y$ = 0.5, and $k̃z$ = $k̃∥$ = 0.09. Here, $ΔŨi=1.57$⁠. ECDI occurs by considering $Ũd=117$⁠, and IITSI is excited assuming $α=0.5$⁠. The red and black contours contain 700 and 1000 line segments, respectively. This method of solving dispersion relations is roughly 1000 times more computationally expensive than the fixed-point algorithm of Cavalier et al.¹⁷ due to the need of root finding using discretized frequencies (both real frequency and growth rate).

Figure 1(a), i.e., the ECDI-only case, contains four closed curves, two of which comprise the $ϵi=0$ contour and two of which comprise the $ϵr=0$ contour. The two solutions of Eq. (10) with α = 0, or, equivalently Eq. (4), denoted by the blue circle in Fig. 1(a), can be seen to be the intersections of $ϵi=0$ and $ϵr=0$ contours. These solutions are complex conjugates, for which the fixed-point iteration works. In addition, it can be seen that the real frequencies are symmetric about the $ωd,r=ωr−kxUi+=0$⁠, i.e., $ωr=kxUi+$ line. Because instability growth exists for $γ>0$⁠, the solution to the ECDI dispersion relation is the root with $ωd,r>0$⁠, leading to a formation of plasma waves in the direction of $k=(kx,ky,kz)$⁠.

However, it is to be noted that the curves encompassing the $ϵi=0$ contour never intersect, and the same is true for the two curves within the $ϵr=0$ contour. Hence, the curves do not intersect at $ωd,r=0$ and γ = 0, meaning that this point is not a solution to Eq. (4). Yet, when discretizing the $(ωd,r,γ)$ space, the two contours seem to intersect and these intersections tend to be detected as the solutions to the dispersion relation with the line segment method. As a result, because these numerical intersections exist for very small $|γ|/ωpi$ values, an assumption is made in the root-finding algorithm that only roots with $|γ|/ωpi>0.005$ are considered as a solution to Eq. (4). This simplification to the root-finding algorithm works well for much of the (k_x, k_y) space, except in between low $k∥$ resonances in which the instability growth rates are considerably small. This cutoff should not affect finding the roots of larger growth rate modes which are more physically relevant.

Figure 1(b) shows the IITSI-only case, where $Ũd$ is set to zero and α is set to 0.5. Both the real and imaginary components of ϵ are modified from the ECDI-only case. The addition of the doubly charged ion contribution leads to another lobe for $ϵr=0$⁠. The pole (note: not a solution) of $ϵr=0$ and $ϵi=0$ at $ωd,r=γ=0$ corresponds to the singly charged ion contribution, while that at $ωd,r=k̃xΔŨi=0.785$ corresponds to the doubly charged ion contribution. Similar to the ECDI-only case in Fig. 1(a), the complex roots exist only when $ϵr=0$ and $ϵi=0$ contours intersect, as shown using the blue circles in Fig. 1(b). Two roots correspond to γ = 0, while the other two roots are complex conjugates. The visualization of the roots is reminiscent of the traditional two-stream instability,¹⁶ which exhibits two real roots and two complex conjugate roots if an instability condition is met. For various combinations of parameters, i.e., various wavenumbers and singly and doubly charged bulk velocities, the fastest propagating real root may also become slightly damped [$γ/ωpi∼O(−10−3)$]. The complex conjugate roots of the IITSI-only case in Fig. 1(b) have the same value of ω_r while the growth rates are same in magnitude but have opposite signs, i.e., symmetric across the γ = 0 line. The ω_r is determined mainly by the IITSI resonance, and $ωr>0$ means that the waves initiated by the IITSI will be in the direction of the k. The real part of the IITSI roots, therefore, lies between the singly and doubly charged ion velocities.¹³ 

Finally, Fig. 1(c) shows the results of a coupled ECDI and IITSI, which have features similar to the ECDI-only [Fig. 1(a)] and IITSI-only [Fig. 1(b)] cases. It can be seen in Fig. 1(c) that $ϵr=0$ contour is similar to the IITSI-only case while the two $ϵi=0$ contours (similar to the ECDI-only case) stem from both the singly and doubly charged ion branches. Such a profile yields four roots for the coupled ECDI-IITSI solution. In addition, it can be seen that the ϵ contours are no longer symmetric due to the coupling of ECDI and IITSI. Most importantly, there are two roots with $γ>0$⁠, and these roots are not complex conjugates. Hence, one needs to quantitatively evaluate both these roots to obtain the most unstable root. Depending on the shape of the ϵ contours, it can be seen from Fig. 1(c) that one unstable root is due to the $ϵi=0$ contour stemming from the singly charged ion branch and the other is due to the $ϵi=0$ contour stemming from the doubly charge ion branch. This illustrates the need of a generalized complex root-finding method for arbitrarily complex dispersion relations. This process provides a method for numerically solving dispersion relations by mapping the problem into one which conventional root-finding algorithms can be leveraged without making an assumption about the structure of the $(ωr,γ)$-solution space. The variation of the dispersion relation as a function of the wavenumber will be shown and discussed later (cf. Fig. 4).

First, the 3D ECDI dispersion relation [Eq. (8)] is solved with only singly charged ions for small $kxλD$ to verify the root-finding algorithm against the fixed-point iteration method developed by Cavalier et al.¹⁷ Twenty-one terms (⁠$n∈[−10,10]$⁠) are kept in the summation given in Eq. (9) such that the maximum truncation error is on the order of $10−6$⁠. The relevant simulation parameters are given in Table I. The angular frequency and growth rate for $kzλD=0.01$⁠, 0.03, and 0.09 are compared.

Figure 2 demonstrates excellent agreement between the fixed-point iteration scheme of Cavalier et al.¹⁷ [see Eq. (8)] and the complex root-finding algorithm as discussed in Sec. II D. In the singly charged ions only case, the discrete resonances are evident for small $kzλD$⁠, approaching the 2D limit, as shown in Eq. (6). These resonances occur for k_y near the resonance condition $kyUd≈nωce$⁠. From the values shown in Table I, $ω̃ce=51.2$ and $Ũd=117.2$⁠, corresponding to resonance conditions of $kyλD≈0.437n$ for $n∈ℕ$⁠. The second, third, and fourth modes visible in Fig. 2 are consistent with this metric, lying within $0.05ωpi$ of the predicted $E×B$ harmonic. The lowest frequency mode, which is particularly visible for $kzλD=0.01$⁠, corresponds to MTSI, as discussed in Eq. (5). Using the MTSI dispersion relation, the maximum growth rate³ is found to be $γmax=ωLH/2$⁠, where ω_LH is the lower-hybrid frequency defined by $ωLH=ωpi(1+ωpe2/ωce2)−1/2$⁠. Using the parameters mentioned above, the normalized lower-hybrid frequency is $ω̃LH=0.104$⁠, resulting in a maximum growth rate for MTSI of $γ̃max=0.052$⁠. The solution obtained from the root-finding algorithm is in good agreement with the theoretical prediction of the MTSI growth rate. Furthermore, Gary and Sanderson³⁹ predicts that in the limit of $Ud≫ωr/k$⁠, the wave number of maximum growth of ECDI is $kyλD=1/2$⁠. When the applied electric field is increased to 10⁴ V/m, the wave number of maximum growth from the 3D dispersion relation is found to be approximately $kyλD=0.68$⁠, which is in good agreement with Gary's theory. Finally, the resonant structure of ECDI observed for small $kzλD$ transitions into a broadband, ion acoustic-like instability as $kzλD$ increases. This is evident through the broadband signature of the growth rate, as shown in Fig. 2(b), over a wide range of k_y values for $kzλD≥0.03$⁠. These metrics, in addition to the benchmarking of the other modules in the 3D dispersion solver, demonstrate that the three-dimensional ECDI roots are being properly evaluated.

The plots of the drift-shifted real frequency $ωd,r=ωr−kxUi+$ and the growth rate γ as functions of k_x and k_y are illustrated in Fig. 3 using a 180 × 180 grid of k_x and k_y values. The real frequency plots are masked such that only (k_x, k_y) values that have a nonzero growth rate are shown. Note that Fig. 2 is identical to the results of Figs. 3(a)–3(c) for $kxλD=0.01$⁠. The characteristics of the three-dimensional ECDI growth rates largely reflect the same trend discussed in Fig. 2. The notable difference between the three cases is the decrease in the resonance strength as k_x increases. As such, the $E×B$ dynamics in the y-direction become less pronounced, and the instability growth rate becomes more broadband as $kzλD$ increases. Similarly, the unstable regions of IITSI, as shown in Fig. 3(d), agree with the work of Hara and Tsikata¹³ in that unstable regions exist for $kxλD<1$⁠, the growth rate increases with k_y for $kxλD∈[0.4,0.9]$⁠, and the shifted real frequency monotonically increases with k_x. The three-dimensional ECDI and IITSI growth rates for $kzλD=0.01$ are both roughly 50% off the two-dimensional $(kzλD=0)$ values evaluated by Hara and Tsikata.¹³ In contrast, the real frequency is similar between the 2D and 3D cases for ECDI and IITSI, further benchmarking the root-finding algorithm for finite k_z.

Using the root-finding algorithm, described in Sec. II D and verified in Sec. III, we now investigate the coupling between the ECDI and IITSI due to the mixture of singly and doubly charged ion streams. The results obtained in this section evaluate Eq. (11) with $α≠0$ and $Ui2+=22.71×103$ m/s in addition to the parameters shown in Table I. Using these conditions, $ΔŨi=O(1)$⁠.

In order to understand the effects of adding doubly charged ions in the presence of ECDI, Fig. 4 shows the contours of $ϵr=0$ and $ϵi=0$ as well as the roots of the 3D ECDI-IITSI dispersion relation for $kyλD=0.9$⁠, which corresponds to one of the 2D ECDI resonance modes, as can be seen from Fig. 3(a2). For all cases, there are two roots that are unstable but not complex conjugates. Here, $kzλD=0.01$ and $α=0.2$ is considered.

For a small $kxλD$ that is close to 0, as shown in Fig. 4(a), the IITSI roots are close to $(γ,ωd,r)=(0,0)$ as the origin of the doubly charged ion branch is $ω̃d,r=k̃xΔŨi+≈0$⁠. The other root, i.e., the ECDI root, is largely unperturbed, similar to the ECDI-only case, as shown in Fig. 1(a). In Fig. 4(b), when $kxλD=0.2$⁠, the singly and doubly charged ion branches separate due to the increasing k_x, and the growth rate for the IITSI root starts to increase, similar to the IITSI-only case as shown in Fig. 1(b). However, the unstable root due to the ECDI mode is still dominant over the IITSI root. The presence of the IITSI lobes increases the instability growth rate of the rightmost root in Fig. 4(b) by about 0.2% compared to ECDI-only growth rate with the same k. Additionally, the rightmost ECDI root in Fig. 4(b) gets pushed toward larger real frequencies by 32% compared to the ECDI-only phase case with the same k.

For $kxλD=0.4$⁠, as shown in Fig. 4(c), the $ϵi=0$ contour transitions from a combination of the ECDI and IITSI features into two separate ECDI type contours stemming from the singly and doubly charged ion branches. It can be seen in Fig. 4(c) that the two roots in the middle are no longer a complex conjugate due to the coupling of the ECDI and IITSI. In addition, here, the unstable root in the middle (at $ω̃d,r≈0.4$⁠) shows a larger growth rate compared to the rightmost root (at $ω̃d,r≈1$⁠). The growth rate of the IITSI branch becomes dominant over the ECDI root. This leads to a continuous change in γ but a discontinuous change in ω_d,r as k_x increases. As $kxλD$ continues to increase through 0.8 and 1.0, as shown in Figs. 4(d) and 4(e), the shape of the contours begins to further change. Due to the separation of the singly and doubly charged ion branches as the k_x increases, two separate $ϵi=0$ contours that are similar to the ones observed for the ECDI case [as shown in Fig. 1(a)] appear around the singly and doubly charge ion branches, while the $ϵr=0$ contours still remain similar to the IITSI case [as shown in Fig. 1(b)]. As the $ϵi=0$ contours separate, the two middle roots are not complex conjugates (e.g., the $ωd,r$ separates and the magnitudes of γ are no longer identical). The unstable root at $0<ωd,r<ωpi$ has a larger growth rate than the rightmost unstable root at $ωd,r>ωpi$⁠. As $kxλD$ further increases, the central loop of $ϵr=0$ stretches until it eventually breaks and forms two ECDI-type contours, as shown in Fig. 4(f), i.e., $kxλD=1.1$⁠. Once the $ϵr=0$ contours finally separate, the instability has transitioned to two separate ECDI contours. Each pair of closed contours in Fig. 4(f) is identical to the ECDI-only contour in Fig. 1(a). It can be considered in Fig. 4(f) that one closed contour is for singly charged ions, as it is centered at $ωd,r=0$⁠, i.e., $ωr=kxUi+$⁠, and the other closed contour is for doubly charged ions, as it is centered at $ωd,r=kxΔUi+$⁠, i.e., $ωr=kxUi2+$⁠.

In summary, the most unstable root due to the coupling of ECDI and IITSI can be characterized by analyzing the ϵ_r and ϵ_i contours as follows:

• At lower k_x, the growth rate of the middle (IITSI) root is smaller than that of the rightmost (ECDI) root. See Figs. 4(a) and 4(b).

• For intermediate k_x, the growth rate of the middle root grows and becomes larger than that of the rightmost root. The $ϵi=0$ contours separate into ECDI-like features stemming from the singly and doubly charged ion branches, while $ϵr=0$ remain similar to the IITSI feature, illustrating the coupling between ECDI and IITSI. See Figs. 4(c)–4(e).

• For larger k_x, two separate contours that correspond to the singly and doubly charged ions are observed. Hence, the dispersion relation is a superposition of two ECDI contours as the doubly charged ion branch $ωd,r=kxΔUi+$ separate from the singly charged ion branch $ωd,r=0$⁠. See Fig. 4(f).

The complex root-finding algorithm discussed above is used to solve the coupled ECDI-IITSI dispersion relation for two values of α: 0.05 and 0.2. A doubly charged ion fraction up to 0.2 is found in HETs.⁴⁰ Similar to Fig. 3, the growth rates and real frequencies for $α=0.05$ and $α=0.2$ cases are plotted as functions of $kxλD$ and $kyλD$ discretized over 180 values from 0 to 2, respectively. Here, three different $kzλD$ values, $kzλD=$ 0.01, 0.03, and 0.09, are shown, similar to the values used in Figs. 2 and 3(a)–3(c).

Figure 5 illustrates the ECDI and IITSI coupling for $α=0.05$⁠, which corresponds to a weak IITSI.¹³ In contrast to the singly charged ECDI-only case for $kzλD=0.01$ in Fig. 3(a) where the ECDI resonances, i.e., $kyUd≈nωce$⁠, are sufficiently separated, Figs. 5(a2), 5(b2), and 5(c2) show that the most unstable roots are more complex, having influence from both the ECDI and IITSI. The colorbar spans the same range of values for real frequencies, as shown in Figs. 5(a1)–5(c1), and growth rates, as shown in Figs. 5(a2)–5(c2), respectively.

It can be seen from the growth rate plots in Figs. 5(a2)–5(c2) that the IITSI mode stays similar for different values of $kzλD$⁠. This is because the IITSI is an electrostatic instability caused by the streams of singly and doubly charged ions in the axial direction. On the other hand, the ECDI modes transition from resonance-type solution to broadband structures due to the increasing $kzλD$⁠, as shown in Figs. 3(a2)–3(c2). The coexistence of the IITSI and ECDI modes results in a complex structure of the most dominant unstable solution. In the limit of the 2D case¹³ (i.e., $kzλD→0$⁠), it was found that the most unstable mode of the ECDI-IITSI solution is at $kxλD=0$⁠, as the 2D ECDI growth rate due to the resonant modes [see Eq. (6)] is larger than the IITSI growth rate. This means that the plasma wave generated due to the coupled ECDI-IITSI occurs in the direction of the $E×B$ drift in the 2D limit. Interestingly, in the 3D cases, shown in Figs. 5(a2)–5(c2), as the 3D ECDI growth rate is smaller than the 2D ECDI growth rate due to the broadband structure as $kzλD$ increases, the maximum growth rate of the 3D ECDI-IITSI dispersion can be found to exist in the region where the IITSI occurs. The maximum growth rate is now at a finite $kxλD$ and a finite $kyλD$⁠, indicating that the plasma wave generated due to the coupled ECDI-IITSI occurs in an oblique direction (in the x–y plane) for the 3D cases.

Figures 5(a1)–5(c1) show the real frequencies of the maximum growth rate, exhibiting more clearly the transition between ECDI and IITSI modes. At smaller $kxλD$ [e.g., $kxλD≤0.3$ for Figs. 5(a) and $kxλD≤0.15$ for Figs. 5(b) and 5(c)], the maximum growth rate is large where the ECDI mode is most unstable. At a larger $kxλD$ (e.g., $kxλD≥0.3$⁠), the ECDI mode becomes less dominant, and the IITSI mode results in the maximum growth rate. As discussed in Figs. 4(b) and 4(c), the real frequency, ω_r, of the ECDI root is larger than that of the IITSI root, which explains the discontinuity in the real frequency of the most unstable mode, as shown in Figs. 5(a1)–5(c1). At a much larger $kxλD$ (e.g., $kxλD≥0.8$⁠), the real frequency and the growth rate do not exhibit any discontinuity as the IITSI root of the singly and doubly charge ions, as shown in Figs. 4(c)–4(e), transitions to the ECDI root of the singly charged ion contribution, as shown in Fig. 4(f). The switching between low and high real frequency roots exemplifies the sensitivity of this coupled instability to the system dimensions, as this determine the largest allowable wavelength modes (and as such the smallest wavenumbers) in each direction. Depending on the smallest allowable wavenumber, the coupled instability can be characterized by different values for ω_r.

Figure 6 shows the real frequency and the growth rate of the most unstable solution for $α=0.2$⁠, which corresponds to a strong IITSI.¹³ While the profile of the most unstable roots is similar to the $α=0.05$ case shown in Fig. 5, it can be seen from Fig. 6 that the IITSI roots become more dominant compared to the ECDI roots. The increase in the IITSI growth rate is consistent with Ref. 13. The maximum growth rate is 33% larger for the coupled $α=0.2$ case (Fig. 6) compared to the coupled $α=0.05$ case (Fig. 5). The maximum growth rate of the 3D ECDI-IITSI case with $α=0.2$ and $kzλD=0.01$ [Fig. 6(a2)] is roughly 7% larger than the IITSI-only growth rates [Fig. 3(d)]. This maximum occurs in the region near $kxλD=0.7$ and $kyλD=1.8$ where the ECDI resonance overlaps with the IITSI region of greatest instability growth. As $kzλD$ increases, the maximum coupled ECDI-IITSI growth rate approaches that of IITSI-only.

One interesting observation can be made in Fig. 6(a1). The singly charged ECDI resonance dominates the growth rate of the doubly charged ECDI resonance for many of the (k_x, k_y) values. However, at $kxλD>1$ and for large enough α, i.e., $α≥0.2$⁠, the growth rate of the doubly charged resonance is able to surpass the growth rate of the singly charged resonance. In contrast to Fig. 5(a1), the drift-shifted frequency of the fastest growing mode in the ECDI resonances discontinuously changes as the ECDI root switches from the singly charged ion branch to the doubly charged ion branch. This can only occur for sufficiently large values of α, as the relative importance of the doubly charged ion term in the ion susceptibility is a function of $kxΔUi$ and $2α/(1−α)$⁠. For $α=0.2$⁠, the quantity $2α/(1−α)$ is 0.5 compared to the $α=0.05$ case, in which the quantity is about 0.1, justifying why this effect is only seen in the $α=0.2$ case. The drift-shifted frequency of the doubly charged resonance, and hence the phase velocity of the wave, is significantly larger than the singly charged resonance. For instance, this can be seen from Fig. 4(f), where there are two unstable roots: one with $ωd,r/ωpi≈0.8$ and the other with $ωd,r/ωpi≈2.2$⁠. In 3D, the resonances are thermally broadened, which can cause the singly charged ECDI root to have a smaller growth rate, resulting in a larger growth rate due to the doubly charged ECDI root.

One of the greatest advantages of considering the 3D ECDI and IITSI dispersion relation is the ability to capture the thermal broadening of the cyclotron resonances as $kzλD$ increases.¹⁷ In 1D and 2D studies, the ECDI roots exhibit resonance conditions, as can be seen in Eq. (6). As a result, the growth rates of the resonances in 2D are larger than the 3D case for a given electric field, magnetic field, electron temperature, and electron density. In the 2D ECDI-IITSI dispersion, i.e., $k∥=0$⁠,¹³ the ECDI growth rates dominate the IITSI growth rates for all values of $kxλD$⁠, i.e., the largest growth rates occur at $kxλD=0$⁠. This means the plasma wave initiated by 2D ECDI-IITSI traverses in the azimuthal direction with the allowable k_y values within the plasma device of interest, e.g., Hall effect thruster discharge channel. However, in the 3D ECDI-IITSI case (i.e., $kz≠0$⁠), the maximum growth rate of the ECDI resonances becomes lower and widened, which resembles a transition from the resonant modes to broadband ion-acoustic modes, as k_z increases. Consequently, the IITSI growth rates can come to dominate over the ECDI growth rates for intermediate $kxλD$⁠. For $kzλD=0.01$⁠, the maximum growth rate occurs in the IITSI region along the ECDI resonances, as the ECDI and IITSI couple. The first three ECDI resonances at $kyλD≈0.5,0.9$⁠, and $1.4$ for $α=0.05$ have the largest growth rates with wavenumbers at $kxλD≈0.64$⁠, corresponding to waves propagating at 48°, 34°, and 27° with respect to the axial electric field, respectively. Thus, if ECDI and IITSI are present in the same region within HETs, the results suggest that the resulting waves driven by the coupled ECDI-IITSI propagate at some oblique angle to the axial electric field, rather than purely in the azimuthal direction as previously predicted by two-dimensional theory. It is important to note that finite ion temperature effects could serve to stabilize the coupled instability,²⁸ which is reserved for future work. Furthermore, because this prediction is made with linear dispersion theory, the exact values of the propagation angle with respect to the axial direction may be altered due to nonlinear effects. Nevertheless, considering the 3D kinetic dispersion relation for coupled ECDI and IITSI may explain the experimental evidence of axially propagating modes observed by Tsikata et al.¹² 

The analysis of dispersion relations using complex root finding enables the study of plasma kinetic instabilities that have multiple unstable roots. The algorithm for solving dispersion relations discussed in this paper can be extended to any arbitrary dispersion relation that has more than one unstable root to evaluate the most unstable mode. As long as the dispersion relation can be numerically evaluated for various $ω=ωr+iγ$ values, the intersections of the contours of $Re(ϵ(ω,k))=0$ and $Im(ϵ(ω,k))=0$ can be found, yielding valid solutions to the dispersion relation. While not as computationally efficient as the fixed-point algorithm of Cavalier et al.,¹⁷ this complex root-finding algorithm is a robust method for analyzing complicated dispersion relations.

Applying the complex root-finding algorithm to the three-dimensional coupling of ECDI and IITSI reveals interesting behavior in ω_r and γ as the wavenumber and direction of propagation changes, including the switching between ECDI and IITSI roots for intermediate values of wavenumber parallel to the electric field. The largest difference between the two-dimensional and three-dimensional study of the coupling of ECDI and IITSI is the location of the largest growth rates. In previous two-dimensional studies,¹³ the ECDI growth rate at k_x = 0 always dominates the IITSI growth rates at finite k_x. However, the present three-dimensional theory indicates that the highest growth rate for a given k_y occurs at finite k_x, resulting in a wave that propagates at an oblique angle to the electric field. This finding gives new insight into possible waves present in the plume regions of Hall effect thrusters in which both ECDI and IITSI are present.

The authors acknowledge the discussions with S. Tsikata, J. C. Adam, and A. Héron. This work was supported by a NASA Space Technology Graduate Research Opportunity, under Grant No. 80NSSC21K1273, NASA through the Joint Advanced Propulsion Institute, a NASA Space Technology Research Institute under Grant No. 80NSSC21K1118, and the U.S. Department of Energy, Office of Science, Office of Fusion. Energy Sciences, under Award No. DE-SC0020623.

The authors have no conflicts to disclose.

Andrew Christopher Denig: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Methodology (lead); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal). Kentaro Hara: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (lead); Methodology (supporting); Supervision (lead); Visualization (equal); Writing – original draft (supporting); Writing – review & editing (equal).

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