© ESO 2018

1. Introduction

Diluted ion beams propagating along the background magnetic field B₀ are widespread in space and astrophysical plasmas, including solar wind (Marsch 2006, and references therein), terrestrial foreshocks (Paschmann et al. 1981, and references therein), supernova remnants (Bell 2005, and references therein), and many other astrophysical environments (Zweibel & Everett 2010, and references therein). As the plasmas are typically quasi-neutral, the background electrons tend to follow the beam ions compensating their current. Depending on particular settings, the compensating currents can also be provided by other plasma components, like co-streaming electron beams injected simultaneously with the ion beams. Plasma instabilities developing in such compensated-current systems not only regulate the plasma and beam parameters keeping them close to the marginally unstable states, but can also be important sources for the background plasma heating, energetic particle acceleration, and amplification of the background magnetic field.

Plasma waves in the compensated-current systems can be driven unstable by resonant (Duijveman et al. 1981; Gary 1985; Vojtenko et al. 1990) and non-resonant (Winske & Leroy 1984; Bell 2004; Achterberg 2013) wave–particle interactions. Resonant kinetic instabilities of various wave modes, driven by the beam ions, have been studied extensively in the past. Parallel-propagating Alfv én and fast waves have been found to be strongly unstable for beam velocities higher than a few Alfvén velocities (e.g., Gary 2005; Marsch 2006, and references therein). Concurrent instabilities of oblique (kinetic) Alfvén waves come into play at lower (but still super-Alfvénic) beam velocities (Voitenko 1998).

The above-mentioned instabilities can be driven directly by the beam ions (Sentman et al. 1981; Winske & Leroy 1984; Gary 1985) or by the electron return currents (Winske & Leroy 1984; Bell 2004; Chen & Wu 2012, and references therein). The non-compensated electron currents flowing along B₀, may also drive both the resonant (Voitenko 1995) and non-resonant (Malovichko & Iukhimuk 1992; Malovichko 2007) instabilities of Alfvén waves. The simplest case of purely parallel propagating Alfvén waves has been considered in application to the current-carrying coronal loops (Malovichko & Iukhimuk 1992), where these waves appeared to always be unstable. Later on, the analysis was extended by accounting for the oblique propagation (Malovichko 2007) and the currents curried by low-density beams (Malovichko 2010), and applied to the terrestrial magnetosphere and coronal loops.

Self-consistent modifications of the background magnetic field by the electric currents, neglected by Malovichko & Iukhimuk (1992), Malovichko (2007); and Voitenko (1995), may reduce or even stabilize current instabilities. This issue does not concern instabilities developing in compensated-current systems. Such systems, formed around supernova remnants by high-energy streaming cosmic rays (CRs), were studied by Bell (2004), who found a new non-resonant Alfvénic instability (hereafter Bell instability). Since then, the Bell instability and its modifications have attracted a lot of interest (see, e.g., Amato & Blasi 2009; Bret 2009; Zweibel & Everett 2010; Schure et al. 2012; Achterberg 2013; Kobzar et al. 2017, and references therein). Following Bell (2004), the primary focus was on the unstable modes with finite k_zV̄_bz/ω_Bi propagating along B₀ (where V̄_bz is a characteristic velocity of the beam ions along the mean magnetic field B₀ || z; k_z is the parallel wavenumber; and ω_Bi is the ion-cyclotron frequency). Compensated currents can also drive an oblique Alfvén instability (Malovichko et al. 2014) for which the perpendicular wave dispersion due to finite k_⊥V_Tb⊥/ω_Bi is essential (k_⊥ and V_Tb⊥ are the perpendicular wavenumber and beam thermal velocity in the plane ⊥ B₀).

Other electrostatic and electromagnetic instabilities may develop in compensated-current systems (see, e.g., Gary 2005; Bret 2009; Brown et al. 2013; Marcowith et al. 2016, and references therein). Which wave modes grow fastest critically depends on the beam and plasma parameters. In the case of cold diluted proton beams propagating along B₀, the electrostatic two-stream and Buneman instabilities are much faster than the electromagnetic Alfvénic instabilities (see, e.g., Fig. 44 in Bret et al. 2010). Nevertheless, as is noted by Bret et al. (2010), these electrostatic instabilities are quickly saturated, and then electromagnetic Alfvénic/Bell instabilities come into play. In the hot beam/plasma systems, where the two-stream/Buneman instabilities cannot develop, the electromagnetic Alfvénic/Bell instabilities dominate.

The Bell instability has a maximum growth rate γ_Bell ≃ 0.5j̄_bω_Bi, where j̄_b = n_bV_b/ (n₀V_A) is the beam current normalized by the Alfvén current. This maximum is attained at the parallel wavenumber |k_zm| V_A/ω_Bi = 0.5j̄_b and the perpendicular wavenumber k_⊥ = 0. These expressions are exactly the same as for the instability studied earlier by Winske & Leroy (1984) in application to the terrestrial foreshock. The difference is that the role of V̄_bz in the setting considered by Winske & Leroy (1984) is played by the bulk velocity of the beam V_b rather than the large velocity spread of CRs. Both the Winske–Leroy and Bell instabilities grow fastest when the wave vector k is parallel to B₀; they are physically the same instability that can be named the compensated-current parallel instability (CCPI).

