Energetic non-thermal particles, or cosmic rays, are a major component of astrophysical plasmas next to magnetic fields, radiation, and thermal gas. Cosmic rays are usually sub-dominant in density but carry as much pressure as the thermal plasma background. In some cases, cosmic rays drift at faster speeds with respect to the normal modes' phase speeds of the background plasma. Because of this, cosmic rays are a strong source of free energy that causes new classes of kinetic or convective instabilities. Recent years have seen the development of intense analytical and numerical efforts to analyze the onset of an instability produced by the motion of these particles at fast bulk speeds: this is the streaming instability. The streaming instability has been applied to different space plasmas and astrophysical contexts like strong shocks, jets, or in interstellar and intergalactic medium studies. Streaming instabilities participate in the production of magnetic turbulence at scales corresponding to the gyroradius of the particles. By scattering off their self-generated waves, cosmic rays are coupled to the background thermal plasma. This mechanism is able to self-confine cosmic rays around sources and launch winds out of the disk of the galaxy, hence impacting galactic matter dynamics and ultimately the galactic star formation rate. We discuss a few science cases, which should be accessible in the near future for analytical calculations and numerical simulations.

## I. INTRODUCTION

Cosmic rays (CRs) are charged particles with suprathermal energies that pervade the interstellar and intergalactic space to finally reach the Earth. CRs have kinetic energies spanning from a few MeV (and perhaps less) to a fraction of Zeta ($1021)$ eV or ZeV. The origin of these particles, that is, the way they are accelerated, is still a mystery more than 100 years after their discovery by V. F. Hess, but it is generally assumed that it involves interactions between the charged particles and local magnetic fields. The streaming instability (SI hereafter) in space and astrophysical plasmas is intimately related to the physics of CR acceleration and propagation in the interstellar (or ISM) or intergalactic medium (or IGM) (Amato, 2011). In effect, CRs carry an important amount of free energy, and as relativistic particles, they tend to stream in the ISM at a speed close to the speed of light c. Yet, as charged particles they are forced to propagate along ambient magnetic field lines and therefore can adopt an anisotropic distribution (Bykov *et al.*, 2013). For all these reasons, CRs can rather easily trigger various plasma instabilities. Dipole terms in the anisotropy expansion produce a drift and lead to different branches of streaming instabilities discussed below, while quadrupole terms lead to pressure effects involved in the firehose/mirror instabilities. These are not discussed in this manuscript. Interested readers can report to Achterberg (2013) and Zweibel (2020) for specific discussions.

The streaming instability (SI) occurs because CRs stream at a mean speed larger than the speed of normal plasma modes in the ambient plasma. In astrophysics, owing to the natural scale fixed by the particles, that is, their Larmor radius, the modes in question can be well described by the magnetohydrodynamic (MHD) theory. Basically, the SI is triggered when the CR streaming speed exceeds the local Alfvén speed. $vA=B/4\pi \rho $ fixed by the ambient magnetic field strength B and the local mass density of the plasma *ρ*, and the background plasma reacts to reduce the streaming speed of CRs. (Notice that all quantities derived in this text are expressed in CGS Gaussian units.) As we will see below, the SI has two main regimes: a resonant regime and a non-resonant regime. Historically, the resonant regime was investigated first (Wentzel, 1969; Kulsrud and Pearce, 1969; Skilling, 1971; Lee, 1972; Skilling, 1975). More recently, the non-resonant streaming regime was developed in astrophysical contexts by Bell(2004) although it had been known in the space plasma community for a longer time (Gary, 1993).

Although this article focuses mainly on SI in *some selected* astrophysical-related problems (see questions below), we also shortly discuss some recent developments for this topic in the context of space plasma environments.

Despite established theoretical frameworks, several issues still remain to be clarified in connection with the streaming of CRs: *In shock physics studies*: What is the saturation magnetic field level expected from SI? (Sec. III C); What are the properties of the turbulence resulting from the triggering of SI? (Secs. III C and V A); How does the relative importance of resonant and non-resonant modes evolve with CR injection rate and shock Alfvénic Mach numbers? (Secs. III C and V A); and What is the related time-dependent maximum CR energy? (Sec. V A). *In jet studies*: How much does the SI contribute to the generation of turbulence necessary to accelerate non-thermal particles in order to compensate their strong radiative losses? (Sec. V B). *In interstellar medium studies*: What are the physical processes responsible for the saturation of the SI? (Sec. III B); How much do the resonant modes contribute to the driving of galactic winds? How much CRs can penetrate into molecular clouds if they can trigger resonant SI (RSI) and, hence, to what extent does this explain the evolution of the ionization rate with the density column? (Sec. V C). All these questions are addressed in this article.

The article has the following layout. Section II provides a general overview of the SI: its causes and the different procedures to find its linear growth rate. Section III discusses recent numerical techniques and setup used to investigate the SI. It also discusses their main findings. Section IV shortly reviews recent studies about the SI in the space plasma community. Section V proposes some possible developments in different astrophysical contexts where the SI is expected to play an important role.

## II. THE STREAMING INSTABILITY LINEAR GROWTH RATE IN THE RESONANT AND NON-RESONANT REGIMES

### A. General statements

As for all kinetic instabilities, the SI growth rate in collisionless plasmas can be derived in a standard way from the Vlasov equation. [One should also mention the derivation obtained in Melrose (1968) based on the resolution of Maxwell equations in the cold fluid theory. In this formalism, waves are treated in a semi-classical formulation as plasmons with momenta $\u210fk\u2192$. Calculations are done in the quasi-linear theory (QLT) limit where particle and wave interactions are treated as perturbations. The growth rate is derived using the quantum detailed balance theory.] We have

where $F(r\u2192,p\u2192,t)$ is the CR distribution function in the phase space ($r\u2192,p\u2192,t$), and $p\u2192$ and $v\u2192$ are the particle momentum and velocity, respectively. The electromagnetic Lorentz forces are included in the particle equation of motion

where *q* is the particle charge and $E\u2192$ and $B\u2192$ are the local electric and magnetic field vectors, respectively.

The growth rate then results from a linearization of the Vlasov equation (Lerche, 1967; Gary, 1993). It consists of splitting *F* as the sum of $F0+F1$ composed of the equilibrium *F*^{0} and the perturbed *F*^{1} distributions. The solutions are searched for proton $F1+$ and electron $F1\u2212$ in terms of the equilibrium distribution $F0\xb1$ (Gary, 1993) (the discussion refers here to an electron–proton plasma, but other particle species can either be added). The plasma is perturbed by circularly polarized perturbations, which must satisfy the Maxwell–Ampère equation. This requires us to insert the perturbed current proportional to $(F1+\u2212F1\u2212)$ into the Ampère equation and to derive the wave dispersion relation. The relation can be put into the form $k2=J(\omega )+iK(\omega )$, where k is the wave number. The functions *J* and *K* are reals, but the wave frequency $\omega =R(\omega )+i\Gamma $ is complex. Under the assumption that the growth (damping) rate $\Gamma \u226aR(\omega )$, it can be written as $\Gamma \u2243\u2212K(R(\omega ),0)/J\u2032(R(\omega ),0)$, where the prime denotes the derivative with respect to $R(\omega )$. It is a growth rate if $K(R(\omega ))<0$. The analytical derivation of the growth rate depends on the specific form of the unperturbed distribution function *F*^{0}. The latter leads to the different branches of the SI either resonant or non-resonant as discussed below. The amplitude of the dipole in the anisotropy expansion fixes the drift speed controlling the growth of the streaming instability. The resonant branch of the SI involves wavenumbers of the order of the inverse of the particle Larmor radius (see Sec. II B), while in the non-resonant branch the wavenumbers are larger than the inverse of the particle Larmor radius (see Sec. II C). We now discuss the resonant and non-resonant regimes of the instability more specifically.

### B. The resonant streaming instability (RSI)

A convenient way to express the growth rate in the resonant regime is to calculate it in the frame where the electric field of the waves in resonant interaction with CRs vanishes (Skilling, 1975). The approach is similar to the one detailed above. It starts from the Vlasov equation, but the electric component of the force is removed. The growth rate now reads

