At Earth’s bow shock, like at other collisionless shocks, ions reflected back upstream may form a beam population in velocity space. The free energy of this backstreaming ion beam drives a range of kinetic instabilities in the foreshock. For relatively tenuous and fast (compared to the background Alfven speed) ion beams traveling parallel to the magnetic field, the fastest growing linear instability is the resonant right-hand instability (RHI) (Gary, 1991). For low beam densities, this mode is a low-frequency wave carried by the background and excited by a cyclotron resonance with the beam ions. As described in another article in this collection (Winske and Wilson, 2022), Peter Gary was a pioneer in the Vlasov theory of electromagnetic ion beam instabilities in space plasmas (Gary et al., 1984; Gary et al., 1985; Gary, 1991). Gary’s work is an important piece of a large body of research on the theory and observation of RHI waves in the foreshock. Here, we re-examine properties and dynamics of the RHI using a modern hybrid particle-in-cell (PIC) code and example data from NASA’s Magnetospheric Multiscale (MMS) mission. The hybrid PIC simulations show that a finite temperature of the background decreases the amplitude of non-linear structures that develop and that a relatively broadband spectrum of modes may be excited with a range of propagation angles oblique to the local magnetic field.

Many early studies of RHI were inspired by data from the International Sun-Earth Explorer (ISEE) spacecraft. The data showed abundant wave activity in the foreshock, and it was recognized that RHI and other low-freqency modes were associated with backstreaming ions (Hoppe et al., 1981). Some effects of these waves include modifying the transport and heating of ions in the upstream plasma (Lee, 1982). The RHI in particular may also drive ultra-low frequency (ULF) waves in the 30-s period range (Watanabe and Terasawa, 1984; Greenstadt et al., 1995), which couple to particles trapped in Earth’s radiation belts. The non-linear evolution of RHI was found to be a possible driver of large-amplitude magnetic pulsations observed in the forshock (Akimoto et al., 1993). More recent simulation studies of RHI have focused on the global context of (ULF) waves and their transport into the magnetosphere (Blanco-Cano et al., 2009; Palmroth et al., 2015; Kajdič et al., 2021; Turc et al., 2022). In addition, very high-resolution field and particle data are now available from NASA’s MMS mission (Burch et al., 2016). Meanwhile, laser-driven laboratory experiments have offered a means of taking detailed measurements of ion beam instabilities in a reproducible environment (Heuer et al., 2020).

Here, we revisit RHI with a hybrid PIC code to further explore the waves and non-linear structures associated with the instability. Because earlier work focused on moderately cool beams with v _thb ≲ v _A (Hoshino and Terasawa, 1985; Winske and Gary, 1986; Akimoto et al., 1993) we also consider the effects of warmer beam and background ion populations. The finite temperatures moderately reduce the RHI growth rate, and they tend to reduce the amplitude of non-linear structures. We also explore the spectrum of oblique waves excited by a parallel ion beam. A statistical analysis of Cluster spacecraft data uptream of the quasi-parallel bow shock showed that the beam-driven waves have a power spectrum peaked at oblique propagation angles (Eastwood et al., 2005). This is at first glance at odds with the fact that the fastest growing mode for the beam and plasma conditions is the purely parallel propagating RHI. We find a range of oblique modes are excited in 2D and 3D simulations, consistent with solutions of the hot plasma dispersion relation. The simulations predict a typical perpendicular scale length of the non-linear RHI structures to be ∼ 100 to 200 d _i , indicating that the structures are planet-sized and may significantly impact the magnetosphere.

In this section, we include a brief review of the resonant right-hand instability (RHI) and define the conventions we use in our analysis. The RHI is a solution of the dispersion relation for magnetized plasmas with a beam ion population traveling along the magnetic field. For the cases of interest here, the fastest growing mode is purely parallel propagating. For purely parallel modes, the relevant dispersion relation for a plasma where each species j has an isotropic Maxwellian velocity distribution is (Gary et al., 1984): 1 − k 2 c 2 ω 2 + ∑ j w p j 2 ω 2 ζ j Z ζ j + = 0 (1) ζ j = ω − k u j / k v thj (2) ζ j + = ω − k u j + ω c j / k v thj (3) where k is the parallel wavenumber, ω p j 2 = n j Z j 2 e 2 / ϵ 0 m j , v thj = 2 T j / m j , u _j is the bulk fluid drift velocity, and Z(x) is the plasma dispersion function. In the quasi-neutral cold plasma limit, this results in a polynomial relation D(ω, k) = 0 that yields four distinct unstable modes (Weidl et al., 2019) including the RHI.