The physical mechanism of the CCPI is related to the fact that for sufficiently small parallel wavelengths and sufficiently high V̄_bz, the beam protons become partially demagnetized (unfrozen off the perturbed magnetic field). The demagnetization reduces the beam contribution to the fluctuating currents δj ⊥ B₀ flowing along the twisted perturbed magnetic field lines, whereas the electron currents remain magnetized, thus providing the non-compensated fluctuating transversal currents. These currents amplify the initial perturbations via the positive feedback loop giving rise to CCPI. This kind of instability is sometimes called reactive.

Surprisingly, despite its importance in astrophysical applications, the CCPI theory is still poorly developed. Many important properties of the instability (the wavenumber dependence of the instability growth rate, behavior of the maximum growth rate in the parameter space, instability boundaries in the parameter spaces, etc.) have not been fully investigated. In the present paper, we study CCPI of Alfvén waves in more detail in the framework of kinetic theory. We consider a simple model of the compensated-current system where the hydrogen plasma is hydrogen plasma is penetrated by the low-density proton beam and the beam current and charge are compensated by the background electrons. Despite its simplicity, this model is applicable to the reactive CCPI driven by compensated currents in many space and astrophysical environments.

2. Plasma model and dispersion equation for Alfvén waves

We consider a three-component plasma consisting of the background steady ion component (i), the low-density ion beam (b) propagating with velocity V_b along z || B₀, and the electron component (e) providing the neutralizing current and charge:(1) (2)

We assume here that the beam ions (b) and the background ions (i) are protons. All plasma components are modeled by the shifted Maxwellian velocity distributions(3)

where n_s, V_s, , T_s, and m_s are the mean number density, parallel bulk velocity, thermal velocity, temperature, and particle mass of the plasma species s, and v = (v_x, v_y, v_z) – velocity-space coordinates. The subscripts z and ⊥ indicate directions parallel and perpendicular to B₀. The plasma model defined by Eqs. ((1)–(3)) has been extensively used in the past (see, e.g., Gary 2005, and references therein). The neutralizing current can also be provided by the co-propagating electron beam (see, e.g., Zweibel & Everett 2010 and references therein); however it does not alter the reactive CCPI for low-density ion beams n_b ≪ n₀ (Amato & Blasi 2009).

The non-trivial solutions to the Maxwell–Vlasov set of equations exist if the wave frequency ω and the wave vector k = (k_x, k_y, k_z) satisfy the following dispersion equation (see, e.g., Alexandrov et al. 1984)(4)

where ε_ij is the dielectric tensor and δ_ij is the Kronecker delta. For the parallel-propagating modes with k_x = k_y = 0, the components of the dielectric tensor given by Alexandrov et al. (1984) reduce to(5)

where ξ_s,n = (ω − k_zV_s + nω_Bs) / (k_zV_Ts), ω_Ps (ω_Bs) is the plasma (cyclotron) frequency. Instead of the plasma dispersion function W (x), we use the function(6)

introduced by Alexandrov et al. (1984). It has the following asymptotic expansions:(7)

and(8)

where η = 0 for Im x > 0, η = 1 for Im x = 0, and η = 2 for Im x < 0.

In the case of parallel propagation, the dispersion Eq. (4) splits into two independent equations,(9)

describing left-hand (sign −) and right-hand (sign +) polarized electromagnetic waves. In what follows we consider the left-hand polarized Alfvén branch undergoing the compensated-current instability. Taking into account quasi-neutrality (2) and current compensation (1), Eq. (9) for Alfvén waves can be written as(10)

where

In the following sections we consider important limits of (10) typical for the reactive CCPI instability.

3. Dispersion relation for parallel-propagating waves

As we are going to analyze the reactive non-resonant instability, we neglect the contribution of the imaginary part of J₊ (ξ_b,−1). Furthermore, we consider a low-frequency instability with |ω/ω_Bi| smaller than other terms in ξ_b,−1, which allows us to neglect the ω-dependent part in the argument of function J₊. In this case (10) reduces to the following quadratic equation with respect to ω/ω_Bi:(11)

where(12)

ζ_b = −V_b/V_Tb − 1/ (k_zρ_Tb), and ρ_Tb = V_Tb/ω_Bi. To avoid misunderstanding, we note that although ρ_Tb looks like the ion beam gyroradius, it is defined by the parallel beam temperature rather than the perpendicular one, and here it has a different physical meaning.

Equation (11) is the second-order eigenmode equation for Alfvén waves modified by the ion beam and return electron current (second and fourth terms, respectively). Its solution is straightforward:(13)

From (13) it is obvious that the instability can be driven by the last term under the square root when k_zV_b > 0. In what follows we assume V_b > 0 considering potentially unstable waves with k_z > 0 (in the case of V_b < 0, the identical instability develops for k_z < 0 ). In the absence of the beam, Eq. (13) reduces to the Alfvén wave dispersion, ω = k_zV_A at n_b = 0.

The wave with dispersion (13) becomes unstable when the last term under the square root dominates. This term represents effects due to the electron current. The growth rate γ = Im[ω] of the corresponding instability is(14)

Here we introduce the cumulative destabilizing parameter(15)