where $\omega c=qB/mc$ is the cyclotron frequency and *m* is the particle mass. We introduce $UB=B2/8\pi $ the magnetic energy density of the background magnetic field of strength B. Here, as a matter of simplification, we have derived the growth rate for forward linear Alfvén waves (so combining two oppositely circular polarizations) moving along the background magnetic field. We can verify that the condition for the instability to grow is indeed $\u2202\mu F<0$. The anisotropy is due to the perturbed distribution *F*^{1}, and it should be small in order to fulfill a necessary condition for the perturbation theory to apply. In Eq. (3), the cyclotron resonance condition encoded in the Dirac function in the integral selects a particular particle pitch-angle cosine $cos\u2009(v\u2192,B\u2192)=\mu $ associated with a given wave number *k*. This is why the instability is called resonant. The resonant pitch-angle cosine $\mu R=1/krg$, where $rg=pc/eB$ is the particle gyroradius. Scattering particles with a small pitch-angle cosine requires high wave numbers, which usually contain much less turbulent power. Most of the scattering effect occurs then at $\mu R\u223c1$, so at wave numbers $k\u223c1/rg$.

Turning back to the observer frame, the linear growth rate can be conveniently expressed in terms of the CR distribution gradient along the local background magnetic field (Wiener, Zweibel, and Oh, 2013). The derivative of *F*^{1} with respect to *μ* [$\u2202\mu F1$ in Eq. (3)] is proportional to the gradient $\u2202zF0$ (the background magnetic field is supposed to lie along z here); in the case, only forward perturbations are produced. Otherwise, a second term dependent on the momentum derivative of *F*^{0} has to be added. The analysis then involves the contribution of backward propagating waves, see Skilling (1975) for details. Expressed in terms of the CR density $nCR$, Eq. (3) leads to

$A(p)\u22641$ is a function dependent on the CR energy distribution, see Wiener, Zweibel, and Oh (2013). The parameter $\eta =\delta B2/B2$ is the turbulence level. The resonant wave number is $k=m\omega c/p$. The form given in Eq. (4) will be useful when discussing problems related to CR propagation in interstellar or intergalactic media where the non-resonant instability discussed below is less likely destabilized. This formalism is discussed in Secs. III and V C. Notice finally that right-handed resonant modes have been observed in a laser experiment by Heuer *et al.* (2018). This opens a new research field for the investigation of the SI in laboratory plasmas that can help in our understanding of space and astrophysical cases.