To study the linear growth rates of the RHI for oblique propagation including finite beam and background ion temperatures, we consider numerical solutions (Montgomery et al., 1975) of the full hot plasma dispersion relation assuming drifting Maxwellian ion distributions (see, for example the appendix of Gary, 1991). We use the open-source New Hampshire Dispersion Solver (NHDS) (Verscharen and Chandran, 2018). We work in the rest frame of the background plasma, where the background ions have no net drift. The background plasma is taken to be uniform and consisting of ions of mass m _i and unit charge. We denote the background density n ₀ and the temperature T ₀. The background magnetic field is also uniform, in the positive x direction, and of strength B ₀. A drifting Maxwellian ion beam population is included with a density n _b , temperature T _b , and a drift speed u _b in the positive x direction. The Alfven Mach number of the beam is M _A = u _b /v _A , were we normalize the drift speed to the background Alfven speed v A = B 0 / μ 0 n 0 m i . For simplicity, we take beam ions to be the same species as the background ions. The electron population is charge- and current-neutralizing.

At Earth’s bow shock, the reflected ion population is characteristically low-density (with relative beam fractions less than a few percent) and fast (with Alfven Mach numbers M _A > 2). Under these conditions, the low-frequency wave spectrum is dominated by the RHI. Figure 3 shows the maximum growth rate of the RHI over a range of propagation angles for a beam of density n _b /n ₀ = .015 and an Alfven Mach number of M _A = 10, similar to typical parameters at Earth’s foreshock. The three curves show the growth rates for three different beam temperatures. The uppermost curve is a relatively cold beam, and the peak growth rate for parallel propagation agrees with the large M _A approximation for cold plasmas (Gary, 1978; Weidl et al., 2019), γ / Ω c i ∼ ( n b / 2 n 0 ) 1 / 3 . When finite temperature effects are included, the growth rate decreases as the beam temperature increases. This is because a smaller fraction of the beam is gyroresonant with the mode. The background temperature has a very minor effect on the growth rates as long as the background thermal speed is relatively low (v _th0 ≪ u _b ) as is typically observed at Earth’s foreshock. We note that while additional beam-driven modes, including oblique Alfven modes Daughton and Gary. (1998) and the non-resonant mode (e.g., Gary. (1991), Chen et al. (2022), exist over this range, the RHI is dominant for the relatively fast (M _A > 2) and tenuous beams that we consider here. In Section 3, we explore the non-linear development of the purely parallel-propagating RHI using a hybrid PIC code.

We begin this section with an example of foreshock waves driven by reflected ions observed by NASA’s MMS mission (Burch et al., 2016). The MMS data are plotted in Figure 1, showing a fairly typical foreshock crossing in an interval upstream of the quasi-parallel bow shock on 30 January 2019. The magnetic field components in Figure 1A show a background field dominated by the radial B _x component along with large-amplitude wave fields (mainly B _y and B _z ). The wave power is peaked near the local ion cyclotron frequency ( ∼ 0.01 H z ) . In the spacecraft frame, the bulk solar wind velocity is u _x ∼ −300 km/s ∼ −10 v _A , and the reflected ions appear as a more diffuse population centered near v _x ∼ 0 (see Figure 1C). The approximate relative density of the reflected population is n _b /n ₀ ≲ .05, and the relative parallel drift of the beam gives an Alven Mach number of M _A ∼ 10. These rough parameters, which are fairly typical of Earth’s foreshock, serve as the basis for the numerical simulations below.

Because the electrons are far from resonance with the RHI and the mode frequencies are well below the electron plasma and cyclotron frequencies, hybrid numerical codes (Lipatov, 2002; Winske et al., 2003) that treat the ions as a kinetic species and the electrons as a massless neutralizing background are suitable for studying the linear and non-linear evolution of the RHI. Here, we use a hybrid version of the particle-in-cell (PIC) code VPIC (Bowers et al., 2008; Le et al., 2021; Keenan et al., 2022) to model the RHI. Numerous earlier studies of ion streaming instabilities have used similar electromagnetic hybrid PIC codes (e.g.,Winske and Quest., 1986; Hada et al., 1987; Gary and Winske, 1990; Winske and Omidi, 1992; Akimoto et al., 1993; Dubouloz and Scholer, 1995; Hellinger and Mangeney, 1999; Wang and Lin, 2003; Heuer et al., 2018; Weidl et al., 2019; Holcomb and Spitkovsky, 2019). Because the simulations in this section are 1D in the x direction, they only allow the growth of the purely parallel propagating modes, which are the fastest growing modes for the parameters we use. Note that for 2D or 3D systems with oblique modes, damping on the electrons can become more important. While this effect is included in the linear dispersion solver, it is not captured by the fluid model of the hybrid code.

As in other hybrid PIC codes, the electron model in Hybrid-VPIC takes the form of an Ohm’s law for the electric field: E = − u i × B − 1 n e ∇ p e + 1 n e J × B + η J − η H ∇ 2 J (4) where quasi-neutrality imposes n = n _e = ∑ _s Z _s n _s (including a sum over species s of ions), the velocity u _i is the charge-weighted ion flow u_i = ∑ _s Z _s n _s u _s /n _e , and the current density is taken in the low-frequency approximation as μ ₀ J = ∇ × B. We use a system of units based on the background magnetic field B ₀ and ion density n ₀, with times normalized by the cyclotron frequency ω _ci = eB ₀/m _i and lengths given in terms of the ion inertial length d i = ( ϵ 0 m i c 2 / e 2 n 0 ) 1 / 2 . Because we use particle shapes that are sums of quadratics in each direction, we choose grid resolutions with Δx ranging from .25 to 1 d _i to avoid an unphysical numerical dispersion that occurs when low spatial resolution and high-order particle shapes are used in hybrid PIC codes (Stanier et al., 2020). For these 1D simulations, we include 2000 particles per cell for each ion population. The normalized resistivity η/(B ₀/n ₀) and hyper-resistivity η H / ( B 0 / n 0 e d i 2 ) are set to small values in the range of 1 × 10^–4 to 5 × 10^–3. For the simulations here, the electron pressure follows a simple isothermal closure, such that the electron pressure is given by p _e = n _e T _e , with T _e a constant (Le et al., 2016). Test simulations with an adiabatic electron closure showed no discernible differences.

As described in Section 2, our 1D simulations contain a uniform background plasma of density n ₀ and temperature T ₀ (we set T _i = T _e = T ₀) and magnetic field B ₀ in the x direction. To this is added a streaming population of beam ions of density n _b (recall an equal number density of electrons is implicit in the quasi-neutral assumption of the hybrid code), drift velocity in the positive x direction u _b , and temperature T _b . Figure 2 shows typical results from an RHI simulation with relatively cold background and beam populations. The three sets of panels in Figure 2 are at three different times over the course of simulation.

The leftmost panels are at time t*ω _ci = 30 when the RHI is nearing the end of a phase growth consistent with the linear instability. The background and beam density profiles as well as the total magnetic field magnitude in Figure 2A are relatively unperturbed. The ion phase space distribution in x−v _x space in Figure 2 shows a weak modulation of the beam ions (the population centered at v _x ∼ 10v _A ). Nevertheless, the RHI is here already strong enough to modulate the beam ions in perpendicular velocity space. The three panels of Figures 2J–L show the perpendicular velocity distribution of the beam ions in v _y −v _z space at the three locations marked by vertical dashed lines in Figure 2A. These distributions exhibit gyrophase bunching, with the beam ions undergoing motion in the plane perpendicular to the magnetic field. As in previous simulations and observations (Hoshino and Terasawa, 1985; Thomsen et al., 1985; Fuselier et al., 1986; Gary et al., 1986), the gyrophase bunched ions are out of phase with the wave magnetic field. In Figures 2J–L, the white line gives the direction of the wave magnetic field (B _y and B _z ), while the black line shows the direction of the local bulk beam velocity in the y−z plane. The field and beam velocity are roughly 90° out of phase.

At later time, the RHI waves steepen into non-linear features. Non-linear structures have been observed with beam populations upstream of the bow shock, and they have been identified as shocklets (Hoppe et al., 1981; Hada et al., 1987) or magnetic pulsations (Akimoto et al., 1993). As in the earlier simulations of Akimoto et al., 1993, the non-linear pulsations driven by the RHI are characterized by correlated magnetic field strength |B| perturbations (see the blue curve in Figure 2B) that are correlated with the density n perturbations (red curve).

Thermal velocity spreads of the background and beam ions that are not large compared to the relative drift speed do not qualitatively affect the linear properties of the RHI, though the finite temperatures moderately reduce the growth rates. The non-linear features that develop, however, are weaker in amplitude in our simulations with higher beam and background temperatures. We show example data from a simulation with a beam thermal spread v _thb = 3v _A and background thermal speed of v _th0 = v _A in Figure 3. For cold beams ions as in Figure 2, practically all of the beam ions can become bunched where the RHI waves steepen. In addition to bunching in gyrophase angle, these resonant beam ions can be dramatically slowed down in the parallel direction, even locally coming to a stop in the background frame (see Figure 2H). For hot beams that are more diffuse in velocity space, on the other hand, a relatively smaller fraction of the beams ions are near exact resonance with the RHI mode. As a result, a smaller fraction of the beam ions in Figure 3H are slowed by the wave fields. This results in a much less spiky beam ion density profile in Figure 3B than for the cold ions case. Furthermore, because the RHI couples to compressional modes, the higher background pressure weakens the amplitude of the non-linear features. For even higher background temperatures with v _th0 = 4v _A (not plotted), there are no discernible pulsations or spikes in the density or magnetic field profiles.

The beam ions do display strong gyrophase bunching in perpendicular velocity space, although the non-linear structures that develop are relatively weak for the warmer beam and background temperatures (see Figure 3). We display side-by-side in Figure 4 ion v _y velocity distributions from the MMS event and from the Hybrid-VPIC simulation of Figure 3. The MMS data in Figure 4A show the v _y distribution of ions over time, which may be taken as a proxy for distance x because the waves rapidly cross the spacecraft. A distribution in x−v _y phase space from the hybrid PIC model is displayed in Figure 3B covering a range with a similar wave phase and amplitude as the MMS data. Note that ∼ 90 ° phase shift characteristic of the RHI between the wave fluctuations (carried by the background ions) and the gyro-bunched beam ions is visible in both the MMS data and the hybrid PIC simulation data.

In this section, we consider 2D and 3D hybrid simulations to study the evolution of the RHI in multiple spatial dimensions. The 2D or 3D geometry allows the development of a spectrum of modes with k vectors oblique to the magnetic field. As visible in Figure 5, the RHI growth rate is relatively insensitive to the propagation angle out to ∼ 35 ° . Interestingly, a statistical survey of 30-s waves observed in the Earth’s foreshock found that the distribution of wave propagation angles was typically peaked at an oblique angle of ∼ 20 ° to the magnetic field (Eastwood et al., 2005). Previous hybrid simulations suggested that refraction of steepening waves driven by the fastest growing parallel propagating modes could explain the presence of oblique modes in observations (Dubouloz and Scholer, 1995). On the other hand, (Strumik et al. 2015) studied the development of ULF waves in the foreshock of a 2D hybrid global magnetosphere model with a quasi-radial IMF and quasi-parallel bow shock. They found that averaging over the spectrum of excited RHI modes at different propagation angles gave a spectrum similar to the observations. We consider this possibility in our simplified uniform beam simulations below.

To examine the spectrum of oblique modes, we consider a 2D simulation in the x−z plane of a uniform beam, building on early hybrid simulation work on 2D ion-ion beam instability growth (Winske and Quest, 1986). The simulation domain is of size L _x × L _z = 1024 d _i × 1024 d _i = 2048 × 2048 cells, the background (v _th0 = v _A ) and beam (n _b /n ₀ = .02, M _A = 10, v _thb = 4.3v _A ) ion populations are each sampled by 400 particles per cell, and the time step is δt = .01/Ω _ci . In Figure 6, we compare the growth rates predicted by numerical solution of the hot plasma dispersion relation in (a) to the growth rates extracted directly from the hybrid PIC simulation in (b). The growth rate in each case is plotted in terms of the parallel (k _‖ and perpendicular k _⊥ wave numbers. The peak growth rate occurs for parallel propagation and corresponds to the usual RHI. Consistent with Figure 5, there is a relatively broad range of oblique wave vectors with growth rates very close to the maximum. Because of this, a wide spectrum of waves with varying propagation angles is excited. The modes plotted are all almost exactly right circularly polarized. The real frequency increases by a factor of few beyond propagation angles of 40°, and those more perpendicular modes may connect to a different wave branch.

The magnetic field and density structures that develop after the RHI saturates are very similar in 2D and 3D. Figure 7 shows comparisons of (a, c) the magnetic field component B _z and the (b, d) the plasma density between a 2D and a corresponding 3D hybrid simulation. These simulations are similar to the above simulation, but with a domain of size of L = 512d _i = 512 cells in each spatial dimension. The wide spectrum of unstable oblique modes produces magnetic fluctuations with oblique phase fronts. To quantify this effect, we show in Figure 8 a Fourier power spectrum of the wave magnetic field component B _z plotted in Figure 7A. The symmetry of the spectrum for k going to −k is simply a result of the Fourier transform of the real function B _z . The peak in the power spectrum is at k _x d _i ∼ .1, corresponding to a characteristic wavelength of ∼ 60 d i . As noted by (Strumik et al., 2015), the fluctuation power density is spread over a range of different k vectors centered at k _z = 0. Averaging over this spectrum can explain the statistics of oblique ULF waves observed by (Eastwood et al., 2005) in Earth’s foreshock.

As another way of displaying the parallel and perpendicular structure of the saturated magnetic field fluctuations, we plot cuts of the magnetic field components in Figure 9 along the (a, c) parallel or x direction and the (b, d) perpendicular or z direction. Again, the characteristic parallel wavelengths here are ∼ 60 d i . While the parallel mode (k _z = 0) is fastest growing, the cuts in z show large variations in the perpendicular direction. The typical perpendicular length scales are ∼ 100 to 200 d _i and are associated with the excitation of a wide spectrum of oblique wave vectors in the original linear modes. Although not studied here, this spectrum of non-linear fluctuations in 3D can contribute to cross-field diffusion of ions (Kucharek et al., 2000).

Using a modern hybrid PIC code, we revisited the resonant right-hand instability (RHI), which is the dominant electromagnetic ion beam instability for parallel-propagating ion beams that are relatively fast and tenuous. The parameters for the simulations were inspired by MMS observations of typical quasi-parallel foreshock fluctuations. RHI is prevalent upstream of the quasi-parallel region of Earth’s bow shock and has been the subject of a large number of previous observational and theoretical studies, including important early works by Peter Gary (Gary et al., 1984; Gary et al., 1985; Gary et al., 1986; Gary and Winske, 1990; Gary, 1991). Here, we focused on properties of the RHI related to finite temperature effects and oblique propagation. The hybrid simulations show that warm (β > 1) background and beam ion populations produce less steep non-linear features than the cold populations assumed in many previous studies. Note the RHI itself scatters reflected beam ions effectively, and it is possible that the RHI scatters an initially cooler beam into the more diffuse beam with a larger velocity spread observed by MMS. In any case, a high-beta background plasma requires additional energy to be compressed, which explains the weaker non-linear compressional features observed in simulations with high plasma beta (particularly β ≫ 1, which is not typical of the solar wind at Earth’s foreshock).

Multi-dimensional (2D and 3D) hybrid simulations demonstrated that a wide spectrum of oblique modes is excited, in agreement with growth rates predicted by numerical solution of the hot plasma dispersion relation. The RHI instability growth rate is a fairly flat function of propagation angle out to ∼ 35 ° . The theoretical growth rates for the RHI agreed with the range of modes excited in the hybrid simulations, and the non-linear stage contained fluctuations with characteristic perpendicular length scales 2–3 times longer than the typical parallel wavelength. At Earth’s bow shock, typical parallel length scales would be 60d _i , which corresponds to .5 to 1.5 R _E (Earth radii) for typical solar wind densities. The corresponding perpendicular lengths scales are 2–3 R _E . These non-linear structures are therefore planet-sized, and they may impact the planetary magnetopshere.