that includes all beam parameters. One can think of it as of product of the normalized current j̄_b = n_bV_b/ (n₀V_A) and velocity spread V̄_Tb = V_Tb/V_A of the beam ions. The growth rate (14) is analyzed below analytically and numerically, and its scalings are found in some important limits. It is interesting to note that the right-hand polarized magnetosonic instability can be obtained from the above equation by changing the sign of the first term under the square root (the magnetosonic instability hence requires k_zV_b < 0).

4. Compensated-current instability driven by hot ion beams

By hot beams we mean the beams whose thermal velocity is significantly larger than the bulk velocity, V_Tb ≫ V_b. For such beams, the growth rate (14) can be simplified by neglecting the small term V_b/V_Tb in ξ_b,−1. The argument of J₊ is then simplified to ξ_b,−1 ≈ −1/ (k_zρ_Tb) ≡ ζ_b. In this case γ_k depends on the normalized parallel wavenumber k_zρ_Tb and two dimensionless bulk parameters, V_A/V_b and α_b. Then the (maximum) instability growth rate γ_m = max_kγ appears to be a function of α_b and V_A/V_b only, whereas the dependence on the general multiplier V_A/V_Tb is trivial and can be excluded by the renormalization of γ_m. We note that the hot-beam condition V_Tb > V_b restricts the applicability range of the analytical results obtained below, but in general it does not restrict the instability range (see also Sect. 6).

4.1. Instability areas in the parameter space

Here we find the instability threshold and the instability area in the parameter space (α_b, V_A/V_b). To this end, we present the growth rate (14) in the following useful form:(16)

From (16), the instability condition is obtained as(17)

Since the right-hand side of (17) is positive, it is obvious that only super-Alfvén beams, V_b > V_A, can trigger instability. Therefore, the absolute threshold for the beam velocity is and the system is stable with respect to reactive CCPI for all V_b < V_A.

Using (17), it is also possible to find the threshold for _b analytically. First, solving (17) with respect to the k_k-dependent term A_k/ (k_zρ_Tb), we find that the unstable wavenumbers k_z should satisfy(18)

When the velocity threshold is exceeded, V_b > V_A, the right boundary of (18) is always larger than the left boundary making the interval between them non-empty. As the function A_k/ (k_zρ_Tb) is limited by the maximum value max_k [A_k/ (k_zρ_Tb)] ≈ 0.411 achieved at , the condition (18) can only be satisfied for sufficiently large α_b. From the left-hand inequality, it immediately follows the instability condition for α_b and the corresponding threshold:(19)

The instability condition is satisfied above the threshold curve defined by (19), which is shown in Fig. 1 by the solid line. The unstable area above this curve in the parameter space (α_b, V_A/V_b) is shaded. The dependence of the threshold on V_A/V_b is rather weak; it grows from the minimum value at V_A/V_b → 0 to the maximum value at V_A/V_b → 1. The absolute threshold for α_b below which the system is stable is . The meaning of the right boundary in (18), shown in Fig. 1 by the dashed line, is clarified in the following subsection.

Fig. 1. Instability threshold in the parameter space (αb, VA/Vb) (solid line); the CCPI develops at all . The dashed line shows the split threshold above which there are two separate ranges of unstable wavenumbers kz

In terms of the normalized beam current j̄_b = n_bV_b/ (n₀V_A) and thermal velocity V̄_Tb = V_Tb/V_A, (19) can be written as . Then the instability condition for the beam thermal velocity reads as(20)

This threshold-like condition is an important new result quantitatively demonstrating the destabilizing effect of the beam temperature. It shows the threshold above which the velocity spread of the beam ions triggers the instability even for weak beams.

Similarly, the threshold condition for the beam current can be written as(21)

which quantifies the range of unstable beam currents. Again, it is seen that even a weak ion beam can activate CCPI if the beam thermal velocity is sufficiently high. In particular, the beam current required for the instability can be many orders of magnitude smaller than the Alfvén current.

We note that varies slowly for fast super-Alfv énic beams and can be approximated as at V_b/V_A > 3. For rough estimations, in all velocity ranges can be replaced by its average value 3.5.

4.2. Unstable wavenumber ranges

Properties of CCPI are illustrated further by Figs. 2 and 3 showing all three terms of the condition (18): the left and right boundaries, and the function (k_zρ_Tb)⁻¹A_k. The unstable ranges where (18) is satisfied are shaded. A regular single-peak behavior of the function (k_zρ_Tb)⁻¹A_k, as seen in Figs. 2 and 3, allows us to investigate how the unstable wavenumber range evolves with α_b.

Fig. 2. Illustration of the condition (18) for VA/Vb = 0.9 and αb = 6. In this case and there is only one unstable wavenumber range (shaded area).

Fig. 3. Illustration of the condition (18) for VA/Vb = 0:9 and αb = 10. In this case and there two unstable wavenumber ranges presented by two shaded areas.

When α_b increases but is still smaller than , the left boundary of (18) decreases remaining above the maximum of (k_zρ_Tb)⁻¹A_k. In this case there are no unstable wavenumbers and the system is stable. Once α_b rises above , the decreasing left boundary of (18 ) drops below the maximum of (k_zρ_Tb)⁻¹A_k and the unstable wavenumber range k_z1 < k_z < k_z2 appears, where k_z1 and k_z2 are the lower and upper roots of the equation(22)