Astrophysical and space plasmas can span a wide range in temperature and can be partially, rather than fully, ionized. These two effects modify the above growth rate solutions. The growth rate is then reduced in the presence of neutrals in partially ionized plasmas (O'C Drury *et al.*, 1996). Neutrals can heavily damp resonant modes in the high-frequency (or wave number) regime because ion and neutral motions are decoupled. Thermal effects (Landau damping) damp oblique propagating MHD waves strongly (Foote and Kulsrud, 1979).

### C. The non-resonant streaming instability (NRSI)

An alternative approach particularly well adapted to shocks solves the instability growth rate in the frame where the CR distribution is isotropic, while the background thermal plasma now *drifts* in at a super-Alfvénic speed $ud$ with respect to it (Zweibel, 2003; Amato and Blasi, 2009; Zweibel and Everett, 2010). The growth rate is as above deduced from a weak damping analysis of the dispersion relation. This formalism naturally includes the non-resonant streaming branch. The modes in this regime grow faster than the resonant modes when the drift speed is high enough. The criterion to destabilize the NRSI is (Zweibel and Everett, 2010) as follows:

where $UCR$ is the CR energy density. Amato and Blasi (2009) give the linear growth rate of the streaming instability in its non-resonant, but also resonant regime in case CR distribution is drifting with respect to the background plasma. If the condition in Eq. (5) is verified, non-resonant modes grow with a maximum growth rate (Bell, 2004)

where the CR current is $JCR=qnCRud$. The maximum growth rate occurs at $kmax=\Gamma g,NR,max/ua\u226brg\u22121$. As the instability develops at small scales, largely outside of the cyclotron resonance, it is called non-resonant. The growth rate at $k\u2264kmax$ scales as $\Gamma g,NR=\Gamma g,NR,maxk/kmax$ (see Fig. 1). In this instability, the source of free energy involves the CR current [so including all the particle distribution and the non-resonant instability can also be investigated using an MHD approach as in Bell (2004)]. The instability is in fact due to the current $J\u2192b$ in the background plasma compensating the CR current. The Lorentz force $J\u2192b/c\u2227B\u2192$ contributes to inflate any seed magnetic loops and then leads to the instability (Milosavljević and Nakar, 2006).

The above estimation has been obtained for a CR current drifting parallel to the background magnetic field. More general dispersion relations can be found in the case of non-collinear current configurations (Bell, 2005). The non-resonant instability can even grow in the case $J\u2192CR\u22a5B\u2192$ unless the perturbed wave number $k\u2192$ is perpendicular to $B\u2192$. Non-aligned CR currents can occur because of CR diamagnetic drifts or in the case of inhomogeneous background magnetic field configurations (Riquelme and Spitkovsky, 2010; Nekrasov and Shadmehri, 2014).

The effects of neutrals have been treated in Bykov and Toptygin (2005); Reville *et al.* (2007) and plasma beta effects in Zweibel and Everett(2010). It is shown that neutrals cannot completely stabilize the NRSI. Reville, Giacinti, and Scott (2021) discuss the effect of the presence of CR streaming over the “evanescent band” of Alfvén wave propagation in partially ionized media. This non-propagating band appears if the ratio of the mass densities in ions to neutrals $\rho i/\rho n<1/8$. It is caused because a disturbance in the magnetic field decays due to ion–neutral friction before the magnetic tension is transferred to neutrals by ion–neutral coupling (Soler *et al.*, 2013). These authors show that if the current imposed by CRs is large enough, such an evanescent zone does not exists anymore. Thermal effects are also usually unable to stabilize the NRSI in most of astrophysical situations [see, however, recent hybrid simulation results by Marret *et al.* (2021)]. In hot media (with temperatures $\u223c107$ K that correspond to the IGM), the maximum growth rate is reduced and large k perturbations are destabilized with respect to lower temperatures cases.

## III. NUMERICAL STUDIES: STATUS AND PERSPECTIVES

Numerical studies of the SI are challenging because non-linear interactions between non-thermal particles—the source of the instability—and the background thermal plasma require a kinetic model for the non-thermal component to be properly described. This is true even in the case of the non-resonant instability because information on the non-thermal particle distribution is necessary to properly calculate their charge and current densities. In addition to this intrinsic difficulty, the numerical technique used to study the instability needs to be selected as a function of the physical problem to be investigated by fixing a proper setup. Then, care has to be taken while extrapolating numerical results to the astrophysical context by doing a proper rescaling of the physical parameters.

### A. Numerical techniques

The numerical techniques adapted to investigate non-thermal particle (or CR) back-reaction over the background plasma can be ordered by increasing physical scales.

#### 1. Thermal electron and proton scales: Particle-in-cell methods

Particle-in-cell (PIC) methods treat the plasma kinetically, describing both ions and electrons as a collection of particles and solve the equations of particle motion by iteration under the effect of the Lorentz force coupled to the Maxwell equations (Birdsall and Langdon, 1991; Pohl, Hoshino, and Niemiec, 2020).

#### 2. Thermal proton scales: Hybrid methods

In hybrid methods, electron dynamics are not treated kinetically but using a fluid model. Ions are treated using PIC techniques (Lipatov, 2002). Non-thermal particles can be included in the PIC solver (Gargaté *et al.*, 2007).

#### 3. Large scales including non-thermal component scales: Magnetohydrodynamics and kinetic

In these approaches, the thermal plasma is treated using MHD equations, whereas non-thermal particles are treated kinetically using PIC techniques (Bai *et al.*, 2015; van Marle *et al.*, 2018). Alternatively, non-thermal particles can be treated using a kinetic model like Vlasov (Palmroth *et al.*, 2018) or Vlasov–Fokker–Planck (Reville and Bell, 2013).

#### 4. Spatial coupling

Finally, it is possible to perform a simulation that combines both kinetic and fluid approaches, not by separating the plasma into two components, but by assigning different methods to different spatial locations. In this fashion, a simulation can use the PIC to simulate the formation and structure of a shock, whereas the large-scale plasma flows are treated through MHD for computational efficiency (Makwana, Keppens, and Lapenta, 2017, 2018b).

### B. The resonant regime

The resonant regime has been investigated using the PIC (Holcomb and Spitkovsky, 2019; Weidl, Winske, and Niemann, 2019b; Shalaby, Thomas, and Pfrommer, 2021), hybrid (Weidl, Winske, and Niemann, 2019a; Haggerty, Caprioli, and Zweibel, 2019), and PIC-MHD (Bai *et al.*, 2019; Lebiga, Santos-Lima, and Yan, 2018; Plotnikov, Ostriker, and Bai, 2021; Bambic, Bai, and Ostriker, 2021) approaches. The final results (linear growth rate, saturation level of magnetic fluctuations) depend on the simulation setup. The most recent advances have benefited from the above simulation techniques coming to a mature stage. They are discussed below. The interested readers should report to extensive reviews (see, e.g., Marcowith *et al.*, 2020) to have a more complete view of this rapidly growing field. Below, we subjectively selected two articles among the first to address linear and non-linear phases of the RSI. In particular, these works investigate the so-called 90° pitch-angle scattering problem faced by the quasi-linear theory (QLT) of particle transport in turbulent magnetic fields (Skilling, 1975; Shalchi, 2009). In effect, the QLT based on resonant wave–particle interaction predicts diverging parallel mean free paths as the particle pitch-angle cosine goes toward zero. This anomaly usually requires the introduction of non-linear corrections to the QLT by adding some perturbation of the particle trajectory through for instance the broadening of the cyclotron resonance in Eq. (3). The two different setups based on two simulation techniques are reviewed in order to cross-check their results.

Holcomb and Spitkovsky consider the case of the triggering of streaming modes by an anisotropic particle distribution composed of protons. For the purpose of their study, they consider two types of distribution: (1) a monoenergetic gyrotropic ring distribution, $fCR\u221d\delta (p\u2212p0)\delta (\mu \u2212\mu 0)$,has been widely adopted in the context of solar wind studies and (2) a power-law distribution, isotropic in a frame moving with a (non-relativistic) drift speed $v\u2192d$ with respect to the observer frame, is more relevant to astrophysical contexts, like supernova remnant shocks or in ISM studies. An advantage of the ring distribution is to isolate the wave–particle resonance in the momentum phase space as both momenta and particle pitch-angle are initially set.

Using 1D3V PIC simulations, the authors follow the growth of resonant modes for both types of particle distributions. While linear growth rates in case 1 are well reproduced, the authors notice that the effect of a drift speed $ud\u226buA$ breaks the degeneracy among right-handed and left-handed modes in the case 2, the former being generated in excess by cosmic rays (positive charges). The fastest right-handed (left-handed) growing modes are shifted toward smaller (larger) wave numbers with respect to linearly polarized case. In general, simulations in the case 2 produce smaller growth rates with respect to analytical estimates. The origin of this discrepancy is possibly numerical or due to the expression of the scattering frequency strictly valid in the QLT limit.

Then, two important aspects of the SI are discussed: the $90\xb0$ problem in the QLT and the magnetic field saturation level. In fact, both issues are related. The ability to cross the pitch-angle cosine *μ* = 0 strongly depends on the density of the injected CR with respect to the background gas, for both distribution types. The initial anisotropy marked by the amplitude of the drift speed is the second main controlling factor related to the *μ* = 0 crossing. For case 2, if the density is too low, the amplitude of self-generated perturbations is low and CRs are unable to cross the *μ* = 0 boundary because of lack of left-handed modes. Then, the isotropization process cannot proceed and the level of the magnetic field stalls. An initial low level of anisotropy favors the generation of left-handed modes and particle isotropization.

The authors propose different ways to evaluate the level of saturation magnetic field, a level that depends on the particle distribution setup. In the case 1, so for a ring distribution, particle trapping is the main saturation mechanism and the simulations agree with Shevchenko, Galinsky, and Ride (2002), as well as the ratio of the saturation magnetic field to the background one scales as $(nCR/nback)2/3$ down to $nCR/nback=10\u22124$. If this ratio is $\u226a1$, the saturation is not associated with an isotropization. In case 2, so for a power-law distribution, the level of fluctuations at a given scale is fixed by balancing the pitch-angle scattering frequency with the linear growth rate. The saturation magnetic field ratio scales as $(nCR/nback)1/2$. However, because isotropization is not effective at low $nCR/nback$ ratios, the drift speed remains always larger than $ua$. In summary, the linear physics is well reproduced. The crossing of the $90\xb0$ barrier and the anisotropy reduction are intimately related. This requires the triggering of appropriate mode polarization with a sufficiently high amplitude.

This work provides important insights into the non-linear physics of the RSI. It would then be worthwhile discussing how it could be transposed to more realistic astrophysical situations where the ratio $nCR/nback$ does not exceed $\u223c10\u22129$ in the ISM or $\u223c10\u22127$ inside or close to CR sources. In particular, in the former case, the initial anisotropy is low with $ud$ close to $uA$, a regime yet unexplored. One may wonder if the isotropization process could at all occur in these conditions. Conversely, it could also be interesting (and technically feasible) to explore the parameter space with high CR density and high drift speeds to test the relation between the saturation magnetic field and $nCR/nback$ for case 2. Longer term simulations (if technically possible) would probably be appropriate to handle the non-linear phase of the RSI.

The PIC simulations have some intrinsic limitations as they need to resolve both thermal electron skin depth scale $c/\omega p$ and CR scale $Rg$. In ISM conditions, these scales are separated by a factor $\u223cc/ua\u226b1$ (typically $\u223c106$). Hence, the PIC techniques can only investigate plasmas with high Alfvén speed. One way to alleviate this issue is to use the PIC-MHD method, as is the case in (Bai *et al.*, 2019; Lebiga, Santos-Lima, and Yan, 2018). PIC-MHD simulations can also follow the evolution of the system over longer timescales.

Using this approach, Bai *et al.* (2019) conducted 1D3V simulations where the mean CR distribution *f*_{0} follows a kappa distribution, which is continuous at all p (it is constant at small momenta and follows a power law above a threshold momentum). They use an offset distribution *δf*, which encodes the CR anisotropy and which is the driving source of the streaming modes. The calculation of this offset permits to reduce the numerical noise in the PIC module inserted in the MHD solver. We notice that using a similar technique the pioneering work of Lucek and Bell (2000) considers the case of a monoenergetic particle distribution drifting at a speed $u\u2192d$ with respect to a background plasma.

In Bai *et al.* (2019), the CR distribution has a small drift with $ud/ua\u22432$. Yet, simulations with a ratio up to eight and even ten as in Lucek and Bell (2000) have also been performed. As discussed above for the PIC simulations, low drift speeds produce equal growth of both right- and left-handed waves. PIC-MHD simulations confirm this result, and the theoretical linear growth rates are well reproduced. The particle isotropization depends on the density of CRs in each momentum bins and hence on the amplitude of the destabilized resonant modes. Low-energy CRs are isotropized more rapidly as for a power-law distribution with a negative index that they are able to trigger waves of both polarizations to higher amplitudes. In higher drift speed cases, the momentum limit where the drift speed is reduced to $ua$ at the end of the simulation increases with $ud$. The authors find that the mean drift speed (over the whole distribution) drops to $ua$ at the end of the simulation, a result in tension with the PIC simulations. It is instructive at this state to compare run Hi1 in Holcomb and Spitkovsky (2019) (initial drift speed $vd=7.9va$, CR density to background gas density $nCR/ni=2\xd710\u22124$, number of particle per cell: 250) with run vD8 in Bai *et al.* (2019) (initial drift speed $vd=8va$, CR to background gas density $nCR/ni=10\u22124$, number of particle per cell: 256). We notice the main difference between the two techniques in the ratio, $va/c$, and the simulation time. In PIC-MHD simulations, the magnetic saturation level taking into account numerical energy dissipation approaches the expected value obtained by balancing the momentum transfer from CR to waves.

Another interesting aspect addressed by the authors is the momentum feedback induced by CRs through wave generation over the background magnetized gas. The simulations clearly show the transfer of momentum from the CR to the gas with a good total energy conservation. This effect is likely contributing to the launch of winds from the galactic disk. Another aspect, which is still beyond the scope of current simulations, is the correct evaluation of the background gas heating through the triggering of the RSI. As stated by Holcomb and Spitkovsky (2019) and Bai *et al.* (2019), this requires to properly balance the wave generation by an appropriate damping: turbulent damping due to the interaction of self-generated waves with background turbulence, non-linear Landau damping in hot ionized plasmas, and ion–neutral damping in partially ionized plasmas. This important physical process will be discussed in Sec. III D.

The recent implementation of ion–neutral damping in PIC-MHD approach was presented by Plotnikov, Ostriker, and Bai (2021). In Fig. 2, we present the selection of results of that study, comparing the evolution of the instability without and with moderate ion–neutral damping rate (blue- vs red-colored lines in different panels). The ion–neutral damping reduces the linear growth rate of the instability (panel a) and the saturation level of waves $\delta B/B0$ (panel c) and, more importantly, removes low-*k* and high-*k* resonances where the growth rate is smaller than the damping rate (compare the blue and red curves in panel b). The first effect reduces the rate at which the most unstable modes reduce the drift velocity for particles resonating with most unstable modes (panel d). The latter effect removes entirely the isotropization of particles that are normally resonant with wavelengths, which are not able to grow with significant damping. The ion–neutral damping is important in some astrophysical environments where the accelerator (e.g., supernova remnant) is close to a partially ionized interstellar cloud (e.g., molecular clouds). The outcome of the instability defines how the CRs penetrate into the cloud (Morlino and Gabici, 2015; Ivlev *et al.*, 2018; Silsbee and Ivlev, 2019).

In summary, PIC and PIC-MHD appear to be appropriate tools to investigate the linear physics of the RSI, as well as non-linear wave particle interaction processes involved in its saturation. The two main parameters controlling the wave growth and the non-linear evolution are the CR to gas density ratio and the CR drift to ambient Alfvén speed ratio. Large speeds and density ratios break the degeneracy in the wave circular polarization production: right-handed waves peak at higher k with respect to the linearly polarized case. Both techniques reproduce linear growth rate well. Yet, the conclusions concerning the magnetic field saturation need still be fully explored and understood. The PIC simulations show incomplete isotropization at high-speed/low-density ratios, contrary to PIC-MHD simulations.

### C. The non-resonant regime

This regime has been investigated using PIC (Riquelme and Spitkovsky, 2009; Kobzar *et al.*, 2017), hybrid (Gargaté and Spitkovsky, 2012; Caprioli and Spitkovsky, 2014b; Haggerty, Caprioli, and Zweibel, 2019), and PIC-MHD (Bai *et al.*, 2015; van Marle *et al.*, 2018, 2019) methods in the framework of CR acceleration at shock waves and in Schroer *et al.* (2021) in the context of ISM studies. The NRSI instability is interesting as it allows to quickly generate magnetic field perturbations, useful to confine particles around the shock and hence to increase their maximum energy. Unfortunately, as stated in Sec. II C, these perturbations are triggered at wave numbers much in excess with respect to the inverse of particle gyroradius. An essential aspect is then to investigate the instability non-linear phase and the different ways that the turbulent energy can be redistributed to larger (resonant) scales. We have thus selected two works, which were among the first to investigate the parameter space under which the RSI and NRSI can be triggered at shock waves. These two works have then also investigated the non-linear phases. As for the RSI case above, we consider two different numerical techniques in order to cross-check their results.

Caprioli and Spitkovsky (2014b) evaluate the respective role of resonant and non-resonant streaming modes at fast strong shocks. They first discuss the parameter domain range where either the resonant or the non-resonant instability dominates, that is, where the linear growth rate of a given instability is the fastest. A criterion for the dominance of the non-resonant branch is (if the CR distribution is isotropic in the shock rest frame and for a CR distribution scaling as $p\u22124$) $Ma\xi CR(mpc)>1$. Here, the Alfvénic Mach number is $Ma=ush/va,\u2009\xi CR(p)=nCR(>p)/ng$, and $ush,\u2009nCR(>p)$ and $ng$ are the CR density above a momentum p, the shock speed, and the ambient gas density, respectively. Then, they evaluate the level of saturation of the non-resonant instability to be $\delta B/B0\u223cM0/2$, and the Mach number *M*_{0} is associated with the maximum CR energy $Emax$ and the energy at which CR is injected into the shock process $Einj$ through $M0=MaEmax/Einj$. The authors find a Mach-number-dependent turbulent spectrum. Shocks with Mach numbers below 30 have a turbulent spectrum scaling as $k\u22121$ consistent with the expectation from the resonant instability. Above this value, magnetic field amplification is enhanced and non-resonant modes are driven more rapidly. The power spectrum is also found to be dependent on the magnetic field obliquity, that is, when the upstream magnetic field is not parallel to the shock normal, so makes an angle $\theta B$ with respect to it. The authors evaluate particle injection and acceleration efficiencies at oblique shocks. Here, the shock drift acceleration (SDA) process due to the electric field carried by the plasma advected toward the shock front starts to play a significant role in the injection. SDA allows supra-thermal particles to cross the shock from downstream to upstream because the acceleration produced by the convective electric field. However, it becomes more and more difficult for the particles to overcome the electrostatic barrier at the shock front as the obliquity angle goes toward 90°. It results that the injection of supra-thermal particles and hence at least proton acceleration ceases above $\theta B\u223c55\xb0$. van Marle *et al.* (2018) propose a Fourier analysis of the non-resonant instability in high-Mach-number shocks and subsequent turbulence triggered in the upstream medium (see Sec. III C). The time-dependent development of the non-resonant instability is well captured by the Fourier analysis (see Fig. 3). It shows at early timescales that the development of the fastest wavenumber $kmax=JCR/2B0$, which is shifted by a factor 4 in the downstream medium due to shock compression. Later times show the development of a turbulent medium up- and downstream. van Marle *et al.* (2019) then investigate 3D parallel shock evolution and find similar results (even if the magnetic field growth is slower in 3D geometry with respect to 2D results and hence a smaller particle acceleration efficiency). The long-term evolution of the shock also shows a strong corrugation of the shock front and hence of the magnetic field obliquity.

Contrary to the above results obtained for an equipartition between upstream gas and magnetic field pressures, van Marle (2020) performed 2D PIC-MHD simulations of low sonic Mach shocks ($M=2\u22124$) (but high Alfvénic ones, $Ma=20\u221240$) in high-beta plasmas, a setup more adapted to galaxy cluster shock waves. The author finds a transition regime in Mach number for particle acceleration rather a critical value: above $Ma$ = 2.5 efficient acceleration occurs and a power-law distribution is built up with an index consistent with the test-particle solution, that is, scaling as $p\u22123r/(r\u22121)$, where r is the shock compression ratio. In parallel, strong instabilities are produced for M = 3–4 shocks. However, the nature of these instabilities has not been certified yet in this study.

PIC-MHD and PIC or hybrid simulations converge in their results in the case of high Alfvénic numbers (up to 30) parallel shocks. Some differences were, however, found in the particle acceleration process at high magnetic field obliquities: van Marle *et al.* (2018) still found some evidence of particle injection and acceleration at $\theta B\u223c70\xb0$ but not above, whereas Caprioli and Spitkovsky (2014b) evaluated this limit at $\theta B\u2009\u223c\u200955\xb0$. Haggerty and Caprioli (2019) resolved this discrepancy by using extensive hybrid simulations, which allow to include relativistic ion dynamics and accurate estimations of the injection rate (which PIC-MHD codes cannot). They confirmed the $55\xb0$ limit. In Sec. V, we discuss further investigations about this issue.

### D. Perspectives in numerical studies

Simulations have already started to explore the driving of SI in a large variety of plasma background conditions.

Although a purely kinetic approach (PIC) remains the most realistic because of the microphysics included in the model, it comes with severe computational penalties in terms of the number of processors and the amount of memory space required to run the simulations. This not only restricts the size of the spatial and temporal domain that can be modeled, but also restricts the physical conditions that can be simulated because issues like signal travel time and shock formation time can become prohibitive. It can also create problems for the accuracy of the simulation because effects with longer wavelengths, effects that are strictly multi-dimensional, or effects that only appear over long periods of time are difficult, if not impossible, to capture and may simply be lost in the model.

The alternatives, whether it be a multi-fluid approach, di-hybrid, or a coupling between PIC and MHD, all sacrifice accuracy in physics for greater computational efficiency. This allows them to explore those effects that can only be observed on scales that are out of reach for the PIC method, given current computational limitations. Whether the physics aspects, that is, lost matter, vary depending on the problem under investigation. For example, for the propagation of CRs through the ISM, the physics of the methods described above are sufficient. However, when looking at the production of non-thermal particles at a shock front, both the multi-fluid approach and coupled PIC-MHD lack the necessary physics and instead rely on an *ad hoc* estimate based on the PIC results, whereas the PIC results should be sensitive to the effect of large-scale fluctuations induced by high-energy particles.

Ideally, a combined approach could be used, where the PIC method is used to calibrate the input for larger scale models; for example., van Marle (2020) performed PIC-MHD simulations of high-*β* shocks using an injection rate for non-thermal particles that was derived directly from the PIC simulations for the same physical parameters done by Ha *et al.* (2018). The results of both sets of simulations appear to be generally in agreement but show quantitative differences brought on by large-scale, long-term side effects that the PIC simulations could not resolve.

Finally, the properties of background plasma, the degree of ionization and temperature, need to be integrated into the parametric survey of the long-term evolution of the instability, so beyond the evaluation of the linear growth rate (see recent efforts in that direction by Marret *et al.* (2021); Reville, Giacinti, and Scott(2021)). In order to include long-term dynamics and large-scale evolution, MHD codes should include neutral species models combined with fluid or kinetic treatment of CRs.

## IV. STREAMING INSTABILITY STUDIES ASSOCIATED WITH SOLAR/STELLAR ACTIVITY

### A. Coronal mass ejections

The resonant SI is of particular relevance in the problem of acceleration and transport of solar energetic particles (SEPs). We notice that solar wind speeds do not likely allow conditions in Eq. (5) to be fulfilled; hence, the SI regime is restricted to the resonant case. SEPs are high-energy particles coming from the Sun. They are composed of protons, electrons, and ions with energies ranging from a few tens of keV to a few GeV. SEP acceleration is associated with two types of eruptive processes: solar flares and coronal mass ejections (CMEs) (Verkhoglyadova, Zank, and Li, 2015). CMEs are believed to produce the gradual events and the highest SEP energies (Kouloumvakos *et al.*, 2019; Reames, 2021). One of the favorite mechanisms to accelerate these particles is the diffusive shock acceleration (DSA) process (Blandford and Ostriker, 1978), also known as Fermi 1 acceleration, where energetic particles gain energy in successive shock crossings induced by scattering off magnetic perturbations in up- and downstream media. The upstream turbulence necessary to scatter particles back to the shock front can be well produced by the particle themselves if they carry enough free energy.

One difficulty is how to explain the delay and the flux of these energetic particles at 1 AU (astronomical unit) after the onset of a solar eruption and prior to the shock arrival. For this to apply, SEPs have to escape upstream the shock, which requires a sufficiently large mean free path. Observations also indicate that the particle mean free path is energy dependent: low energy particles show a reduced, energy-independent mean free path, while high-energy particles show larger, energy-dependent mean free paths (Palmer, 1982). This may suggest that accelerated particles need to have a minimum amount of anisotropy to trigger a SI. In order to investigate this complex problem, several models have been proposed, which combine particle acceleration and transport around the shock front and the generation and transport of magnetic perturbations, namely, Alfvén waves self-triggered by the particles (Lee, 1983; Gordon *et al.*, 1999; Vainio, 2003). The basic assumption of the models is comprehensively detailed in Lee (1983): the shock front is treated as planar, the upstream medium (the solar wind) has a bulk flow *u* with respect to the shock front, energetic particles are drifting parallel to the background magnetic field at speeds $v\u226bu$, and they trigger left- and right-handed circularly polarized waves along the background magnetic field. Then, a time-independent system of two coupled equations is solved for the particle distribution *F*(*p*, *z*) as a function to the distance z to the shock front upstream and for the wave amplitude *I*(*k*, *z*) defined by $\u27e8\delta B\u27e92=\u222bdkI(k)$. The quasi-linear theory is further assumed in order for Eq. (4) to be applied. The system reads as

where $\theta B$ is again the angle between the shock normal and the background magnetic field, $uA$ is again the local Alfvén speed, and *κ* is the parallel (spatial) diffusion coefficient of SEP to the background magnetic field.

However, time-dependent models are necessary to understand time profiles of the particle fluxes for different species associated with CMEs. Time-dependent effects can, in general, only be investigated using numerical simulations (Li, Zank, and Rice, 2003) or assuming specific simplified analytical models. Lee and Ryan (1986) generalized the calculations by Lee (1983) including time effects but with momentum-independent diffusion coefficients in spherical geometry. Lee (2005) further generalized the theory to oblique shock geometries. Another approach has been proposed in the framework of the onion-shell model (Bogdan and Volk, 1983; Zank, Rice, and Wu, 2000). In this model, analytical spherical solutions are used to derive upstream particle and self-generated wave distributions (which in the QLT scale as 1/x, where x is the distance to the shock front). In the test-particle limit, the particle distribution at the shock is a function of the shock compression ratio, *r*(*t*), at a given time *t* as $F(p)\u221dp\u22123r(t)/(r(t)\u22121)$. This solution is used to derive the downstream solution, accounting for convection and adiabatic expansion losses of the plasma. The particle maximum energy is obtained by competing dynamical and shock acceleration timescales. The onion-shell model serves as a basis to numerical model of SEP acceleration and transport (e.g., the PATH code Rice *et al.*, 2003) where particles escaping from the shocked medium are tracked using a Monte Carlo approach. The model has provided satisfactory fits of large SEP events with the simulated time intensity signals and event spectra (Verkhoglyadova *et al.*, 2010). The code has been since extended to account for a 2D geometry (Hu *et al.*, 2017). However, Alfvén wave-generated dynamics are not treated in the above models. This aspect is important to account for the escape process of the highest SEP energies (Zank, Rice, and Wu, 2000).

Vainio (2003) proposes an explicit treatment of the time evolution of the self-generated turbulence spectrum. Vainio and Laitinen (2007) present a Monte Carlo simulation software (CSA code), which does include wave generation through particle streaming based on the previous analytical treatment. The code solves a system of coupled 1D kinetic equations equivalent to Eqs. (7) and (8) but time-dependent and using the stochastic differential equation formalism (hence, the Monte Carlo character of the code) (Kruells and Achterberg, 1994; Marcowith and Casse, 2010; Strauss and Effenberger, 2017) in order to reconstruct the evolution of the particle distribution function. These studies find that the amplitude of self-generated waves is the greatest close to the Sun, whereas, as the CME moves outward, the wave intensity decreases and the particle mean free path increases.

We notice finally the work done by Lange *et al.* (2013), which uses an hybrid technique (the GISMO code) combining incompressible MHD and a particle-in-cell module tracking the energetic particles (see Sec. III for a discussion about this technique). The turbulence is driven anisotropically mimicking the effect of particle streaming. In this study, the authors compare the derivation of the pitch-angle diffusion coefficient using the non-linear simulations with the results deduced from the QLT. They find a good agreement between the two approaches at low and moderate turbulence levels but an increasing difference as the ratio of the turbulent to background field increases because of resonance broadening effects. The investigation of self-generated waves is one of the scientific objectives of the ongoing inner heliospheric missions solar orbiter (Zouganelis *et al.*, 2020) and Parker Solar Probe (Fox *et al.*, 2016).

One may also wonder if the heliospheric termination shock (TS) and the heliopause region (HR) may be sites of particle streaming effects. Thanks to the Voyager probes measurements, it appears that several observations could be investigated in such a framework: (1) highly anisotropic particle spikes have been identified by the Voyager 1 probe (Decker *et al.*, 2005) and (2) Voyager 2 observations show that anomalous cosmic rays (ACRs) (ACRs are produced in the inner heliosphere due to charge exchange between interstellar neutrals and the solar wind ions). These can have dynamical impacts on the TS dynamics (Richardson *et al.*, 2008). However, because the TS is almost perpendicular and the wave growth timescale should scale in the QLT as the cosine of the angle of the background magnetic field line with respect to the shock normal, it is hence not expected to have strong wave generation through the SI at the TS (le Roux and Webb, 2009). Nevertheless, dynamical effects of ACRs may necessitate to include more non-linear effects in the modeling of the particle transport. Indeed, Guo, Florinski, and Wang (2018) show that Voyager data are better reproduced with parallel diffusion coefficient values that can induce TS shock position variation. If the ACR pressure is high enough to induce such dynamical effects, then a modeling including a calculation of the diffusion coefficients levels may require the triggering of self-generated instabilities.

### B. Perspectives in stellar physics

One can note here that recent advances in the modeling of active stars can now provide us with refined views of their magnetosphere dynamics. These objects differ from our Sun in terms of activity level, magnetic field amplitude, rotation speed, stellar wind temperature, and stellar wind speed. Yet, their activity is likely connected with episodes of magnetic reconnection as for the Sun. Stellar EPs should be produced either during stellar flares or during CME expansions. EPs' acceleration and propagation modeling are important in the context of proto-star or proto-planetary accretion disk ionization (Fraschetti *et al.*, 2018) and for their direct impact over the atmosphere of planets (Fraschetti *et al.*, 2019). To our knowledge, no work transposing CME particle acceleration models to active magnetic dwarf or T Tauri stars has been yet proposed.

## V. THE STREAMING INSTABILITY IN ASTROPHYSICAL PLASMAS

The SI is very relevant in a large variety of astrophysical phenomena [see reviews by, e.g., Amato (2011); Marcowith *et al.* (2016, 2020)]. We will consider three main (the list is nonexhaustive) subjects where this instability is anticipated to play an important role in astrophysical plasma systems: (i) particle acceleration at collisionless shocks, (ii) instabilities driven by fast beams in jets issued from compact objects, and (iii) cosmic ray transport in the interstellar/galactic medium.

### A. Particle acceleration at collisionless shocks

Analytical and numerical studies presented above show that DSA is likely operating in fast shocks through the development of self-generated turbulence driven by strong anisotropy of the CR distribution imposed by the motion of the shock. Yet, it has been for long known that DSA on pre-existing background turbulence cannot produce CR energies much in excess to a few tens of TeV in supernova remnants (Lagage and Cesarsky, 1983). This aspect is really problematic as signal-to-noise ratios (SNRs) are anticipated to be the main CR sources in our Galaxy. Hence, the turbulence issued from the streaming instability seems absolutely essential to allow a longer particle confinement around the shock front and hence to allow more Fermi cycles to operate: stronger turbulence levels mean higher CR energies. The price to pay is that an extensive work is necessary to understand the time evolution of the instability and the process of turbulence generation. In addition to theoretical considerations, a strong magnetic field amplification seems to be required by the detection of x-ray filaments associated with synchrotron radiation produced by relativistic electrons in all historical supernova remnants (Parizot *et al.*, 2006). Finally, let us insist that the necessity to trigger self-generated turbulence at scales smaller than the particle gyroradius is also essential to unlock particle trajectories in relativistic shocks and to allow the Fermi acceleration cycles to proceed (Pelletier, Lemoine, and Marcowith, 2009). While progress has been made in the understanding of the injection process and how the CR distribution is built ahead of the shock front, still some issues remain to be clarified and are discussed below.

#### 1. The injection efficiency

Even if the process of injection is not directly connected with the ignition of the streaming instability, but a good knowledge of how CRs are injected into the DSA process is mandatory for a proper understanding of the longer term evolution of the shocked gas–CR system and then how high energy CRs can inject turbulence in this system through the interplay of the SI. The efficiency of particle injection into the Fermi acceleration process is an essential ingredient of the DSA process. Its amplitude *η* can only be evaluated using a microscopic approach because space and time scales involved in this process can be as short as the thermal electron scales. The injection efficiency can be defined in different ways. For simulation purposes, one can define it as the fraction of the bulk incoming energy flux converted into particles with energy larger than some threshold $Eth$. For instance, Caprioli and Spitkovsky (2014a) use $Eth=10\u2009m\u2009ush2/2$. These authors found that *η* is dependent on the shock speed $ush$, the magnetic field obliquity $\theta B$ (the angle between the shock normal direction with the upstream magnetic field), and the ambient plasma magnetization *σ* (the magnetization can be defined as the ratio of magnetic energy to kinetic flow energy).

A more complete survey of injection efficiency depending on the three above-mentioned parameters using the PIC simulations seems essential to conduct in complement to the results obtained with hybrid methods. Moreover, some obvious and in fact ongoing simulations (Crumley *et al.*, 2019) can extend PIC and hybrid simulations in electron–proton (e–p) plasma to multi-dimension for various magnetic field obliquity and shock speed including the mildly relativistic shock speed regime. To properly conduct these, injection efficiency deduce from PIC/hybrid simulations can be inserted as a recipe to account for the dependence with the local magnetic field obliquity.

Another important aspect is the time evolution of particle injection into the DSA process. The injection rate is, as we outlined below, dependent on the high-energy particle feedback as well as the time evolution of the shock. For instance, during late SNR evolution stages, which cover most of the object lifetime, the injection of pre-existing CRs is in competition with the thermal leakage process where non-thermal particles are directly injected from the thermal pool plasma (Cristofari and Blasi, 2019; Brose *et al.*, 2020).

#### 2. High-energy particle feedback

As acceleration proceeds, high-energy particles while triggering resonant and/or non-resonant perturbations can induce a corrugation of the shock front and hence modify the local shock magnetic field obliquity. This is a major non-linear feedback not yet covered in PIC and hybrid simulations [see the most extended simulations obtained in Haggerty and Caprioli, 2019] and still rarely in MHD studies (van Marle *et al.*, 2019). In addition, the perturbations and the turbulence triggered at large scales by the high-energy particles modify the propagation of less-energetic particles (Drury, 2011) [this is the so-called magnetic bootstrap process, see Blandford and Funk (2007)].

As the problem is clearly multi-scale, one should continue to develop combined PIC (Vlasov) and MHD approaches where large-scale shock dynamics is handled by PIC (Vlasov)-MHD codes, but in the meantime where some sub-spaces of the simulation box are studied using PIC simulations to evaluate particle injection dependence imposed by large-scale dynamics. A similar approach has been recently adopted in the context of magnetic reconnection in magnetic island coalescence (Makwana, Keppens, and Lapenta, 2018a).

#### 3. Maximum CR energy

A related subject is the way the maximum CR energy increases with time. This depends on the diffusion regime. Caprioli and Spitkovsky (2014b) find that for low-Mach-number shocks, the turbulent power spectrum $\delta B2(k)/B$ scales as $k\u22121$ in the shock precursor in agreement with expectation that the resonant SI dominates the generation of turbulence. Holcomb and Spitkovsky (2019) find that once in the relativistic regime the maximum CR energy $Emax(t)$ scales as t, a result consistent with particles being scattered with a Bohm diffusion coefficient as expected from the turbulence power law. At higher Mach numbers, a clear transition is seen and the turbulent spectrum appears softer as the dominant SI shifts from the resonant to the non-resonant branch. In this regime, the time dependence of $Emax$ is not yet been tested. It seems then important to extend the simulations to the high-Mach-number regime including 3D (van Marle *et al.*, 2019) and obliquity effects. The next step would be then to extract some general self-similar behavior for the shock and particle acceleration dynamics and evaluate the effective maximum energy expected in different CR sources.

#### 4. Ambient medium effects

The impact of the properties of the ISM of IGM over the development of the SI has still rarely been addressed until now. The ambient medium in which the shock front propagates is likely in-homogenous, with density, temperature, and magnetic stratification over scales comparable to the shock radius. Some analytical calculations using a WKB analysis have tried to handle these effects (Pelletier, Lemoine, and Marcowith, 2006), but more refined numerical setup including a turbulent ambient medium seems now accessible to numerical simulations. In that line of research, we can mention recent kinetic-MHD simulations (Inoue, 2019), which report on the production of non-resonant modes during supernova remnant (SNR) shock–molecular cloud interaction produced by the current of CR having escaped the remnant shock. In addition, partially ionized media require the inclusion of a neutral component in a multi-fluid kinetic model (Morlino *et al.*, 2013). The inclusion of a neutral component in the above-mentioned techniques still needs to be operational. Finally, the inclusion of a pre-existing turbulent spectrum generated from large-scale simulations should be of interest to mimic the impact of perturbations triggered by more energetic CRs (see Sec. V A 2).

#### 5. Streaming instability studies in other astrophysical objects

Another promising line of research is the investigation of the generation of SI in different astrophysical objects such as supernovae (Marcowith *et al.*, 2018), novae (Martin and Dubus, 2013), or young stellar jets (Padovani *et al.*, 2016; Araudo, Padovani, and Marcowith, 2021) in order to derive a more refined estimation of the expected magnetic field strength at these fast shocks as well as time-dependent and multi-zone non-thermal emission models. In particular, it should also be possible to investigate the production and propagation of high-energy particles in the vicinity of high mass binary stars. The collision regions between the stellar winds of hot massive stars contain high-Mach shocks with potentially high compression rates (depending on whether the shock is adiabatic or radiative [e.g., Stevens, Blondin, and Pollock, 1992; Pittard, 2009; van Marle *et al.*, 2011]) and evidence for the production of high-energy particles in such an environment has been observed (Abdo *et al.*, 2010). Dubus, Lamberts, and Fromang (2015); Balbo and Walter (2017) combined a hydrodynamical model of the colliding wind shocks in the *η*-Carinae system with an *ad hoc* description for the DSA to predict the time-dependent nature of the particle acceleration process as a function of the orbital position of the two stars. Using current numerical methods, this model could be extended to include the influence of the non-thermal particles on the local shock conditions.

One can also expect the triggering of the SI in relativistic flows. This aspect is discussed in Sec. V B.

### B. Fast beam-driven instabilities in compact objects

Relativistic non-thermal particle beams or jets are expected to be the rule in flows issued from compact objects, that is, in active galactic nuclei, gamma-ray bursts, x-ray binaries, or in pulsar winds. Microphysics is often included in a phenomenological way in the models where non-thermal particles and magnetic turbulence energy densities are parametrized as fractions of the jet/wind kinetic power. Going toward more microphysics requires to include a series of complex processes in source dynamics. A first essential physical ingredient in compact source emitting gamma-ray photons is the photon–photon electron–positron pair process. The photon–photon pair creation threshold is satisfied if the energy of the two photons exceeds the electron mass energy in the rest frame of the incoming electron, a criterion in practice easily satisfied in compact objects. It seems hence not unrealistic to include an electron–positron component into the matter content of the relativistic beam or jets ejected from compact sources. Another unavoidable aspect of jet/beam physics is that local radiation fields are intense. Non-thermal electron–positron pairs released from the central object are then subject to strong synchrotron, Inverse Compton, or in the densest regions Bremsstrahlung losses. One important aspect is first to account for these radiative losses in the time evolution of the non-thermal particle distribution (Saugé and Henri, 2004). Proper dynamics of the non-thermal component will have to include these radiation effects at some stage. Still, the way the particles are re-accelerated along the jet (by Fermi-like processes or magnetic reconnection) is a missing piece of the models. A possible way to ensure this re-acceleration is to trigger a streaming instability due to the relative motion of the pair beam with respect to the background electron–proton plasma. We discuss below two astrophysical configurations in which the SI may supply the necessary amount of free energy.

#### 1. Jets from compact objects

Jets in compact objects are structured in layers, composed of slower outflows pervaded by fast jets. These multi-layers outflows are modeled in the two-flow or spine-sheath frameworks (Sol, Pelletier, and Asseo, 1989; Henri and Pelletier, 1991; Ghisellini, Tavecchio, and Chiaberge, 2005). The shearing motion between the central beam with respect to the background slower plasma must trigger a series of instabilities (Marcowith *et al.*, 1997) among which the SI should be one of the most important.

PIC simulations of mildly relativistic or relativistic electron–positron beams have started to address the physics of beam/background interaction, the setup of shocks, and the development of beam-driven instabilities (Sironi and Giannios, 2014; Bret and Dieckmann, 2017). Of particular relevance for non-thermal radiation from jets is the case of a pair beam drifting inside a background plasma composed of electrons and protons (Dieckmann *et al.*, 2019; Dieckmann, 2020). In these simulations, the pair plasma has a temperature of 400 keV, while the background electrons and protons are cold with a temperature of 2 keV. The simulations show that the central beam is rapidly separated from the background plasma by the formation of an electromagnetic piston. The protons of the background plasma are entrained at non-relativistic speeds and heat by a filamentation instability, which develops with the positrons. Entrained protons form a cocoon around the pair beam. The cocoon is separated from the background plasma by a fast magnetosonic shock.

The next steps will have to focus more closely on non-thermal (with temperatures in excess to the electron rest mass) particle dynamics. First, using hybrid or PIC-MHD techniques it is relevant to investigate beam-driven instabilities and the subsequent feedback of self-generated turbulence on the particle acceleration cycle. It would be interesting to test to which extend SI can be triggered in case the beam is composed of relativistic particles and evaluate its dominance over other types of kinetic/convective instabilities. Then, MHD simulations will be necessary to capture the structure of the shearing flows. Sub-grid models can be included to mimic the effect of the kinetic physics (self-generated turbulence, radiation).

#### 2. Gamma-ray halos around pulsars

Around pulsars recent TeV gamma-ray observations by the HAWC (High-Altitude Water Cherenkov) on-ground gamma-ray observatory (see https://www.hawc-observatory.org/) show extended TeV halos likely composed of energetic electron–positron pairs (Abeysekara *et al.*, 2017). From the extension of these halos, it is possible to deduce using a simple diffusion model that particles propagate from the nebula in the interstellar medium with a reduced diffusivity with respect to ISM standards. The origin of this suppression of the particle's mean free path is not clear yet: it can be due to self-generated turbulence triggered by escaping electron–positron pairs through the triggering of the SI (Evoli, Linden, and Morlino, 2018) or to specific properties of the background (ISM or associated with a SNR) turbulence (López-Coto and Giacinti, 2018).

One key question if we want to explain reduced particle diffusivity in pulsar gamma-ray halo by invoking the SI is the effective ability of electron–positron pairs to trigger it (Fang, Bi, and Yin, 2019). If, in the context of pulsar halos, as it is debated, pure electron–positron beams do not carry enough energy density to trigger a SI (Fang, Bi, and Yin, 2019), then the addition of a ionic component may allow it. Acceleration of ions in the magnetosphere of compact objects is also a potential way to explain ultrahigh energy CRs (Guépin, Cerutti, and Kotera, 2020).

### C. Cosmic ray transport in the interstellar/galactic medium

While trying to escape from their acceleration sites, the most energetic CRs can trigger SI and then turbulence in the interstellar medium surrounding a CR source. This self-generated turbulence can confine the particles and have an impact over several CR propagation-related issues. A first question is how much should we expect of source contribution to the CR grammage calculation (D'Angelo, Blasi, and Amato, 2016)? We remind here that the grammage is the quantity of matter intercepted by a CR during its journey in the galaxy. This quantity is inferred from direct observations of secondary to primary element abundances in the CR spectrum. A second problem already mentioned above is linked with the recent detection of gamma-ray halos around Geminga and Monogem pulsars (see Sec. V B). Should we expect with the advent of new gamma-ray facilities like the Cherenkov telescope array (CTA) (see https://www.cta-observatory.org/) or the Large High Altitude Air Shower Observatory (LHAASO) (see http://english.ihep.cas.cn/lhaaso/) to discover more of these gamma-ray halos around other pulsars, around some SNRs (Malkov *et al.*, 2013; D'Angelo *et al.*, 2018; Inoue, 2019; Nava *et al.*, 2019) or even massive star clusters (Bykov *et al.*, 2020)? This second question could be rephrased as follows: How do CRs escape from their sources and how much is their propagation close to their sources different from the propagation in the large-scale ISM?

In effect, the anisotropy produced by CRs escaping from a source conducts to the generation of streaming modes. Depending on the intensity of the escaping current, the non-resonant regime may even take over from the resonant regime for CRs propagating in the ISM (Zweibel and Everett, 2010; Schroer *et al.*, 2021) or CR escaping the precursor of SNR in interaction with a molecular cloud (Inoue, 2019). For CRs to escape from SNRs, the ability of the resonant SI to trigger a sufficient amount of turbulence to reduce the particle diffusion coefficient depends strongly on the properties of the surrounding ISM and on the CR energy regime (Nava *et al.*, 2016, 2019). In hot, fully ionized ISM, the CR self-generated turbulence is damped mostly because of its interaction with the background turbulence generated at large galactic scales. Ion–neutral collisions reduce the development of the resonant SI in partially ionized media below TeV energies. Moreover, resonant waves produced by multi-TeV particles are less affected because ion and neutral motions are coupled at low frequencies (Brahimi, Marcowith, and Ptuskin, 2020). It should be of great interest to investigate using either kinetic or fluid approaches the back-reaction effects of the CRs and the self-generated turbulence over the different ISM phases. Not only SNR are anticipated to be CR sources and it would be very valuable to redo these studies in the environment of massive star clusters and OB star associations (Bykov *et al.*, 2020). The analysis of the SI triggering can be based on similar setup as proposed by Lucek and Bell(2000), where some prescribed anisotropic CR distribution can be injected into the ISM including or not a specific model of background turbulence. This setup should help to constrain the way CR can escape from their sources.

In addition to the way CRs escape from their sources, the way CRs propagate in the different phases of interstellar gas during their journey toward Earth is still poorly known. For instance, of particular relevance is the impact of propagation of CRs over the ionization rate at different density columns in molecular clouds (Padovani *et al.*, 2020), so at different depth in the molecular medium. The ionization rate controls the coupling between gas and magnetic fields and hence is an essential control parameter of the molecular cloud collapse process. Beyond this, it seems important to investigate the process of collapse depending if the cloud is isolated or under the influence of a CR source. The resonant SI appears as one of the main instability (but likely not the only one), which controls the penetration of CRs in the molecular cloud (Cesarsky and Volk, 1978; Morlino and Gabici, 2015; Ivlev *et al.*, 2018; Phan, Morlino, and Gabici, 2018; Silsbee and Ivlev, 2019; Phan *et al.*, 2020; Silsbee and Ivlev, 2020).

ISM dynamics is strongly controlled by a series of feedback at all scales (Chevance *et al.*, 2020). CRs contribute to this feedback by their pressure (or the gradient of their pressure) and, as mentioned above by ionization. At large galactic scales, CRs escaping the galactic disk can drive SI instabilities and contribute to the launching of galactic winds (Recchia, Blasi, and Morlino, 2017; Girichidis *et al.*, 2018; Blasi and Amato, 2019) and the acceleration of cold clouds in the intra-galaxy cluster medium (Wiener, Zweibel, and Ruszkowski, 2019). Depending on their diffusivity, CRs can modify the development of the thermal instability in the ISM (Shadmehri, 2009; Commerçon, Marcowith, and Dubois, 2019) or in halos (Kempski and Quataert, 2020), this diffusivity is expected to be reduced in the case of triggering of the SI. The pressure supported by CRs is a source of feedback by producing a buoyancy of galactic magnetic field lines, which participates to the galactic dynamo, and this is the Parker instability (Parker, 1992). Only until a recent study, the effect of CR streaming has been included in the Parker instability analysis (Heintz and Zweibel, 2018).

Finally, the self-generated turbulence produced by CRs has another implication. It can be damped in the interstellar (or intergalactic) medium basically in two ways: first, it can be partly used to reaccelerate sub-relativistic CRs, and second, it can be damped in collisional (like ion–neutral collisions) or collisionless (like Landau damping) processes. The latter contributes to heat the ISM (Zweibel, 2020) or the IGM (Ruszkowski, Yang, and Reynolds, 2017). Hence, it appears important to include such wave generation aside the contribution of CRs (Thomas and Pfrommer, 2019; Hopkins, Squire, and Butsky, 2021) and aside the contribution of large-scale injected turbulence (Drury and Strong, 2017) in dynamical studies.

Hydrodynamical or magnetohydrodynamical models calculating ISM dynamics from the scales involved in cloud collapse to scales involved in galactic winds have to include at some terms the physics of the SI by the addition of a CR component. It seems important to properly include the effect of SI into fluid simulations including heating/cooling and thermal conduction. Treating a kinetic instability triggered by low-energy CRs with Larmor radii smaller than typical grid resolution scales requires in this formalism to do some sub-grid physics (Thomas and Pfrommer, 2019; Hopkins, Squire, and Butsky, 2021) or to use mesh-free simulations (Hopkins, 2017; Hopkins *et al.*, 2021; Thomas, Pfrommer, and Pakmor, 2021) if one wants to catch large-scale galactic dynamics. Inclusion of the physics of multi-CR components in fluid models can help to include some kinetic effects, in particular, the impact of the different turbulence scales over CR propagation (Girichidis *et al.*, 2020).

## VI. CONCLUSION

CRs due to their high carried pressure and/or strong anisotropy produced by their streaming are able to trigger different types of kinetic and convective instabilities. In this article, we only focus on the different regimes of the streaming instability. This instability has a resonant branch and a non-resonant branch, which have been both subjects of strong research activities during recent years in various scientific contexts: CR acceleration at space plasma and astrophysical shocks, and CR propagation from their sources and in the interstellar/intergalactic gas. As these instabilities can be triggered by GeV or multi-GeV CRs, they should also participate to a feedback process induced by CRs over multi-scale dynamical processes in our galaxy, from molecular cloud collapse to the production of galactic winds or in galaxy clusters. The inclusion of CR streaming and the damping of self-generated turbulence as strong heating process into fluid models seem important to account for the feedback processes necessary to properly account for the galactic star formation rate. This should pass by theoretical and numerical studies of the impact of CR streaming over the thermal instability and the Parker instability. In another context, the study of the dynamics of relativistic electron–positron beams in jets issued from compact objects is also an opportunity to investigate the possible role of the SI in mediating turbulence production and particle stochastic acceleration. Microscopic kinetic processes associated with non-thermal components are now entering into the game of the multi-scale physics of astrophysical complex systems.

## ACKNOWLEDGMENTS

This work is supported by the ANR-19-CE31-0014 GAMALO project. The authors thank the referees for their valuable comments which helped to improve the manuscript substantially.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.