As long as α_b is not far from the threshold , there is a single unstable wavenumber interval surrounding . This situation is illustrated in Fig. 2, where V_A/V_b = 0.9, , and . However, when α_b increases further, the right boundary of (18) also drops below the maximum of (k_zρ_Tb)⁻¹A_k, which happens at(23)

In this case, shown in Fig. 3 for α_b = 9, the right-hand side inequality of (18) is not satisfied in the range , where and are the lower and upper roots of the equation(24)

Instead of unstable, we have now a prohibited wavenumber range around . As a result, the unstable wavenumber range splits into two: the first unstable range is and the second .

The split threshold (23) is shown in Fig. 1 by the dashed line. For parameter values above this line, the instability develops in two wavenumber ranges, as mentioned above. These unstable ranges are shown in Fig. 3 by the shaded areas.

Furthermore, Fig. 4 shows the α_b dependence of the unstable wavenumber ranges, where the outer and inner boundaries are respectively defined by the left-hand and right-hand margins of (18). It is seen that below there is no instability, at there is a single unstable range of k _z, and above there are two unstable ranges.

Fig. 4. Unstable wavenumber ranges in the (αb, kz) plane for VA/Vb = 0.9. The outer boundary is defined by the left-hand side and the inner boundary by the right-hand side of condition (18). It is seen that below there is no instability, at there is a single unstable range of kz, and above there are two unstable ranges.

From Fig. 3 it is obvious that k_z1ρ_Tb and are located between k_zρ_Tb ≈ 0.77, where (k_zρ_Tb)⁻¹A_k is zero, and , where (k_zρ_Tb)⁻¹A_k is maximum. This wavenumber range corresponds to −1.3 < ζ_b < −0.65, where ReJ₊ (ζ_b) can be approximated by the liner numerical fit(25)

Using this in (22) and (24), we find k_z1 and as(26) (27)

From these expressions we see that with increasing α_b the difference between k_z1 and decreases, , and the first unstable range becomes very narrow.

On the other hand, the roots k_z2ρ_Tb and bounding the second unstable range, are located above , where ζ_b > −0.65. Then, using the small argument expansion (7) for ReJ₊ (ζ_b), we find(28) (29)

At large values of α_b, both the width of the second unstable range and the gap between the unstable ranges grow linearly with α_b.

In summary, the most important analytical result obtained here is the instability boundary in the parameter space (V_A/V_b; α_b), which can be used directly to analyze observational data. The compensated-current systems with V_A/V_b < 1 and are unstable. The unstable area in the parameter space (V_A/V_b; α_b) is divided further by into two unstable sub-areas: with one unstable wavenumber range, and with two unstable wavenumber ranges.

4.3. Instability growth rate

Once α_b rises above , an unstable range between k_z1 and k_z2 appears. The instability growth rate (14) as function of k_z is shown in Fig. 5. The plasma parameters α_b and V_A/V_b in this figure are chosen in such a way as to illustrate the behavior of CCPI in the unstable wavenumber ranges found above. So, the case α_b = 6 with one unstable wavenumber range is shown by the dashed line and the case α_b = 10 with two unstable wavenumber ranges is shown by the solid lines. The dotted curve in Fig. 5 is for the case α_b = 8, which is close to the splitting threshold. It is seen that when the right instability boundary in (18) approaches the maximum of function (k_zρ_Tb)⁻¹A_k, the valley and the second peak in γ_k appear. This happens at , where(30)

Fig. 5. Wavenumber dependence of the instability growth rate driven by super-Alfvénic ion beams with VA/Vb = 0.9 for three values of αb: αb = 6, 8, and 10. For larger αb, the unstable area and the maximum growth rate extend to larger kzρTb.

is the value of α_b at which a local “plateau” in γ_k occurs at the wavenumber where ∂γ_k/∂k_z = 0 and . For all , the secondary peak of γ_k exists at . Since , the secondary peak arises before the interval of prohibited wavenumbers appears.

It is seen that CCPI is stronger and the most unstable wavenumbers are larger for larger α_b. The secondary peak that appears at is lower than the main peak. These trends are confirmed below analytically.

The most unstable wavenumber and the corresponding maximum growth rate γ_max can be found by maximizing (14) with respect to k_z, γ_max = max_k (γ_k), which we call the CCPI growth rate. The normalized CCPI growth rate γ_max/ω_Bi as a function of n_b/n₀ and V_b/V_A is shown in Fig. 6 for hot beam with V_Tb/V_A = 10². It is seen that γ_max increases rapidly, roughly proportional to both n_b/n₀ and V_b/V_A, which means it is proportional to the current n_bV_b. This behavior agrees with the current nature of CCPI confirmed below analytically by (34) and (35).

Fig. 6. Normalized growth rate γmax/ωBi as a function of nb/n0 and Vb/VA for hot beam with VTb/VA = 102. γmax increases regularly with both nb/n0 and Vb/VA once the threshold is exceeded.

The threshold for n_b/n₀ (V_b/V_A) is lower for smaller V_b/V_A (n_b/n₀), in agreement with (19). In particular, the velocity threshold decreases with n_b/n₀ and reaches the minimum value when n_bV_Tb/ (n₀V_A) → 1.

Dependence of γ_max on the thermal velocity V_Tb is somewhat different (see Fig. 7). First, near the threshold, γ_max increases very rapidly with V_Tb. But this increase quickly slows down as V_Tb departs further from the threshold. Already at , γ_max becomes virtually independent of V_Tb.

Fig. 7. Normalized growth rate γmax/ωBi as a function of VTb/VA for nbVb/ (n0VA) = 0.05 (solid curve), 0.1 (dotted curve), and 0.15 (dashed curve). Starting from zero, the growth rate γmax increases rapidly with VTb, but this increase is quickly saturated. Larger currents nbVb/ (n0VA) result in larger γmax for all VTb.

To understand this behavior, we proceed with the analytical analysis. Here we take into account that in the wavenumber range , where the growth rate attains its maximum, the lowζ_b approximation ReJ₊ (ζ_b) ≈ −(k_zρ_Tb)⁻² is valid. Thus, using A_k = 1 + ReJ₊ (ζ_b) ≈ 1 − (k_zρ_Tb)⁻² in (16), we find the following approximation for the maximum of γ_k:(31)

This maximum occurs at(32)

The last term in the square parentheses in (31) is adjusted by replacing the approximate numerical value by to make it compatible with the exact α_b threshold (19). We verified numerically that the approximation (31) is good for arbitrary α_b, both near the threshold and far from it. In general, with the larger beam velocity and/or temperature, the smaller beam density is needed for instability.

The explicit dependence of the instability growth rate on the beam thermal velocity V̄_Tb follows from (31):(33)

It is seen that γ_max increases quickly with V̄_Tb once the threshold is overcome, . The fast increase of γ_max reflects the instability response to the progressive demagnetization of the beam ions as their velocity spread increases above the threshold.

However, when V̄_Tb becomes large enough, , the term containing it becomes negligibly small and γ_max becomes virtually independent of V̄_Tb. In this high-temperature regime the beam ions are fully demagnetized and the further increase of V̄_Tb no longer affects the instability. This regime corresponds to the asymptotic over-threshold limit , where γ_max simplifies to(34)

The familiar threshold velocity of the beam, , is still present in (34), but the temperature dependence is already missed, as can be observed in Fig. 7 at large V_Tb.

The maximum growth rate (34) simplifies further for the fast beams with V_b/V_A > 3,(35)

with the most unstable parallel wavenumber . The asymptotic scaling (35) recovers the scaling obtained by Bell (2004). As is seen from Fig. 7, expressions (34) and (35) provide good estimations for γ _max at , which also quantifies the meaning of “asymptotic regime” in terms of V_Tb. It appears that the expressions found by Bell are only valid in this asymptotic regime.

For , the secondary peak arises at k_zρ_Tb < 1.54, where we can use approximation (25). Then for this peak we obtain the local maximum(36)

attained at(37)

where and we took into account that . The ratio of this peak to the main peak is(38)

Taking into account that (at V_A/V_b → 1), we see that the peak γ_m2 is always significantly smaller than the main peak γ_max. The maximum ratio γ_m2/ γ_max ≈ 0.45 is achieved at and V_A/V_b ≲ 1.

We note that the unstable fluctuations also have a small oscillatory part Re[ω] = 0.5 (n_b/n₀) ω_Bi. For most unstable wavenumber , the real frequency is smaller than the frequency of the normal Alfvén mode . Since , the instability is aperiodic.

5. Parallel Alfvén instability in particular compensated-current systems

Let us consider two feasible applications of CCPI. First we apply our results to the solar wind upstream of the quasi-parallel terrestrial shock where the plasma conditions are relatively well documented. Then we extend the analysis to the interstellar medium around supernova remnants, assuming the similar scalings of the beam parameters as in the terrestrial foreshock.

5.1. Quasi-parallel terrestrial foreshock

Hot ion beams with V_Tb > V_b > V_A are regularly observed in the solar wind upstream of the terrestrial bow shock where the shock normal is quasi-parallel to the interplanetary magnetic field B₀ (Paschmann et al. 1981; Tsurutani & Rodriguez 1981). This ordering of characteristic velocities suggests that CCPI driven by hot ion beams can develop in the quasi-parallel foreshocks.

More specifically, we use the following scalings for characteristic beam velocities, V_b ≲ V_shock and V_Tb ~ 3V_shock, where the shock velocity is equal to the solar wind speed, V_shock = V_SW. These scalings are compatible with observations reported by Paschmann et al. (1981) and Tsurutani & Rodriguez (1981). Yet another beam parameter, number density n_b, does not vary much around n_b = 0.1 cm⁻³ (Paschmann et al. 1981). In terms of the background solar-wind density n₀ ~ 5–10 cm⁻³, this gives n_b/n₀ ~ 0.01–0.02. Taking the typical value of Alfvén velocity, V_A ≈ 0.1V_SW, we obtain the cumulative destabilizing parameter α_b ~ 2.5–5, which is slightly over-threshold depending on the particular value of V_b. Such proximity of the system to the CCPI threshold can be a signature of CCPI operating in the foreshock and relaxing the beam parameters towards the threshold.

On the other hand, as is seen from Fig. 6, even slight deviations of α_b from the threshold can make CCPI strong. So, for V_A/V_b ~ 0.1 and α_b = 6 the maximum growth rate is already high, γ_max ≈ 0.07ω_Bi, with the most unstable wavenumbers k_zmρ_Tb ≳ 2. Narita et al. (2006) and Hobara et al. (2007) analyzed properties of electromagnetic fluctuations observed around terrestrial bow shock. Most straightforwardly, our results can be compared with the wavenumber distribution of the fluctuations in the quasi-parallel foreshocks shown in Fig. 9 by Narita et al. (2006), where the measured wavenumbers are normalized by the ion gyroradius. In terms of the background ion gyroradius ρ_Ti, with the typical temperature of the diffuse ions T_b/T_i = 4 × 10², our most unstable wavenumbers k_zmρ_Ti ~ k_zmρ_Tb/20 ~ 0.1 map upon the major peak observed at k_zρ_Ti = 0.1 (see upper panel of Fig. 9 in Narita et al. 2006).

In the quasi-parallel foreshock region, Narita et al. observed another, subdominant peak at k_zρ_Ti = 0.6. To explain this peak by CCPI we need a significantly lower beam temperature, T_b/T_i ~ 10, which is more typical for quasi-perpendicular foreshocks. We can speculate that CCPI can also generate this second peak. First, the CCP instability develops in the quasiperpendicular foreshock region where the beams have required temperatures T_b/T_i ~ 10, which is supported by the observed enhancement at k_zρ_Ti ≈ 0.4. Then the unstable fluctuations are convected in the quasi-parallel foreshock region where their observed wavenumbers are k_zρ_Ti ≈ 0.6.

The above estimations suggest that CCPI can contribute to electromagnetic fluctuations observed in the quasi-parallel terrestrial foreshock and can impose limitations on the parameters of the beams formed by reflected ions. Further direct confrontations of observed values of α_b with the stability diagram in Fig. 1 are needed to clarify the role of CCPI in the regulation of ion-beam parameters in the foreshock.

5.2. Foreshock regions around supernova remnants

Supernova remnants expanding in the interstellar medium develop bow shocks at their boundaries. These shocks propagate with high velocities V_shock ~ 2 × 10⁹ cms⁻¹ providing a feasible source of energy for the cosmic ray acceleration, and also for the magnetic field amplification. By analogy with the terrestrial bow shock, we assume that the reflected ions also occur in the supernova foreshocks setting up a compensated-current system. CCPI can develop in supernova foreshocks if parameters of reflected ions (subscript b) satisfy , defined by (19).

For reasonable background density n₀ = 10⁻²− 1 cm⁻³ and magnetic field B₀ ~ 10⁻⁷ − 10⁻⁵ G (Zweibel & Everett 2010), the Alfvén velocity varies in the range V_A = 2 × 10⁴−2 × 10⁷cm/s. Then the resulting Alfvén Mach number in supernova remnants M_A = V_shock/V_A = 10²−10⁵ is much larger than in Earth’s bow shock. For the similar scalings as in the terrestrial foreshocks, n_b/n₀ ~ 0.01, V_b ~ 0.5V_shock, and V_Tb ~ 2V_shock, even with the most unfavorable V_shock/V_A = 10² the destabilizing parameter α_b ~ 10² is much larger than the threshold . In this far over-threshold range, the CCPI operates in the asymptotic regime (35) with very high growth rate γ_max/ω_Bi ~ 0.5. We note that this value is already at the edge of applicability of our low-frequency approximation. Such a high growth rate suggests that the instability modifies the beam parameters strongly, in particular reducing the beam velocity towards the local Alfvén velocity, V_b ≳ V_A.

Let us compare the instability driven by the reflected ions with the similar instability driven by cosmic rays around supernova remnants (Bell 2004; Zweibel & Everett 2010). Taking the background magnetic field B₀ ≳ 10⁻⁶ G and the cosmicrays flux n_CRV_b ~ 10⁴ cm⁻² s⁻¹ (Zweibel & Everett 2010), we estimate the normalized current and the corresponding growth rate around supernova remnants. With ω_Bi ≃ 0.03 s⁻¹, we get in absolute numbers.

The above estimations show that the CCPI instability driven by reflected ions is much stronger than the instability driven by cosmic rays. Therefore, the former instability can be a more efficient amplifier for magnetic fields around supernova remnants. On the other hand, a fraction of the beam ions can be scattered back to the shock by electromagnetic fluctuations generated by CCPI, thus providing a seed population for the further Fermi acceleration to high cosmic-ray energies.

6. Discussion

A number of competing electrostatic and electromagnetic instabilities may arise when different plasma species move with respect to each other (see Gary 2005; Bret 2009, and references therein). The hierarchical structure of these instabilities depends on many parameters and remains an open question (see further discussions in Bret et al. 2010; Brown et al. 2013; Marcowith et al. 2016).

In our setting with hot ion beams, the fast two-stream/ Buneman instabilities are quenched by the large thermal velocities, which are larger than the streaming velocities. Inspection of Fig. 3.20 by Gary (2005) shows that the thresholds of electrostatic ion-acoustic and ion-cyclotron instabilities are significantly higher than the Alfvénic threshold for V_Ti/V_A ~ T_e/T_i ~ 1 typical in the terrestrial foreshock. Among them, the electron/ion cyclotron instability has the lowest threshold velocity, which is still very high, for n_b < 0.1n_e. The ion/ion acoustic instability is suppressed further by large beam temperatures, as is seen from Fig. 3.15 by Gary (2005). Therefore, these high-frequency electrostatic instabilities cannot compete with CCPI in the wide range of beam velocities 1 < V_b/V_A < 10². At higher beam velocities, V_b/V_A > 10², the ion-acoustic and ion-cyclotron harmonic waves can be generated by the electron-ion relative motion. However, even in this velocity range CCPI can develop independently as long as the mean parameters reside in the unstable area (Fig. 1), whereas the kinetic instabilities are quickly saturated by the local quasi-linear plateaus.

Parallel-propagating left-hand and right-hand polarized instabilities have been studied by Gary et al. (1984) and Gary (1985). Using numerical solutions of the dispersion equation, it has been observed that the left-hand polarized Alfvénic instability becomes competitive or even dominant when the beam ions are sufficiently hot (see Fig. 8 by Gary et al. 1984). The condition |ξ_b,−1| < 1 was used by Gary et al. to categorize this instabilas ion-beam resonant, i.e., resulting from the direct resonant coupling of the unstable mode with the beam ions. However, kinetic and reactive effects have not been distinguished for this mode, which did not allow us to realize that above threshold (19) the instability transforms from kinetic resonant to reactive non-resonant (see Fig. 8 and related discussions below). In the reactive regime, the meaning of the condition |ξ_b,−1| ≈ |k_zρ_Tb|⁻¹ < 1 is reversed: here it indicates that the unstable perturbations are on small enough scales to decouple from the beam ions by the demagnetization effect. The resulting Alfvén instability is then driven not by the resonant interactions with the beam ions, but by the bulk return current of the magnetized electrons. The current nature of this instability is similar to the nature of related current instability (Malovichko & Iukhimuk 1992) that can develop in the absence of any beams.

Fig. 8. Contribution of the reactive CCPI growth rate (dashed curve) to the total growth rate (dotted curve) for VTi/VA = Te/Ti = 1, Vb/VA = 10, VTb/VA = 25, and nb/n0 = 0.02. It is seen that the reactive destabilizing effects dominate the instability growth rate for this set of parameters. The wave frequency is shown by the solid line.

Interplay of the reactive and resonant left-hand Alfvénic instabilities also needs further investigations. Our preliminary estimations indicate that the relative importance of the reactive versus kinetic destabilizing effects increases quickly once α_b rises above the threshold . In Fig. 8. we show the contribution of the reactive CCPI to the total growth rate for reference plasma parameters that may occur in foreshocks: V_Ti/V_A = T_e/T_i = 1, V_b/V_A = 10, V_Tb/V_A = 25, and n_b/n₀ = 0.02. The corresponding total growth rate in Fig. 8 is given by Eq. (14) with the imaginary part of J₊ taken into account. It therefore includes both the reactive effects due to the bulk currents and the resonant wave-particle interactions. It is seen that the destabilizing reactive response becomes stronger than the resonant wave response when α_b is still not far from the threshold ( in Fig. 8). The instability is thus driven mainly by the reactive effects and can be analyzed ignoring kinetic resonant effects, as we did in the present study. The same approach can also be applied in the immediate vicinity of the reactive threshold if the quasi-linear plateaus or other local deformations of the velocity distributions weaken destabilizing kinetic effects. Analytical treatment becomes more tangled when reactive and kinetic effects are of comparable efficiency and have to be accounted for simultaneously, in which case the evolution of the system becomes more complex (cf. Yoon & Sarfraz 2017).

There are also left- and right-hand polarized instabilities driven by cold ion beams in the ion-cyclotron frequency range (Mecheri & Marsch 2007). These instabilities are strong when the velocity spread of the beam ions is so small that all the beam ions (and hence the beam as a whole) are resonant. In our settings with hot ion beams these instabilities are quenched similarly to the two-stream/Buneman instabilities.

In the considered case of hot ion beams, V_b/V_Tb < 1, the analytical treatment of wavenumbers k_zρ_Tb < V_Tb/V_b is simplified by neglecting the term ~V_b/V_Tb in ξ_b,−1. As the most unstable wavenumber scales as k_zρ_Tb ≈ α_b/2 (32), this restriction is not stringent:(39)

This condition is the opposite of the firehose instability condition (see Eq. (14) by Malovichko et al. 2014), which means that the CCPI can operate in a wide range of parameters below the firehose threshold. For cooler beams, where the condition V_b/V_Tb < 1 is violated (e.g., in the quasi-perpendicular foreshock regions), the analysis should be extended by accounting for corresponding terms.

The compensated-current parallel instability can also affect other processes in space. For example, it can limit the field-aligned currents generated by Alfvén-wave fluxes in the inner magnetosphere and plasma sheet boundary layer (Artemyev et al. 2016). In the solar wind, CCPI can contribute to the regulation of relative motion of different plasma species. It was found that many states of beaming structures in the solar wind are close to the thresholds of magnetosonic and Alfvén instabilities (Marsch & Livi 1987; Gary et al. 2000) and the firehose instability (Chen et al. 2016). Since CCPI can operate close to these thresholds (and sometimes below them), a refined analysis is needed to decide its role in the solar wind, as compared to the magnetosonic and firehose instabilities. These are subjects for future studies.

7. Conclusions

We investigated reactive non-resonant compensated-current parallel instability (CCPI) of left-hand polarized Alfvén waves in compensated-current systems established by hot diluted ion beams. Ion-beam demagnetization due to finite k_zρ_Tb is crucial for CCPI (ρ_Tb is based on the parallel beam temperature, and thus does not represent the beam ion gyroradius). New analytical expressions for the instability growth rate (31) and threshold (19) are found and analyzed.

The most important new properties of CCPI can be summarized as follows:

1. Reactive non-resonant CCPI depends on all bulk parameters of the beam: beam density n_b, bulk velocity V_b, and thermal velocity V_Tb. These parameters all increase the instability growth rate and can be combined in the single destabilizing parameter α_b = (n_b/n₀) (V_b/V_A) (V_Tb/V_A). The instability develops at , where the instability threshold (19) varies from at V_A/V_b → 0 to at V_A/V_b → 1. The analytical threshold (19) can be directly compared with satellite data to analyze the stability of beam-plasma systems in space.

2. CCPI is strongly affected by the velocity spread of the beam ions V_Tb. It defines the range of unstable beam currents, , with the current threshold varying in the range .

3. The instability growth rate γ_max (33) increases sharply with V_Tb once the threshold (20) is overcome (Fig. 7). This fast increase is caused by the fast demagnetization of the beam ions, in which case they cannot compensate for the perturbed currents of fully magnetized electrons. In a more distant overthreshold range the temperature dependence weakens because of the nearly saturated demagnetization.

4. From the growth rate γ_max (31) it follows that the instability can be strong, γ_max ≳ 0.1ω_Bi, even for modest not far from the threshold. The most unstable wavenumber k_zρ_Tb ≳ 1.54 near the threshold , but increases with α_b quickly approaching the asymptotic scaling k_zρ_Tb ~ _b/2. In this asymptotic regime, our growth rate reduces to (35), the same as was obtained by Bell (2004).

5. Two particular applications to the terrestrial foreshocks and supernova remnants show that the reactive CCPI can operate there. An analysis of Sect. 5.2 suggests that the ions reflected from the shocks around supernova remnants can drive stronger instability than the cosmic rays. In the terrestrial foreshock, CCPI can regulate beam parameters generating electromagnetic fluctuations observed at k_zρ_Ti ≈ 0.1.

Our results complement and extend previous studies on electromagnetic instabilities and their role in space and astrophysical plasmas. CCPI can develop around supernova remnants expanding into interstellar medium, participating in the braking process, and heating and redistributing energy in the supernova shocks. The same concerns the solar-wind regions upstream of the terrestrial bow shock, and other heliospheric shocks, where CCPI can bound the beam parameters and contribute to the low-frequency electromagnetic turbulence. Similarly, CCPI can affect other space and astrophysical environments containing super-Alfvénic ion beams and return currents.

Acknowledgments

This research was supported by the Belgian Science Policy Office (through Prodex/Cluster PEA 90316 and IAP Programme project P7/08 CHARM).

References

All Figures

	Fig. 1. Instability threshold in the parameter space (αb, VA/Vb) (solid line); the CCPI develops at all . The dashed line shows the split threshold above which there are two separate ranges of unstable wavenumbers kz
In the text

	Fig. 2. Illustration of the condition (18) for VA/Vb = 0.9 and αb = 6. In this case and there is only one unstable wavenumber range (shaded area).
In the text

	Fig. 3. Illustration of the condition (18) for VA/Vb = 0:9 and αb = 10. In this case and there two unstable wavenumber ranges presented by two shaded areas.
In the text

	Fig. 4. Unstable wavenumber ranges in the (αb, kz) plane for VA/Vb = 0.9. The outer boundary is defined by the left-hand side and the inner boundary by the right-hand side of condition (18). It is seen that below there is no instability, at there is a single unstable range of kz, and above there are two unstable ranges.
In the text

	Fig. 5. Wavenumber dependence of the instability growth rate driven by super-Alfvénic ion beams with VA/Vb = 0.9 for three values of αb: αb = 6, 8, and 10. For larger αb, the unstable area and the maximum growth rate extend to larger kzρTb.
In the text

	Fig. 6. Normalized growth rate γmax/ωBi as a function of nb/n0 and Vb/VA for hot beam with VTb/VA = 102. γmax increases regularly with both nb/n0 and Vb/VA once the threshold is exceeded.
In the text

	Fig. 7. Normalized growth rate γmax/ωBi as a function of VTb/VA for nbVb/ (n0VA) = 0.05 (solid curve), 0.1 (dotted curve), and 0.15 (dashed curve). Starting from zero, the growth rate γmax increases rapidly with VTb, but this increase is quickly saturated. Larger currents nbVb/ (n0VA) result in larger γmax for all VTb.
In the text

	Fig. 8. Contribution of the reactive CCPI growth rate (dashed curve) to the total growth rate (dotted curve) for VTi/VA = Te/Ti = 1, Vb/VA = 10, VTb/VA = 25, and nb/n0 = 0.02. It is seen that the reactive destabilizing effects dominate the instability growth rate for this set of parameters. The wave frequency is shown by the solid line.
In the text