# Kinetic simulations of electron–positron induced streaming instability in the context of gamma-ray halos around pulsars

Illya Plotnikov<fnref n="1" />, Allard Jan van Marle<fnref n="2" />, Claire Guépin<fnref n="2" />, Alexandre Marcowith<fnref n="2" />, and Pierrick Martin<fnref n="1" />

Received 19 February 2024 / Accepted 6 June 2024

## ABSTRACT

**Context.** The possibility of slow diffusion regions as the origin for extended TeV emission halos around some pulsars (such as PSR J0633+1746 and PSR B0656+14) challenges the standard scaling of the electron diffusion coefficient in the interstellar medium.

**Aims.** Self-generated turbulence by electron–positron pairs streaming out of the pulsar wind nebula was proposed as a possible mechanism to produce the enhanced turbulence required to explain the morphology and brightness of these TeV halos.

**Methods.** We perform fully kinetic 1D3V particle-in-cell simulations of this instability, considering the case where streaming electrons and positrons have the same density. This implies purely resonant instability as the beam does not carry any current.

**Results.** We compare the linear phase of the instability with analytical theory and find very reasonable agreement. The non-linear phase of the instability is also studied, which reveals that the intensity of saturated waves is consistent with a momentum exchange criterion between a decelerating beam and growing magnetic waves. With the adopted parameters, the instability-driven wavemodes cover both the Alfvénic (fluid) and kinetic scales. The spectrum of the produced waves is non-symmetric, with left-handed circular polarisation waves being strongly damped when entering the ion-cyclotron branch, while right-handed waves are suppressed at smaller wavelength when entering the Whistler branch. The low-wavenumber part of the spectrum remains symmetric when in the Alfvénic branch. As a result, positrons behave dynamically differently compared to electrons. The final drift velocity of positrons can maintain a larger value than the ambient Alfvén speed $V_A$ while the drift of electrons can drop below $V_A$. We also observed a second harmonic plasma emission in the wave spectrum. An MHD-PIC approach is warranted to probe hotter beams and investigate the Alfvén branch physics. We provide a few such test simulations to support this assertion.

**Conclusions.** This work confirms that the self-confinement scenario develops essentially according to analytical expectations, but some of the adopted approximations (like the distribution of non-thermal particles in the beam) need to be revised and other complementary numerical techniques should be used to get closer to more realistic configuration.

**Key words.** instabilities – plasmas – methods: numerical – stars: neutron

## 1. Introduction

The recent detection of gamma-ray halos around PSR J0633+1746 and PSR B0656+14 (the Geminga pulsar and the pulsar in the Monogem ring) by the HAWC Collaboration (Abeysekara et al. 2017a,b) demonstrates that the propagation of high-energy cosmic rays around their sources can be strongly modified, involving a reduced diffusivity of the particles by more than two orders of magnitude with respect to what is expected from direct measurements of secondary to primary cosmic-ray nuclei ratios. The presence of a gamma-ray halo around Geminga has been since confirmed by Fermi-LAT in the GeV range (Di Mauro et al. 2019). Again, at these energies, a reduction of the diffusion coefficient by two orders of magnitude seems necessary to explain the halo extension. Di Mauro et al. (2021) confirm this trend using Fermi data analysis for three supplementary sources from the HAWC Collaboration catalogue. More recently, the LHAASO collaboration reported on a similar structure around the middle-aged pulsar J0622+3749 in the 10–100 TeV range (Aharonian et al. 2021). Finally, the H.E.S.S Collaboration reported on a detection of a gamma-ray halo around the Geminga pulsar with similar properties deduced from HAWC data (H.E.S.S. Collaboration 2023).

As pulsars are anticipated to be sources of cosmic-ray leptons (electrons and positrons), the low particle diffusivity at the origin of gamma-ray halos can have several strong implications. Firstly, the suppression of the diffusion around Geminga and Monogem has called into question the contribution of these sources to the lepton flux detected on Earth (Profumo et al. 2018; Fang et al. 2018). It appears that considering the age of these pulsars (about 340 kyrs for Geminga and 110 kyrs for Monogem) the contribution of the confinement zone could provide a local enhancement of the lepton flux which can – at least – explain a fraction of the positron excess observed by the AMS-02 (Aguilar et al. 2014) and the PAMELA (Adriani et al. 2009) experiments (Fang et al. 2018; Tang & Piran 2019). The same trend is found by Di Mauro et al. (2019) considering a time-dependent release of positrons and accounting for pulsar motion. The authors find that Geminga/ Monogem pulsars can at most contribute to 20%/a few % of the positron flux at high energy detected on Earth. Secondly, if the pulsar gamma-ray halo is a general phenomenon, then the cosmic ray content and gamma-ray flux are expected to be enhanced in the galactic disc following the source distribution (Jóhannesson et al. 2019). Schroer et al. (2023) invoke the possible effect due to radiation losses to produce softer

<footnote n="1">¹ IRAP, CNRS, OMP, Université de Toulouse III – Paul Sabatier, Toulouse, France
e-mail: illya plotnikov@irap.omp.eu</footnote>

<footnote n="2">² Laboratoire Univers et Particules de Montpellier (LUPM), Université de Montpellier, CNRS/IN2P3, CC72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France</footnote>
spectra for electron–positron pairs escaping a gamma-ray halo. The loss-softening effect seems necessary to explain the positron spectrum and fraction observed on Earth (Evoli et al. 2021). Manconi et al. (2020), using a two-zone model combined with the ATNF pulsar distribution, can explain AMS-02 lepton data with a reasonable fraction of pulsar spin-down luminosity converted into electron–positron pairs of about 10%. However, this possibility is questioned by Martin et al. (2022a). These authors also used a two-zone diffusion model and show that assuming that gamma-ray halos are common around middle-aged pulsars may lead to overshooting the local positron flux measured with the AMS-02 experiment. Martin et al. (2022a) estimate that at most a fraction of 5–10% of nearby middle-aged pulsars could develop a gamma-ray halo. In total, in the Galaxy the upcom-ing Cherenkov Telescope Array is estimated to be able to detect about 100 gamma-ray halos (Martin et al. 2022b).

The origin of the slow particle diffusion remains a subject of active debate. Several interpretations of the HAWC observations have been put forward so far. Recchia et al. (2021) interpreted the HAWC gamma-ray intensity profiles as a combination of bal-listic and diffusive particle motions, in a model that does not invoke any specific reduced diffusion to explain the gamma-ray data. Yet, accounting for the observed level of gamma-ray emis-sion requires a power input that exceeds the available spin-down luminosity in some cases (Bao et al. 2022). Alternatively, per-perpendicular diffusion could be the origin of the slow diffusion, and favourable magnetic field line orientation with respect to the line of sight produces the observed halos, but the scenario seems more and more unlikely as we discover a growing num-ber of such objects (De La Torre Luque et al. 2022). López-Coto & Giacinti (2018) argue that the reduced diffusivity is due to specific properties of the magnetic turbulence surrounding the pulsars. In particular, the gamma-ray data are well-fitted if the turbulence coherence length is ~1 pc, hence reduced by one or two orders of magnitude with respect to its expected value in the interstellar medium, as deduced from Faraday rotation and depolarisation data (Haverkorn et al. 2008). One possible ori gin of this reduced coherence length is that the pulsar still finds itself in a medium under the influence of the associated super-nova remnant turbulence (Fang et al. 2019). Finally, Evoli et al. (2018) investigate the triggering of the resonant streaming insta-bility by an electron–positron ($e^–e^+$) beam injected at the pulsar wind termination shock. The authors considered the injection of an $e^–e^+$ power-law momentum distribution as $p^{–3.5/-3.2}$ from the pulsar termination shock and solved a system of coupled kinetic equations for cosmic rays and triggered waves includ-ing turbulent wave damping. They find a reduced diffusivity within a radius of about 10–20 pc around the source for time scales of the order of 100 kyr. These results were revised (and corrected) by Mukhopadhyay & Linden (2022), who showed that a reduced diffusion region can be sustained over larger dis-tances and longer duration, encompassing the properties of the Geminga system, even when accounting for ion-neutral wave damping and assuming a flatter particle injection spectrum.

The latter work reveals that the self-confinement scenario relies on several assumptions for it to be relevant. Essentially, it requires that maximum power is channelled into the pair stream to provide a growth of the turbulence up to a sufficient level. This implies: (i) a beam starting right after the pulsar's birth; (ii) one-dimensional injection along a flux tube; (iii) no proper motion of the pulsar; and (iv) an injection spectrum peaking around 10–100 TeV, the energy where reduced diffusion is the most solidly constrained. In realistic conditions, it is not guaranteed that these favourable requirements are met: for (i), injection initially occurs

in the pulsar wind nebula, and pairs may not be immediately available to the ambient medium until the system has evolved to a bow-shock nebula stage, about several thousands or tens of thou-sands of years later, at which time the pulsar has lost much of its spin-down power (Gaensler & Slane 2006); for (ii) and (iii), real-istic conditions such as three-dimensional expansion and proper motion of the pulsar imply beam dilution, hence less power for exciting the instability; and for (iv), particle spectra produced by pulsars and their nebulae are observationally constrained to peak around 100 GeV (Torres et al. 2014).

Nevertheless, Mukhopadhyay & Linden (2022) show that the self-confinement option remains a viable explanation in some conditions, all the more so that it can be combined with other scenarios (pre-existing fluid turbulence imparted by the expan-sion of the remnant<fnref n="1" />, or kinetic turbulence triggered by the acceleration of protons at its forward shock). For that reason, we consider that the underlying physics requires a more detailed understanding.

In this work, we investigate the case of electron–positron driven instabilities in a similar framework as proposed by Evoli et al. (2018) but we extend it to the case of higher electron– positron to background plasma density ratios, thus discarding the test-particle assumption<fnref n="2" />. This regime could be applicable to the instability development in Pulsar wind nebula during early evolutionary phase when the density of released pairs is the high-est. Our study concentrates on the microphysics of both linear and non-linear phases of the pair-driven streaming instability by means of kinetic particle-in-cell (PIC) numerical simula-tions. Contrary to a recent study by Gupta et al. (2021) we consider only the case where the density of streaming electrons and positrons is equal. Thus, driving only resonant instability and discarding faster-growing non-resonant instability. The main objective of this work is to explore the microphysics of the self-confinement of a pair stream (general review on the streaming instability in astrophysical plasmas can be found in Marcowith et al. 2021). Exporting that to a realistic model of a pulsar halo is deferred to future studies.

The layout of this paper is as follows: in Sect. 2 we derive the general expression for the linear growth rate and give estimates in the case of a Maxwell–Jüttner distribution. We also give an estimate of the magnetic field saturation level in the case of the resonant streaming instability. In Sect. 3 we present the results of PIC simulations. We further provide some interpretation. Section 4 proposes a discussion before concluding.

## 2. Streaming instability linear growth rates

Evoli et al. (2018) consider a particular case of a beam com-posed of electrons and positrons. Here, we do not first assume any particular composition of the drifting beam. This allows us to present below the general formalism necessary to derive the growth rate of the streaming instability. This physical set-up is different from the case of a mildly-relativistic electron– positron spine injected into a magnetised electron-proton plasma

<footnote n="1">¹ Pre-existing turbulence could help because it may provide a supplementary source of particle scattering if acting on an appropriate scale. Yet, this effect has to be assessed with a dedicated study.</footnote>

<footnote n="2">² In this work, we use the test-particle terminology following Evoli et al. (2018). It is valid as long as the electron–positron beam is not modifying the dispersion relation of Alfvén waves, that is when the amplitude of excited waves is small: $\delta B \ll B_0$. When this limit is not satisfied, the instability drives either strongly modified Alfvén waves or waves or different nature and scale (e.g., kinetic Alfvén waves, Whistler, ion-cyclotron, etc.).</footnote>
<table><caption>**Table 1.** Notations used in this work.</caption><thead><tr><th>Quantity</th><th>Notation</th></tr></thead><tbody><tr><td>Light speed, elementary charge</td><td>c, e</td></tr><tr><td>Frequency, wavenumber</td><td>ω, k</td></tr><tr><td>Magnetic field strength</td><td>B</td></tr><tr><td>Particle species<br>(electron, proton, CRs)</td><td>α (e, i, cr)</td></tr><tr><td>Mass, charge, momentum</td><td>m<sub>α</sub>, q<sub>α</sub>, p<sub>α</sub></td></tr><tr><td>Lorentz factor</td><td>γ<sub>α</sub> = √(1 + (p<sub>α</sub>/m<sub>α</sub>c)<sup>2</sup>)</td></tr><tr><td>Gyrofrequency</td><td>Ω<sub>α</sub> = q<sub>α</sub>B/γ<sub>α</sub>m<sub>α</sub>c = Ω<sub>cα</sub>/γ<sub>α</sub></td></tr><tr><td>Gyroradius</td><td>R<sub>L,α</sub> = c/Ω<sub>cα</sub></td></tr><tr><td>Background density</td><td>n<sub>i</sub> = n<sub>e</sub></td></tr><tr><td>Beam density</td><td>n<sub>cr</sub></td></tr><tr><td>Alfvén velocity</td><td>V<sub>A</sub> = B/√(4πn<sub>i</sub>m<sub>i</sub>)</td></tr><tr><td>Drift velocity</td><td>V<sub>D</sub></td></tr><tr><td>Electron/ion temperature<br>(background)</td><td>T<sub>e,i</sub> = θ<sub>bg</sub>m<sub>e</sub>c<sup>2</sup>/k<sub>B</sub></td></tr><tr><td>Electron plasma frequency</td><td>ω<sub>pe</sub> = √(4πn<sub>e</sub>e<sup>2</sup>/m<sub>e</sub>)</td></tr><tr><td>Electron skin-depth</td><td>d<sub>e</sub> = c/ω<sub>pe</sub></td></tr></tbody></table>

(Dieckmann et al. 2018, 2019, 2022). A situation where our formalism does not apply, because both plasmas are initially separated by hydrodynamical structures. The formalism developed below could however be of interest, but at later stages if some mixing applies between the two plasma components. Hence, potentially, it could also be applied to compact object jets.

## 2.1. General form of the dispersion relation

We adopted CGS units and all notations used in this work are reported in Table 1. The dispersion relation for modes of frequency ω and wavenumber k can be written as (Amato & Blasi 2009):

$$ \frac{k^2 c^2}{\omega^2} = 1 + \chi. \qquad (1) $$

The plasma response function (or susceptibility) χ has two components, one for the background χ<sub>bg</sub> and one for the CR beam χ<sub>cr</sub>. We conduct the calculation in the beam rest frame in which the CR distribution is assumed to be isotropic. Each beam component is assumed to drift at the same speed V<sub>D</sub> with respect to the background plasma assumed to be cold. We first calculate the CR contribution χ<sub>cr</sub>. Assuming an isotropic distribution in the momentum space g(p) given by f<sub>cr</sub> = n<sub>cr</sub>g(p)/4π we have

$$ \chi_{\text{cr}} = \sum_{\alpha} \frac{\pi q_{\alpha}^2}{\omega} n_{\text{cr},\alpha} \int_{p_0}^{p_{\text{max}}} dp \, p^2 \nu(p) \frac{dg_{\alpha}(p)}{dp} \\ \times \int_{-1}^{1} \frac{(1 - \mu^2)}{\omega + k\nu(p)\mu \pm \Omega_{\alpha}} d\mu, \qquad (2) $$

where μ is the particle pitch-angle cosine. We only consider the fastest growing modes propagating parallel to the background large-scale magnetic field, namely we set k<sub>‖</sub> = k in the above expression. Each particle species has a distribution g<sub>α</sub>(p) such that ∫ dpp<sup>2</sup>g<sub>α</sub>(p) = 1; different forms of particle distribution are discussed below. Finally, n<sub>cr,α</sub> is the CR density. If the calculation is carried in the background plasma frame, a supplementary term (1 − kV<sub>D</sub>/ω) must come in front of Eq. (2). This term is necessary to quench the instability if the drift speed matches the

phase speed of the destabilised modes. For a particular species α the integral over μ can be derived using the Plemelj formula:

$$ \int_{-1}^{1} \frac{(1 - \mu^2)}{\omega + k\nu\mu \pm \Omega_{\alpha}} d\mu = \mathcal{P} \int_{-1}^{1} \frac{(1 - \mu^2)}{\omega + k\nu\mu \pm \Omega_{\alpha}} d\mu \\ -i\pi \int_{-1}^{1} (1 - \mu^2) \delta(\omega + k\nu\mu \pm \Omega_{\alpha}) d\mu, \qquad (3) $$

where the particle speed is a function of the momentum p and can be written as v(p). We find:

$$ \chi_{\text{cr}} = \frac{\pi e^2}{k\omega} \sum_{\alpha} N_{\text{CR},\alpha} \int_{p_0}^{p_{\text{max}}} dp \, p^2 \frac{dg_{\alpha}(p)}{dp} \\ \left[ \left(1 - \mu_{\text{R},\alpha}^2\right) \log \left| \frac{1 - \mu_{\text{R},\alpha}^{-1}}{1 + \mu_{\text{R},\alpha}^{-1}} \right| - 2\mu_{\text{R},\alpha} - i\pi \left(1 - \mu_{\text{R},\alpha}^2\right) \right], \qquad (4) $$

where the wave-particle resonance selects a particular pitch-angle cosine

$$ \mu_{\text{R},\alpha} = \mp \frac{\Omega_{\alpha}}{k\nu} - \frac{\omega}{k\nu}. \qquad (5) $$

We have at this stage considered that the beam is composed of particles with charges q = ±e. Equation (4) has three terms under the p-integral: the two first terms are connected to the CR current and are associated with the non-resonant branch of the streaming instability, the last term is associated with the resonant branch.

The background electron-proton (e<sup>-</sup>p<sup>+</sup>) plasma contribution is given by Eq. (14) in Evoli et al. (2018). Considering that each beam component is cold and drifts at the same speed V<sub>D</sub> we have:

$$ \chi_{\text{bg}} = - \frac{4\pi e^2}{\omega^2} \left( \frac{n_i}{m_i} \frac{\omega + kV_D}{\omega + kV_D \pm \Omega_{\text{ci}}} + \frac{n_e}{m_e} \frac{\omega + kV_D}{\omega + kV_D \pm \Omega_{\text{ce}}} \right). \qquad (6) $$

It is easy to add any ion component (e.g., Helium) to this equation. The background is considered as cold enough that the thermal motion does not induce any modification of the linear growth rate. In the interstellar medium (ISM) this would correspond to warm or cold media. For pulsar propagation into hotter media, like superbubbles, thermal background effects are expected to alter the growth of oblique modes (Foote & Kulsrud 1979). The non-resonant branch of the instability is also modified in that case (Zweibel & Everett 2010; Marret et al. 2021).

## 2.2. Low frequency limits

We now assume that (ω + kV<sub>D</sub>) ≪ Ω<sub>ci</sub> ≪ Ω<sub>ce</sub>. Hence, if we have n<sub>e</sub> = n<sub>i</sub>, the contribution from the background e<sup>-</sup>p<sup>+</sup> plasma is

$$ \chi_{\text{bg}} \approx \frac{(\omega + kV_D)^2 c^2}{\omega^2 V_A^2}. \qquad (7) $$

For cosmic rays contribution, we need to specify the distribution function g(p). In this work we considered a Maxwell–Jüttner (MJ) distribution unless otherwise specified. The results for the power-law and mono-energetic cases are reported in Appendix B. The solutions in the case of a MJ distribution are used in the numerical experiments in Sect. 3. We restrict ourselves to the case of a pure electron–positron beam with the same density for each species. More complex cases involving different electron and positron densities or momenta and the addition of a hadronic component deserve future works.
We assume a beam composed of electrons and positrons with $n_{cr,+} = n_{cr,-} = n_{cr}$ (and $g_+ = g_- = g$) and $\omega \ll \Omega_{e,+}$. The non-resonant term vanishes and we recover the results by Evoli et al. (2018) (their Eq. (17)), namely:

$$ \chi_{cr} = -2 \frac{i\pi^2 e^2}{k\omega} n_{cr} \int_{p_{\text{min}}}^{p_{\text{max}}} dp \frac{dg(p)}{dp} (p^2 - p_{\text{min}}(k)^2), \quad (8) $$

where the maximum particle momentum is $p_{\text{max}}$ and $p_{\text{res}}(k) = \frac{m_e|\Omega_{\text{ce}}|}{k}$. The factor 2 in front of the RHS part of Eq. (8) comes out as we have $n_{cr,+} + n_{cr,-} = 2n_{cr}$. For simplicity and tractability, we adopt the ultra-relativistic limit $E \approx pc \gg m_e c^2$. The MJ distribution is then defined as:

$$ g(p) = \frac{1}{m_e^3 c^3} \frac{1}{\theta K_2(\frac{1}{\theta})} e^{-\gamma(p)/\theta}, \quad (9) $$

where $\theta = k_B T_{cr}/(m_e c^2)$, $k_B$ is the Boltzmann constant, $T_{cr}$ is the temperature of CR electrons (same for CR positrons), and $K_2$ is the modified Bessel function of the second kind. The particle Lorentz factor is $\gamma(p)$. The distribution for positrons is the same with a density $n_{cr,+}$.

The choice of MJ distribution instead of a power-law is motivated by three main aspects. The first one is that MJ is physically motivated and allows to restrict the range of particle Larmor radii, hence to better isolate the resonant modes in our simulations. The second one is that for broad distributions extending to very high energies – such as power-law distributions – excessive numerical resources are required with the PIC method to cover sufficient dynamical scales in energy or momentum. The third reason is that MJ distribution is analytically tractable. This allows direct comparisons of the growth rate with numerical results (control of spurious non-physical effects).

Evaluating the integral in Eq. (8) we find

$$ \begin{gathered} \chi_{cr} \approx i \frac{4\pi^2 e^2}{\omega k} \frac{n_{cr}}{m_e c} \frac{\theta}{K_2(1/\theta)} I(k, \gamma_{\text{min}}), \\ I(k, \gamma_{\text{min}}(k)) = \exp(-\frac{\gamma_{\text{min}}}{\theta}) \left(\frac{\gamma_{\text{min}}(k)}{\theta} + 1\right), \\ \gamma_{\text{min}}(k) = \sqrt{1 + \left(\frac{\Omega_{\text{ce}}}{kc}\right)^2}, \end{gathered} \quad (10) $$

and $\gamma_{\text{min}} = \gamma(p_{\text{max}})$. We have assumed that $\gamma(p_{\text{max}}) \gg \theta$. The dispersion relation then reads as

$$ k^2 V_A^2 = \omega^2 \left(\frac{V_A}{c}\right)^2 + (\omega + kV_D)^2 + i\pi \frac{\omega \Omega_{\text{ci}} n_{\text{cr}}}{k R_{\text{L,e}}} \frac{\theta}{n_i K_2(1/\theta)} I(k, \gamma_{\text{min}}), \quad (11) $$

where $R_{L,e} = c/\Omega_{\text{ce}} = R_{L,cr}/\theta$ and $R_{L,cr} = k_B T_{cr}/eB$. The instability growth rate is given by the system

$$ k^2 V_A^2 \approx \tilde{\omega}^2 + i\omega \frac{n_{cr}}{n_i} \Omega_{ci} I_2(k), \quad (12) $$

with $\tilde{\omega} = \omega + kV_D$, $I_2(k, \theta) = \frac{\pi}{k R_{L,e}} \frac{\theta}{K_2(1/\theta)} e^{-\frac{\gamma_{\text{min}}}{\theta}} \left(1 + \frac{\gamma_{\text{min}}}{\theta}\right)$. The system is of second order as we assume $\Re(\tilde{\omega}) \ll kV_D$ (Evoli et al. 2018). In the test-particle limit (TPL) $\sigma I_2(k, \theta) \ll k^2 V_A^2$, $\sigma = k v_D \frac{n_{cr}}{n_i} \Omega_{ci}$, the linear growth rate $\Gamma \equiv \Im(\tilde{\omega})$ reads as

$$ \Gamma \approx \frac{\pi}{2} \Omega_{ci} \frac{1}{k R_{L,cr}} \frac{V_D}{V_A} \frac{n_{cr}}{n_i} \frac{\theta^2}{K_2(1/\theta)} e^{-\frac{\gamma_{\text{min}}}{\theta}} \left(1 + \frac{\gamma_{\text{min}}}{\theta}\right). \quad (13) $$

and the dispersion relation is almost unmodified $\Re(\tilde{\omega}) = kV_A$. In the relativistic regime $\gamma_{\text{min}}(k) \approx \theta/(kR_{L,cr})$ and $\theta^2/K_2(1/\theta) \to 0.5$ (for $\theta \gg 1$). We note that the dispersion in Eq. (13) can be derived using alternative approaches as it is the case in Bai et al. (2019); Plotnikov et al. (2021).

### 2.3. Higher frequency solutions: polarisation-dependent growth rates

As the simulations presented below depart from the TPL and have limited dynamical range in wavenumber, it is useful to relax the assumption $(\omega + kV_D) \ll \Omega_{ci}$. In this case, processes at frequencies $\Omega_{ci} < \omega < \Omega_{ce}$ are now included (e.g., ion-cyclotron resonances, excitation of Whistler waves). Yet, the two following approximations remain valid: $\Omega_{ci} \ll \Omega_{ce}$ and $\omega \ll \omega_{pe}$. The former is always satisfied for $m_i \gg m_e$ and the latter avoids looking for the resonances on the light-wave branch. In this case, the general contribution from the background plasma

$$ \chi_{bg} = \frac{c^2 \tilde{\omega}^2}{V_A^2 \omega^2} \frac{1 - \Omega_{ci}/\Omega_{ce}}{[1 \pm \tilde{\omega}/\Omega_{ci}][1 \pm \tilde{\omega}/\Omega_{ce}]}, \quad (14) $$

leads to a third-order dispersion relation

$$ k^2 V_A^2 = \omega^2 \frac{V_A^2}{c^2} + \frac{\tilde{\omega}^2}{[1 \pm \tilde{\omega}/\Omega_{ci}][1 \pm \tilde{\omega}/\Omega_{ce}]} + \frac{i\omega C}{kV_D}, \quad (15) $$

where $C \equiv \pi V_D \frac{\Omega_{ci}}{R_{L,cr}} \frac{n_{cr}}{n_p} \frac{\theta^2}{K_2(1/\theta)} I(k, \gamma_{\text{min}})$. After reorganising the terms

$$ \begin{aligned} & \tilde{\omega}^3 \frac{iC}{kV_D} + \tilde{\omega}^2 \left( \Omega_{ce} \Omega_{ci} - (kV_A)^2 - iC \pm \frac{iC \Omega_{ce}}{kV_D} \right) \\ & + \tilde{\omega} \Omega_{ce} \left( \mp(kV_A)^2 \mp iC + \frac{iC \Omega_{ci}}{kV_D} \right) \\ & - \Omega_{ce} \Omega_{ci} \left[ (kV_A)^2 + iC \right] = 0. \end{aligned} \quad (16) $$

Solutions of this equation can be calculated by considering the intersection of the real and imaginary planes solutions.

### 2.4. Magnetic field saturation

The resonant streaming instability has different saturation criteria. Holcomb & Spitkovsky (2019) discuss the saturation process in some details. The saturation is a two-step process. First, the linear growth rate at the fastest-growing wavenumber saturates. Then, in a second time, the drop of the growth over all wavenumbers induces a reduction of the mean drift speed of the CR distribution down to the local Alfvén speed. This second step properly corresponds to the instability saturation. In order to evaluate the total magnetic field strength at saturation we adopt the criterion proposed by Holcomb & Spitkovsky (2019); Bai et al. (2019) which balances the total momentum density lost by CRs $\Delta P_{cr} \approx (4/3)[(V_D - V_A)/c]n_{cr}\langle p \rangle$ with the momentum density acquired by forward propagating Alfvén waves $\Delta P_w \approx \rho V_A(\delta B^2/B_0^2)$. Here, $n_{cr}$ is the CR density and $\langle p \rangle$ the mean CR momentum, $\rho$ is the local gas (thermal proton) density. We find:

$$ \frac{\delta B_{\text{sat}}^2}{B_0^2} \approx \frac{8}{3} \frac{(V_D - V_A)}{V_A} \frac{n_{cr}}{n_i} \frac{m_e}{m_i} \frac{\langle \gamma v \rangle}{c}. \quad (17) $$

The above criterion likely provides an upper limit to the perturbed magnetic energy as it assumes that all momentum from
CR anisotropy is transferred to Alfvén waves. Some amount of momentum is transferred to the background gas (mean momentum and heating), forming the background fluid ‘wind’ (Holcomb & Spitkovsky 2019). Besides, numerical simulations always have a minimum of dissipation.

### 3. Particle-in-cell simulations

#### 3.1. Simulation method and setup

In this section, we use the PIC code SMILEI (Derouillat et al. 2018) to perform 1D3V simulations, where we retain one spatial direction, but all three components of velocities and electromagnetic fields (1D in space and 3D in momentum). The use of 1D geometry is justified by the expectation that the most dominant wavemodes driven by the instability are parallel to the external magnetic field (slab-type Alfvén modes with $k \parallel B_0$, see Sect. 2.1).

We set the simulation in the following way. The 1D box is filled with background electron-ion cold plasma at rest and a streaming beam of hot electron-positron plasma. The electron-positron pairs are inserted randomly in the simulation box. The number density ratio of the beam to background is set to be small, $n_{cr}/n_i \ll 1$. The external magnetic field is parallel to the simulation box $B = B_0e_x$. The typical setup and simulation parameters are listed as follows:

- The drifting beam (CRs) is composed of pairs and the background static plasma is composed of electrons and ions. Both components are charge and current neutral.

- We employ periodic boundary conditions for fields and particles.

- The box size is $L_x = 11264c/\omega_{pe} = 11264d_e$. This ensures that the grid length corresponds to at least a dozen of fastest-growing wavemodes.

- The electron skin depth, $d_e = c/\omega_{pe}$, is resolved with 6 grid points. The time step is $\Delta t = 0.95\Delta x/c$. (corresponding to CFL number equal to 0.95). We also tested twice smaller time step with no observed difference in simulation results.

- There are 100 macro-particles per cell per species, and we use 4-th order macro-particle shape factors<fnref n="3" />. This ensures a reasonably low level of Poissonian noise. As the physical density of the beam pairs is much smaller than the density of the background (i.e., $n_{cr} n_i \ll 1$) we use different macro-particle weights for different species. This effectively decreases numerical noise and increases the statistical significance of the distribution of beam particles (see Holcomb & Spitkovsky 2019).

- The typical mass ratio $m_i/m_e = 100$. But we also explored 1, 200, 300, 400, and 1000.

- The background plasma temperature $T_e = T_i = 10^{-3} m_e c^2$, that is $\theta_{bg} = 10^{-3}$ (in order to get the cell size not too far from the Debye length). The temperature is the same for ions and electrons.

- Number density ratio of CR (positrons) to background ions $n_{cr}/n_i$ between $10^{-3}$ and $10^{-2}$. As the beam and background plasma are both neutral, this implies that the ratio of CR electrons to background electrons ($n_{cr}/n_e$) is the same. The reference run – also called Fid run – corresponds to $n_{cr}/n_i = 10^{-2}$.

- The initial distribution function of the beam is drifting Maxwell–Jüttner distribution (see Eq. (9)) with temperature

<table><caption>**Table 2.** PIC simulation parameters presented in this work.</caption><thead><tr><th>Run ID</th><th>n<sub>cr</sub>/n<sub>i</sub></th><th>m<sub>i</sub>/m<sub>e</sub></th><th>Ω<sub>ce</sub>/ω<sub>pe</sub></th><th>V<sub>D</sub>/V<sub>A</sub></th><th>θ<sub>cr</sub></th><th>t<sub>sim</sub>ω<sub>pe</sub></th></tr></thead><tbody><tr><td>PIC-1<br>(Fid)</td><td>0.01</td><td>100</td><td>0.1</td><td>10</td><td>10</td><td>10<sup>6</sup></td></tr><tr><td>PIC-2</td><td>0.007</td><td>100</td><td>0.1</td><td>10</td><td>10</td><td>3 × 10<sup>5</sup></td></tr><tr><td>PIC-3</td><td>0.005</td><td>100</td><td>0.1</td><td>10</td><td>10</td><td>5 × 10<sup>5</sup></td></tr><tr><td>PIC-4</td><td>0.004</td><td>100</td><td>0.1</td><td>10</td><td>10</td><td>5 × 10<sup>5</sup></td></tr><tr><td>PIC-5</td><td>0.003</td><td>100</td><td>0.1</td><td>10</td><td>10</td><td>10<sup>6</sup></td></tr><tr><td>PIC-6</td><td>0.002</td><td>100</td><td>0.1</td><td>10</td><td>10</td><td>10<sup>6</sup></td></tr><tr><td>PIC-7</td><td>0.001</td><td>100</td><td>0.1</td><td>10</td><td>10</td><td>10<sup>6</sup></td></tr><tr><td>PIC-8</td><td>0.01</td><td>200</td><td>0.1</td><td>14.1</td><td>10</td><td>3 × 10<sup>5</sup></td></tr><tr><td>PIC-9</td><td>0.01</td><td>300</td><td>0.1</td><td>17.3</td><td>10</td><td>10<sup>6</sup></td></tr><tr><td>PIC-10</td><td>0.01</td><td>400</td><td>0.1</td><td>20</td><td>10</td><td>6 × 10<sup>5</sup></td></tr><tr><td>PIC-11</td><td>0.01</td><td>1000</td><td>0.1</td><td>31.6</td><td>10</td><td>10<sup>6</sup></td></tr><tr><td>PIC-12</td><td>0.005</td><td>50</td><td>0.1</td><td>7.1</td><td>10</td><td>3 × 10<sup>5</sup></td></tr><tr><td>PIC-13</td><td>0.005</td><td>200</td><td>0.1</td><td>14.1</td><td>10</td><td>10<sup>6</sup></td></tr><tr><td>PIC-14</td><td>0.005</td><td>300</td><td>0.1</td><td>14.1</td><td>10</td><td>10<sup>6</sup></td></tr><tr><td>PIC-15</td><td>0.005</td><td>1000</td><td>0.1</td><td>31.6</td><td>10</td><td>10<sup>6</sup></td></tr><tr><td>PIC-16</td><td>0.01</td><td>100</td><td>0.05</td><td>20</td><td>10</td><td>10<sup>6</sup></td></tr><tr><td>PIC-17<br>(e<sup>±</sup> bg)</td><td>0.001</td><td>1</td><td>0.01</td><td>14.1</td><td>1</td><td>10<sup>6</sup></td></tr></tbody></table>

**Notes.** The fiducial (or Fid) run is PIC-1.

of $T_{cr}$ (typically $\theta = 10$) and mean speed in x-direction is $V_D = 0.1c$. This sets the typical Larmor radius $R_{L,cr}$ of beam electrons (and positrons).

The simulation can be evolved up to $t\omega_{pe} = 10^6$, hat generally corresponds to $t\Omega_{ci} = 10^3$ (for the case of Fid run: $m_i/m_e = 100$ and $\omega_{pe}/\Omega_{ce} = 10$). In some cases, limited computational resources imposed the simulation to stop earlier. In the Table 2 we provide the typical parameters for different simulations presented in this study.

We note that in the Fid run, plasma beta parameter is $\beta \approx 8\pi n_i k_B T_i / B_0^2 = (v_{th,i}/V_A)^2 = 0.2$. This value remains unchanged for different mass ratios $m_i/m_e$ and different $n_{cr}/n_i$, explored in PIC simulations here. Hence, $\beta = 0.2$ is maintained for simulations PIC-1 to PIC-15.

#### 3.2. Simulation results

In this section, we describe the results of PIC simulations of the instability. We start with a description of the overall time-evolution of the Fid simulation, and then show how the linear growth rates compare to the analytical predictions. We then discuss the saturation of the instability and some properties of the self-generated wave spectrum.

##### 3.2.1. Overall evolution

The simulation starts with an anisotropic distribution of CRs. This anisotropy is small ($\delta f/f_0 \ll 1$) and is equivalent to an average drift velocity in x-direction that is small compared to c, $V_D \ll c$. If the initial drift is super-Alfvénic, $V_D > V_A$, then the instability is triggered. As a result, magnetic field grows exponentially during the so-called linear phase of the instability. At some point this growth ends (start of the post-linear phase) and,

<footnote n="3">³ In this case, the charge and current of one macro-particle is distributed over 5 spatial grid-points (Derouillat et al. 2018).</footnote>
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Overall evolution of the instability in time, for the Fid run. Top panel: wave energy evolution. Blue and red lines depict the evolution of transversal perturbations $\delta B_y^2/B_0^2$ and $\delta B_z^2/B_0^2$, respectively. The yellow line corresponds to the level of electrostatic fluctuations $E_x^2/B_0^2$. Bottom panel: mean velocity evolution of different species. Blue and red lines correspond to the background ions and electrons, respectively. The yellow and magenta lines show the evolution of the drift velocity of beam positrons and electrons, respectively.</figcaption></figure>

eventually, the level of magnetic fluctuations reaches a maximum (start of the saturated phase).

In Fig. 1, we show the evolution of the magnetic wave energy $\delta B^2/B_0^2$ (top panel) and the bulk velocities of different species of the beam-plasma system (bottom panel). This figure shows the result of the Fid simulation. The linear phase of the instability is clearly identifiable between $t\Omega_{ci} = 10$ and 100, where the wave energy growth is exponential (the scale of the y-axis is logarithmic). Immediately after the linear phase, the amplitude of magnetic fluctuations becomes large enough to scatter beam electrons and positrons efficiently. As a result, the streaming velocity of beam species decreases rapidly from $V_D = 10V_A$ to $V_D \approx 2V_A$ (see orange and magenta lines in the bottom panel). At some point (around $t\Omega_{ci} = 150$) the wave growth ends, and the instability enters the saturated phase. During the saturated phase (at $t\Omega_{ci} > 150$) the evolution becomes slow as the isotropization process of remaining particles takes more time. At the final time of the simulation ($t\Omega_{ci} = 10^3$) the level of magnetic fluctuations is $\delta B^2/B_0^2 = (\delta B_y^2 + \delta B_z^2)/B_0^2 \approx 0.06$, while the drift velocity of beam positrons decreased to $\sim V_A$ (i.e., isotropic in the wave-frame). The velocity of beam electrons is slightly below $V_A$. We impute this effect to the observation that some part of

self-generated waves in non-Alvénic (see detailed analysis of wave properties later in the article). We note that some momentum was also transferred to the background fluid, both electrons and ions reach $\langle V_i \rangle_{bg} \approx \langle V_e \rangle_{bg} \approx 0.07V_A$.

More detailed insight into the evolution of the instability can be gained by considering the wave spectrum and the distribution function anisotropy of beam species at different characteristic times. In Fig. 2, we present the wave spectrum $I(k)$ (top panels), spatial profile of magnetic field fluctuations $\delta B_{y,z}/B_0$ (middle panels), and the anisotropy of the distribution function of beam populations $f(V_x)$ (bottom panels). From left to right, we present three representative times: early linear phase ($t\Omega_{ci} = 7.6$), in the middle of the linear phase ($t\Omega_{ci} = 76$), and the late post-saturation time ($t\Omega_{ci} = 186$), respectively. The wave intensity in Fourier space for different polarisations $I_{\pm}(k)$ was done by transforming real-space magnetic fields $\delta B_{x,y}$ and according to $\int I(k)dk = \langle \delta B^2/B_0^2 \rangle$. Details of the method are presented in Appendix C.

During the initial time of the evolution, at $t\Omega_{ci} \sim 7.6$ (panels a-c) the wave spectrum is flat for both left-handed and right-handed wavemodes, the amplitude of magnetic fluctuations is negligible, and the distribution function of beam populations conserve its initial anisotropy. $f(V_x)$ conserves its initial linear slope, except only slight change near the pitch angle $\mu = 0$ (i.e., near $V_x = 0$). In the middle of the linear phase of the instability, at $t\Omega_{ci} = 76.0$, wave-modes with $0.1 \le kR_{L,0} \le 10$ are rapidly growing. The fastest growth occurs around $kR_{L,0} \approx 0.7$ as expected from the analytical calculation in Eq. (13). Panel e shows that the wave fluctuations reach $\delta B_{y,z}/B_0 \approx 0.1$ level, and the distribution function of beam populations starts to be noticeably affected (panel f). We note that the asymmetry between left-handed and right-handed modes (blue and red lines in panel d, respectively) is mirrored by a different behaviour of beam positrons and beam electrons (blue and red lines in panel f, respectively). During the late time when the instability saturated (panels g-i) the wave spectrum is nearly stationary pertaining some asymmetry between left- and right-handed modes. There is a lack of right-handed modes in a narrow band of wavenumber $5 \le kR_{L,0} \le 20$. The amplitude of spatial fluctuations reached $\delta B/B_0 \approx 0.3$, and the anisotropy in the distribution function of beam electrons and positrons is nearly erased. The latter is characterised by the quasi-flatness of $f(V_x)$.

An interesting feature is a line-like feature in the wave spectrum at short wavelengths at $kR_{L,0} \approx 170$ or $kd_e \approx 1.7$ (panels d and g). While being energetically subdominant to the bulk of excited wavemodes at $kR_{L,0} \le 10$ we, however, devote a subsection to the investigation of physical origin of this line feature at the end of this section (Sect. 3.2.5).

### 3.2.2. Linear phase

We now turn to more quantitative diagnostic of the linear phase of the instability. As expected, the accuracy of our numerical simulations (code and setup) can be tested versus the analytical expectations derived in the previous section. By performing spatial Fourier transform of the time-stacked magnetic field data FFT($B_{y,z}(x, t) = \tilde{B}_{y,z}(k, t)$ we can follow the time evolution of intensity in different wave-modes $k$. We then select the approximate time when the growth is exponential ($\Delta t_{lin}$) and then perform the linear fit of $|\log(\tilde{B}_{y,z}(k, \Delta t_{lin}))|$. This procedure provides the numerical growth rate of different wavemodes, available in the simulation. The smallest wavemode is $k_{min} =$
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Wave spectrum $I(k)$ of right-handed and left-handed polarisation waves (top panels: a, d, and g), spatial profile of magnetic field fluctuations $\delta B_{y,z}/B_0$ (middle panels: b, e, and h) and distribution function of beam electrons and positrons $f(V_x)$ and different times of the simulation (bottom panels: c, f, and i). Left column corresponds to early time ($t\Omega_{ci} = 7.6$), middle column to the linear phase ($t\Omega_{ci} = 76$) and right column to the saturated phase of the instability ($t\Omega_{ci} = 186$).</figcaption></figure>

$2\pi/L$ and the largest is $k_{max} = \pi/\Delta$, where $\Delta$ is the spatial grid cell size<fnref n="4" />.

The result for three different simulations is presented in Fig. 3. In this figure we show the growth rates derived from simulations (solid blue and red lines), the basic analytical growth rate in the test-particle limit (black dotted line) and the growth rates with no low-frequency approximation (orange and green lines). The blue-coloured line corresponds to the right-handed modes ($I_+$) while the red-coloured lines correspond to the left-handed modes ($I_-$). Different panels present results from simulations with different beam-to-background number density ratio: $n_{cr}/n_i = 0.01, 0.005, 0.002$ from top to bottom panel, respectively.

We can note that the maximum linear growth rate is obtained at $kR_{L,cr} \approx 0.7$ in agreement with analytical expectations from Equation (13), for all of three simulations. The growth rates are reasonably well reproduced up to $kR_{L,cr} \sim 5-6$. At higher $k$ some differences appear in the different mode polarisation. Circular right-handed modes always produce higher growth rates than circular left-handed ones. This can be explained by adopting the high-frequency relation derived from Eq. (16). However, it is difficult to provide a firm conclusion because of the increasing effect of numerical noise with $k$. The agreement with analytical expectations is better for lower $n_{cr}/n_i$, as can be seen by comparing the top panel with the bottom panel. This is expected since analytical calculations have been derived in the quasi-linear limit.

### 3.2.3. Saturated phase

As previously discussed in Sect. 3.2.1 it can be identified that after the linear phase, where the wave growth is exponential, the wave-saturation occurs rapidly (see top panel in Fig. 1). It corresponds to the phase when the wave intensity reaches its maximum and remains nearly constant afterwards. This level defines the saturation of wave intensity as $\delta B_{sat}^2 = \delta B_{y,end}^2 + \delta B_{z,end}^2$.

In Fig. 4, we show the dependence of saturation level of the magnetic field, $\delta B_{sat}^2/B_0^2$, on the beam density $n_{cr}/n_i$ (panel a) and on the mass ratio $m_i/m_e$ (panel b). The blue circle symbols present the results from different simulations where $n_{cr}/n_i$ varies from 0.001 to 0.01; all other parameters being kept the same ($m_i/m_e = 100, V_D = 0.1c, V_A = 0.01c$. See Table 2). The dashed black line follows the level derived from the criterion assuming that all of the available momentum lost by the decelerated CR beam (from $V_D$ to $V_A$) is transferred to the waves. This criterion was defined in Eq. (17)<fnref n="5" />. The values measured in the simulations are slightly below the expectation but show very good agreement with the expected trend. We note that such a good agreement is remarkable for a PIC simulation and owing to the simplicity of the analytic criterion<fnref n="6" />.

<footnote n="4">⁴ Waves at $k > \pi/\Delta$ are not resolved, as the Nyquist wavenumber is $2\pi/\Delta$.</footnote>

<footnote n="5">⁵ We used the fact that for an ultra-relativistic Maxwell–Jüttner distribution (see Eq. (9)) with $\theta = k_B T_{cr}/m_e c^2 \gg 1$ the average particle momentum is $\langle p \rangle = 3k_B T_{cr}/c$</footnote>

<footnote n="6">⁶ The exception is the simulation with the lowest $n_{cr}/n_i = 10^{-3}$. Here, the simulation was too short to reach saturation. We directly verified that the linear phase was still ongoing at the final time of the run.</footnote>
<figure><img src="image_3.png" /></figure>

**Fig. 3.** Linear growth rate of the instability $\Gamma$ as function of the wavenumber $k$ for two polarisation modes: circular right-handed (blue solid line) and left-handed (red solid line). Results for the simulations with three different values of CR density are presented: top panel corresponds to $n_{cr}/n_i = 0.01$, middle panel corresponds to $n_{cr}/n_i = 5 \times 10^{-3}$, and bottom panel corresponds to $n_{cr}/n_i = 2 \times 10^{-3}$. In each panel, the solid orange and green lines correspond to the theoretical growth rate using the solution of the full dispersion relation Eq. (16). Black dotted line follows the expectation in test-particle limit from Eq. (13).

We also derived the saturation in simulations with different $m_i/m_e$ from 50 to 1000 (panel b of Fig. 4). The same criterion provides a good approximation of the saturation level for both $n_{cr}/n_i = 5 \times 10^{-3}$ and $n_{cr}/n_i = 10^{-2}$ cases. The ‘lack’ of waves is small and can probably be imputed to other channels of dissipation (e.g., background plasma heating, some amount of wave damping).

It is interesting to note that the acquired final velocity of background $e^-$-ion plasma is roughly equal to the saturated wave intensity, $\langle V_i \rangle_{bg}/V_A \approx \delta B_{sat}^2/B_0^2$. This was found in all simulations where both quantities could be measured with good confidence. The wave saturation criterion assumes that the momentum lost by CRs when decelerating from $V_D$ drift to $V_A$

<figure><img src="image_4.png" /></figure>

**Fig. 4.** Saturation level of the magnetic fluctuations. Panel a: wave saturation intensity as function of $n_{cr}/n_i$. Blue circles present the results from 1D3V PIC simulations. The red dashed line follows the scaling from Eq. (17) that balances the momentum loss by CRs with momentum gained by waves. The data point at $n_{cr}/n_i = 10^{-3}$ is a lower limit as the simulation ended before a saturated state was reached. Panel b: same but varying $m_i/m_e$ at fixed $n_{cr}/n_i = 10^{-2}$ (blue symbols) and $n_{cr}/n_i = 5 \times 10^{-3}$ (red symbols).

is transferred to waves. Yet, the Alfvén waves are carried by the background fluid that seems to be dragged to the mentioned equality in the process of saturation of the instability.

### 3.2.4. Properties of the wave spectrum

Overall intensity is only one aspect of the self-generated waves. The properties of wave spectrum may require more subtle analysis. The wave spectrum after saturation in the Fid simulation, $I_{sat}(k) = \delta B_{sat}(k)^2/B_0^2$, is shown in Fig. 5, top panel. The blue and red solid lines correspond to the final spectrum of right- and left-handed circularly polarised waves, respectively. Previously presented asymmetry between two polarisation modes is observed at $kR_{L,cr} > 6$. A line-like feature is clearly seen at $kR_{L,cr} \sim 170$.

In order to understand different cut-off positions for two modes we consider in more detail the proper wavemodes and theoretical growth rates. The corresponding result is presented in Fig. 5, bottom panel. The solutions of the full dispersion relation are obtained by solving Eq. (16). Both real and imaginary (unstable) parts are shown using solid and dashed lines, respectively. The light-wave is shown with a black solid line. Alfvén wave corresponds to low frequency and low-$k$ part of red and blue solid lines. At $kR_{L,cr} \ge 4-6$ the two branches separate into Whistler (blue solid line) and Ion-cyclotron (red solid line). Alternatively, these two wave branches can be analytically derived with more ease by setting $C=0$ (CR term) in the dispersion relation. Solving the equation without the CR term then provides accurate solutions for Whistler and Alfvén/ion-cyclotron (AIC) dispersion relations (e.g., Vainio 2000); which are identical to $\mathcal{R}$-labelled solutions shown in Fig. 5 (bottom panel). By comparing the dashed red/blue lines and solid red/blue lines, one can clearly see that the modes at $kR_{L,cr} > 10$ are growing into the Whistler (right-handed) and ion-cyclotron (left-handed) wavemodes. This
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Wave spectrum at saturation and corresponding plasma wave-modes. Top panel: wave spectra in Fid simulation. The blue line for right-handed waves and red line for left-handed waves. Bottom panel: analytic dispersion relation of waves obtained by solving Eq. (16). The real part of $\omega$ (solid lines) and instability growth rate for the two polarizations that is the imaginary part of $\omega$ (dashed lines). The light-wave is shown with black solid line. The red star symbol at the top left side locates the intersection with $\omega = 2\omega_{pe}$ line. Alfvén wave corresponds to low frequency and low-$k$ part of red and blue solid lines. At $kR_{L,cr} > 1$ the two branches separate into Whistler (blue solid line) and ion-cyclotron (red solid line). Vertical arrows in bottom panel locate the transition into ion-cyclotron regime (at $kd_i = 1$) for left-handed waves, and electron-cyclotron regime $kd_e = 1$ for right-handed waves. The lower red star symbol locates the intersection between the ion-cyclotron branch (red solid line) and $\Omega_{ce} - 3k\nu_{th,i}$ (red dotted line).</figcaption></figure>

is the first hint to understand the asymmetry between $I_+$ and $I_-$ spectra.

We now test the idea that the left-handed polarisation waves ($I_-$) are affected by ion-cyclotron resonance with background plasma. In Fig. 6, we present the saturated spectra of right-handed waves ($I_+$; top panel) and of left-handed waves ($I_-$; bottom panel). Lines of different colours correspond to different mass ratio adopted in the simulations: blue line corresponds to $m_i/m_e = 1000$ (Fid run), red line corresponds to 400, and orange line corresponds to 100. As can be seen in the top panel, $I_+$ spectrum is only weakly affected by the variation of $m_i/m_e$, while the $I_-$ branch (bottom panel) exhibits a cut that transits to lower

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Saturated spectra of right-handed waves ($I_+$; top panel) and of left-handed waves ($I_-$; bottom panel). Lines of different colour correspond to different mass ratio adopted in the simulations: blue line corresponds to $m_i/m_e = 1000$ (Fid run), red line corresponds to 400, and orange line corresponds to 100. Inset: cut-off wavenumber of left-handed branch for different values of $m_i/m_e$. Blue circles correspond to values derived from simulations with $n_{cr}/n_i = 10^{-2}$ and red circles are from simulations with $n_{cr}/n_i = 5 \times 10^{-3}$. The dashed black line follows $k_{cut}d_i = 0.5$.</figcaption></figure>

$kR_{L,cr}$ with increasing mass ratio. The position of the cut is highlighted using vertical dashed lines. For $m_i/m_e = 100$ the cut is located at $kR_{L,cr} \approx 6$, while for $m_i/m_e = 1000$ it is located at $kR_{L,cr} \approx 1.7$.

More systematic analysis of $k_{cut}$ value in large set of simulations is presented in the inset of the bottom panel. Blue circles correspond to values derived from simulations with $n_{cr}/n_i = 10^{-2}$ and red circles are from simulations with $n_{cr}/n_i = 5 \times 10^{-3}$. All other parameters are the same as in Fid run. The black dashed line follows $k_{cut}d_i = 0.5$ scaling. We observe that all measured values cluster around this scaling. The basic interpretation is that all of $I_-$ waves are suppressed for any $kR_{L,cr} > 0.5$ that falls in the regime where the Alfvén waves transition into Ion-cyclotron branch.

An accurate location of resonant wavenumber can be found by considering the intersection between AIC branch (red solid line in bottom panel of Fig. 5) and $\omega = \Omega_{ci} - 3k\nu_{th,i}$ (red dotted line in the same panel)<fnref n="7" />. This gives the cut-off at $k_{cut}d_i = 0.47$,

<footnote n="7">⁷ We found that the factor $3\nu_{th,i}$ in this relation provides better agreement with simulation, instead of more standard $\nu_{th,i}$. It was argued in Engelbrecht & Strauss (2018) that this factor captures better the range of resonantly interacting background particles.</footnote>
that is basically the same as the result plotted as the dashed line in the inset of Fig. 6. In general, we found that the cut-off wavenumber agrees well with the scaling (Engelbrecht & Strauss 2018; Strauss & le Roux 2019):

$$k_{\text{cut},-} = \frac{\Omega_{\text{ci}}}{V_{\text{A}} + 3\nu_{\text{th},i}}. \qquad (18)$$

Left-handed waves with $k > k_{\text{cut},-}$ are expected to be strongly damped by background plasma. An improved calculation of the damping term following Strauss & le Roux (2019) (Appendix B) shows that the resonance with first two ion-cyclotron harmonics suppresses all left-handed waves at $k > k_{\text{cut},-}$. Numerical solutions of the dispersion relation accounting for thermal effects in background plasma – using plasma dispersion function – shows as well that $k_{\text{cut},-}$ is an accurate estimate of both the measured values and of full solution of the dispersion relation. To convince ourselves that this is not a numerical effect, we note that the damping of ion-cyclotron waves is a well-observed process in solar wind plasma (Woodham et al. 2018; Bowen et al. 2022, for recent studies). The experimental cut-off scale agrees remarkably well with what we find in the present simulations.

Concerning the $I_+$ spectrum, the value of wavenumber ($k_{\text{cut},+}$) at which it drops to noise floor level should be regulated by the interaction with the Whistler branch. Using the physical parameters in Fid simulation, we found that the intersection of Whistler branch with $\omega = \Omega_{\text{ce}} - 3\nu_{\text{th},e}$ occurs at $kd_e \approx 0.56$. This value is in good agreement with Strauss & le Roux (2019) (see their Fig. 1 and Eqs. (61)–(64)). However, Fid simulation shows that the cut-off of right-handed modes occurs at lower $kd_e \approx 0.1$ (i.e., $kR_{L,cr} \approx 10$), see top panel of Fig. 5. This value is sensitively lower than expected $k_{\text{cut},+}$. By inspecting in more detail, we found that there are indeed growing waves up to $kd_e \sim 0.3-0.5$ earlier in the simulation. However, the part of spectrum at $kd_e > 0.1$ is damped during the saturated phase. Alternatively, in the simulations with lower CR number density ($n_{cr}/n_i \le 5 \times 10^{-3}$) the waves at $0.1 < kd_e < 0.5$ emerge more from the background level and remain at relatively low-intensity level. We have no satisfactory explanation for this different damping effect that increases for larger $n_{cr}/n_i$. Yet, no wave growth is observed for $kd_e > 0.5$ in any simulation, that is consistent with expected $k_{\text{cut},+}$ value.

As a side check, we performed a simulation where the background plasma is composed of pairs ($m_i = m_e$, run PIC-17 in Table 2). The detailed presentation is reported in Appendix A. To summarise, we found that there is no asymmetry between two polarisation modes of self-excited waves. The initial drift of CRs (beam population) reduced from initial $V_D = 14.1 V_A$ to exactly $V_A$ for both CR electrons and positrons in an identical way. This confirms the idea that the asymmetry is due to the ionic nature of background plasma.

### 3.2.5. Line feature in the wave spectrum

A minor effect of the interaction between CR beam and background plasma produces a high-$k$ line-like feature (e.g., feature at $kR_{L,cr} \approx 170$ in Fig. 2). It was observed in the wave spectra of all simulations except when background plasma is composed of pairs (see Appendix A). As the spectra are only in $k$ space we possess limited information on the frequency of this wave. In order to clarify the nature of this line, we performed dispersion analysis to derive $\omega(k)$ from the simulation. This was done by storing the EM fields on the spatial grid at every time-step between $t\Omega_{\text{ci}} = 95$ ant $t\Omega_{\text{ci}} = 96.6$. It corresponds to the

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Wave dispersion analysis. Top panel a: wave dispersion $\omega(k)/\omega_{\text{pe}}$ at the final times of the linear phase. From figure simulation (see Table 2 for parameters). The diagnostic was performed on the $\delta B_z$ component between $t\Omega_{\text{ci}} = 95$ ant $t\Omega_{\text{ci}} = 96.6$ (at the end of the linear phase of the instability). Electromagnetic light-wave branch with left and right-handed polarisation is easily identifiable. At low $\omega < \Omega_{\text{ce}}$ the Whistler branch is clearly seen as well (theoretical dispersion line over plotted with dashed line). Plasma emission at $\omega = \omega_{\text{pe}}$ is observed and the second harmonic can be observed as well, especially in the brighter zone where it intersects the EM wave. Bottom panel b: $I_+(k)$ spectrum at saturation from Fid simulation. It is the same spectrum as in Fig. 5 but plotted using linear x-axis scale and wavenumber plotted in units of $d_e^{-1}$. This improves the visibility of the high-$k$ part of the spectrum.</figcaption></figure>

end of the linear phase of the instability. By stacking the snapshots of the stored fields and performing 2D FFT on a given field component we derive the dispersion relation $\omega(k)$.

In Fig. 7, we present the $\omega(k)$ dispersion spectrum performed on $B_z$ component of the magnetic field. The light-wave (EM wave) is well-seen with dispersion $\omega^2 \approx k^2c^2 + \omega_{\text{pe}}^2$. As the plasma is magnetised, at low-$k$ the EM wave separates into right- and left-handed circular polarisation branches that have slightly different cut-off frequencies. The right-handed mode has the cut-off frequency at $\omega_{\text{RH}} \approx \omega_{\text{pe}} + \Omega_{\text{ce}}/2$ while the left-handed branch has the cut-off frequency at $\omega_{\text{LH}} \approx \omega_{\text{pe}} - \Omega_{\text{ce}}/2$ (Kulsrud 2005). Exact – textbook – solutions for both branches are over-plotted using black dashed lines. It can be seen that there is a region with enhanced intensity on this branch. This corresponds to a resonance between the EM wave and second harmonic $\omega = 2\omega_{\text{pe}}$. In terms of the wavenumber, theoretically, it corresponds to $kd_e = \sqrt{3}$. That is in good agreement with the measured position of the instantaneous spectrum $I(k)$, shown in the (b) panel of the figure. This spectrum is the same as presented in Fig. 2 and Fig. 5 but using linear x-axis scale and the wavenumber normalized to $d_e^{-1}$ instead of previously used $R_{L,cr}^{-1}$. In this panel,
one can identify the dominant contribution of $kd_e < 0.1$ waves and a much less intense region slightly below $kd_e = \sqrt{3}$. Both panels combined lead to a firm conclusion that this line feature is the second-harmonic emission. We verified in a large set of simulations that the position of this line is always at $kd_e \approx \sqrt{3}$, independently of $n/n_i$, $m_i/m_e$, and $T_{cr}$ values explored.

There is extensive literature on PIC simulations of second-harmonic plasma emission in unmagnetized plasma (e.g., Kasaba et al. 2001; Henri et al. 2019; Krafft & Savoini 2021) and magnetized plasma (e.g., Lee et al. 2022). In the present work, we do not focus on this interaction and limit on basic observation. This is also justified by the fact that this line-emission is largely subdominant to the main instability generated waves in Alfvén branch. It does not seem to have an important effect on the dynamics of the pair beam-driven instability and, more generally, on particle scattering properties.

For completeness, let us describe other noticeable features in Fig. 7. At low $\omega < \Omega_{ce}$ the Whistler branch is clearly seen as well (theoretical dispersion line over plotted with dashed line). The dynamically important part (and by far the brightest spot) of the spectrum corresponds to the region $\omega/\omega_{pe} < 0.1$ and $kd_e < 0.1$. It is not resolved in this figure as the time window for diagnostics was too small.

## 4. Discussion

This article addressed the pair-driven instability when the beam propagates through ionic background plasma at super Alfvénic speed. We used fully kinetic particle-in-cell simulations in order to understand the non-linear evolution of the interaction. Combined with detailed analytical derivation in appropriate approximation we verified the accuracy of linear phase of the instability and confirmed some intuitive features, such as a decrease of the beam velocity to near-Alfvénic speed. The intensity of waves at saturation was found to be in good agreement with the momentum transfer criterion: momentum lost by beam particles when they decelerate (erase the anisotropy in the distribution function) is transferred to Alfvén waves. Both Alfvén and kinetic scales were covered by the wavemodes triggered by the instability. This leads to important damping of waves belonging to the ion-cyclotron and Whistler branches. All short-scale waves with $kd_i \gg 1$ are strongly suppressed by ion-cyclotron resonances with background plasma. We also observed a second harmonic plasma emission. These results might be of interest as input for detailed modelling of particle propagation in TeV halos.

Let us note, however, that by construction our analysis is limited to parallel waves as our simulations were 1D in space. There is no contribution from $k_\perp$ related phenomena. Yet, it is known that at kinetic scales with $k_\perp d_i > 1$ there is a dominant contribution of kinetic Alfvén waves (KAW), that is well-studied in multidimensional kinetic simulations (e.g., Cerri et al. 2016, 2017). A study of the interplay between self-generated parallel waves by the streaming instability and external turbulence is a topic for future work (recently, a theoretical framework to address this problem was proposed by Cerri 2024).

### 4.1. Simulations rescaling in comparison with realistic halos

Here, we discuss the relevance of our work to TeV halos, limitations, and possible extensions using other numerical techniques. To start, let us examine how realistic the parameter set chosen in our simulations is. We consider the case of the Geminga pulsar to provide a rough estimate of the ratio $n_{cr}/n_i$ in a realistic halo.

The calculation is based on the parameters adopted in Schroer et al. (2023), namely a typical pulsar age of 342 kyrs, a present spin-down luminosity of $L_{age} = 3.26 \times 10^{34}$ erg/s and an initial spin-down timescale $t_0 \approx 12$ kyrs. Under the assumption of pulsar spin down from magnetic dipole radiation (breaking index of 3), the luminosity at time $t$ is

$$L(t) = L_{age} \frac{(1 + t_{age}/t_0)^2}{(1 + t/t_0)^2}. \qquad (19)$$

This implies an initial spin-down luminosity $L_0 = 2.9 \times 10^{37}$ erg/s at early times $t \ll t_0$, three orders of magnitude above the present-day value. The particle injection spectrum follows a broken power law<fnref n="8" />

$$ \begin{aligned} Q(E, t) &= Q_0(t) \left(\frac{E}{E_b}\right)^{-a_L} & \text{for } E \le E_b \\ &= Q_0(t) \left(\frac{E}{E_b}\right)^{-a_H} & \text{for } E > E_b, \end{aligned} $$

where we use $E_b \sim 200$ GeV, $a_L = 1.6$, $a_H = 2.4$, in agreement with the ranges of values considered by Schroer et al. (2023). The particle population is energized by the conversion of a constant fraction $\eta \approx 1$ of the pulsar spin-down power, such that the normalisation $Q_0$ of the injection spectrum is obtained by

$$ \int_{E_{\min}}^{\infty} dE E Q(E, t) = \eta L(t), \quad \eta \le 1. \qquad (20) $$

The minimum particle energy is taken as $E_{min} = 1$ GeV. The average particle energy in the injection spectrum can be approximated as

$$ \begin{aligned} \langle E \rangle &= H(a_L, a_H, r) \times E_b \\ &\approx \frac{(2 - a_L)^{-1} + (a_H - 2)^{-1} - (2 - a_L)^{-1} r^{2-a_L}}{(a_H - 1)^{-1} - (a_L - 1)^{-1} + (a_L - 1)^{-1} r^{1-a_L}} \times E_b, \end{aligned} \qquad (21) $$

where $r = E_{min}/E_b$. For the adopted parameters, this yields $H \approx 0.12$ hence $\langle E \rangle \approx 24$ GeV. If $E_{min} = 100$ MeV then $H \approx 0.03$ and $H$ goes as $r^{0.6}$ as $r \to 0$. The factor $H$ is weakly sensitive to $a_L$ and $a_H$, for instance keeping $E_{min} = 1$ GeV with $a_L = 1.8$ (1.4) and $a_H = 2.2$ (2.6) we find $H \approx 0.1$ (0.17).

We assume that, at each time $t$, the particle injection rate $Q(E, t)$ is directly released across the surface of the nebula $S_{PWN}$. The present-day nebula can be approximated as a sphere with radius $R_{pwn} = 0.1$ pc (Posselt et al. 2017), but it was likely much larger in the past, when the pulsar was more powerful. The typical size of a bow-shock nebula is the so-called stand-off distance, where pressure in the spherically expanding pulsar wind equals the ram pressure exerted by the ambient medium on the moving pulsar (see Eq. (16) in Gaensler & Slane 2006). This stand-off distance scales as $\sqrt{L}$ and so $R_{pwn}(t) = R_{pwn}(t_{age}) \sqrt{L(t)/L_{age}}$ and $S_{PWN} \propto L$.

The number density of pairs injected by the nebula at time $t$ can therefore be written as

$$ n_{cr}(t) = \frac{\eta L(t)}{\langle E \rangle S_{PWN} V_D}. \qquad (22) $$

<footnote n="8">⁸ The injection spectrum is also expected to feature a time-varying cut-off at an energy $E_c$, for instance to reflect the progressive decrease of the maximum potential drop, which would lead to $E_c(t) = 1.7 \sqrt{L(t)/10^{36}} \text{ erg s}^{-1} \text{ PeV}$. We neglect this here because the high typical values for this cutoff energy have no impact on our particle density estimates.</footnote>
For a meaningful comparison with the situation implemented in the simulations, we assume that the particles stream across the surface at a velocity $V_D = 10V_A$, where $V_A$ is the Alfvén speed in the surrounding medium. This numerical choice seems reasonable since $V_D$ is basically a free parameter. The typical background gas density around Geminga can be estimated from the absence of $H_\alpha$ emission surrounding the Geminga bow shock to be in the range between $n_i \sim 10^{-3} \text{ cm}^{-3}$ (Caraveo et al. 2003) and $n_i \sim 10^{-2} \text{ cm}^{-3}$ (Posselt et al. 2017). We use the latter higher bound together with a magnetic field strength of $1\mu\text{G}$ to compute $V_A = 22 \text{ km s}^{-1}$.

According to Eq. (22), the particle density immediately beyond the surface of the nebula does not depend on time because $S_{\text{PWN}} \propto L$ and all other parameters are considered constant. Using present day values for $L = 3.26 \times 10^{34} \text{ erg s}^{-1}$, $R_{\text{pwn}} = 0.1 \text{ pc}$ and other parameters being specified above (i.e., $\eta = 1$, $V_D = 10V_A = 220 \text{ km s}^{-1}$, $\langle E \rangle = 24 \text{ GeV}$), we get $n_{\text{cr}} \sim 3 \times 10^{-8} \text{ cm}^{-3}$. Such a value may be considered as a lower limit because the above calculation assumes instantaneous transfer of the pulsar power to the escaping particle flux, whereas there is most likely an accumulation of particles within the nebula for centuries or millennia or more before they can escape or lose their energy. On the other hand, the values of $n_i$ and $B_0$ adopted above favour a small Alfvén speed in the surrounding medium, hence a small escape velocity $V_D$ and therefore a higher particle density. Overall, the beam-to-background density ratio $n_{\text{cr}}/n_i$ that can be expected close to the injection zone, that is close to the pulsar wind nebula, is likely in the range $10^{-6} - 10^{-5}$ (with the caveat that the background density at early times might have been higher than $10^{-2} \text{ cm}^{-3}$). The PIC simulations presented here are therefore not representative of the conditions in the Geminga system.

Another important difference with the astrophysical case resides in our choice of particle distribution. In our simulations, we consider a Maxwell–Jüttner distribution and most of the energy is carried by particles around 10 MeV instead of 200 GeV. Hence, a realistic simulation should have a beam-to-gas density ratio in the range $10^{-6} - 10^{-5}$ and a beam with a bulk energy of $E_b \sim 10^5 m_e c^2$. It seems that if the resonant streaming instability can be triggered in the environment of evolved pulsars it should be in a region close to the pulsar wind nebula, otherwise the surface term in Eq. (22) produces too much beam dilution. Actually, these values require an unrealistically large numerical box as the PIC method has to resolve the electron skin depth (ideally Debye length). Such high energies would require the use of larger-scale simulations, including a magnetohydrodynamic (MHD) description of the plasma. We therefore discuss briefly in the next subsection an alternative approach based on combining PIC and MHD but leave a more complete study of this regime for future work.

We note that we did not provide a direct extrapolation of the results from our PIC experiments to astrophysical scales for the following reasons: (i) the convergence of simulations with low values of $n_{\text{cr}}/n_i \le 10^{-3}$ was not reached, (ii) the saturation level of waves confirms intuitive scaling that is already used in the literature in general, and (iii) employing specific non-periodic boundary conditions is needed for direct comparison with a given astrophysical system.

## 4.2. Large-scale simulations: MHD-PIC

As discussed above, the main limitation is the scale separation that we are able to achieve with full-PIC simulation. This was reflected by the adopted distribution of CRs with an average

<figure><img src="image_8.png" /></figure>

**Fig. 8.** Comparison of the results from full-PIC and MHD-PIC simulations using identical parameters as in run *PIC-3* given in Table 2. The top panel presents the average magnetic energy evolution in time for full-PIC (blue solid line) and MHD-PIC (red solid line) techniques. The inset shows the full time-evolution from MHD-PIC simulation, as it runs over a much longer timescale than any full-PIC run. The bottom panel shows the comparison of wave spectra at the time that corresponds to the final time of the full-PIC simulation. Blue and red solid lines show the spectra from full-PIC for right- and left-handed modes, respectively. Green and orange solid lines show the same for MHD-PIC simulation. Note: the spectra from MHD-PIC simulation were upshifted by a factor of ≈ 2 in order to get a better visual comparison.

energy of $E \sim 10$ MeV. Using MHD-PIC method (Bai et al. 2015; van Marle et al. 2018; Mignone et al. 2018), we can reach higher energies and lower beam densities, which reflect more closely the type of beam-plasma interaction found in astrophysical circumstances. To demonstrate the possibilities of the MHD-PIC method, we performed two simulations in a 2-D box measuring $50 \times 0.5 R_L$, with $R_L$ the Larmor radius of the peak of the beam energy distribution. This box is subdivided into a grid of 2000 by 20 grid cells. All boundaries are set to periodic. The box is filled with a stationary thermal plasma and a magnetic field aligned with the x-axis in which the Alfvén speed is set at $0.01 c$, the plasma being in equipartition. Through this plasma, we run a beam of non-thermal particles with an advection speed of $0.1 c$, and a thermal particle distribution according to a Maxwell–Jüttner prescription with the partial motion being isotropic in the co/moving frame. This beam is made up of an average of 50 particles per particle species per grid cell. We run two MHD-PIC simulations. One with a temperature of $T_{\text{cr}} = 10 m_e c^2$ and a relative particle density of $n_{\text{cr}}/n_i = 5 \times 10^{-3}$ for comparison with the PIC simulations. The second simulation has a beam temperature of $T_{\text{cr}} = 100 m_e c^2$ and a relative particle density of $n_{\text{cr}}/n_i = 10^{-3}$ to demonstrate MHD-PIC’s ability to reach
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Magnetic energy evolution of a MHD-PIC simulation with $T_b = 100m_e c^2$ and relative beam density $n_{cr}/n_i = 10^{-3}$.</figcaption></figure>

a more realistic astrophysical parameter space. The simulation is allowed to run while the gas and particles relax.

For the simulation with $T_{cr} = 10 m_e c^2$ and $n_{cr}/n_i = 5 \times 10^{-3}$ (same parameters as full-PIC simulation PIC-3). The result of this comparison is presented in Fig. 8. As can be seen in the upper panel of the figure (magnetic energy evolution), during the linear phase the slope in both PIC and MHD-PIC is very close. Meaning that the fastest growth rates are in good agreement with both techniques. We found that the peak magnetic wave intensity value was lower in MHD-PIC by a factor of ≈ 2 (difference in peak values for blue and red lines). Yet, the wave spectra (bottom panel) were comparing reasonably well, except for the absence of the cut in wave spectra of left-handed modes at $kd_i > 0.6$. This is expected, as the resonant ion-cyclotron absorption effect is absent if the background plasma is treated using an MHD formalism.

For the simulations with a temperature of $100 m_e c^2$, the resulting magnetic wave energy evolution in our simulation is shown in Fig. 9. The maximum value $\delta B_{sat}^2/B_0^2 \sim 0.05$ is reached at $t\Omega_{ci} \approx 400$. The level of magnetic fluctuations at the saturation state is consistent with the expected value from Eq. (17) which gives $\delta B_{sat}^2/B_0^2 \sim 0.05$ for the adopted parameters (we also used $m_e/m_i = 100$). From Fig. 9, we can also estimate the growth rate of the magnetic field instability for the linear phase, which lasts from $t\Omega_{ci} \approx 90$ until $t\Omega_{ci} \approx 180$. From our simulation, we find a growth rate over that period of approx. $\Gamma/\Omega_{ci} \sim 3 \times 10^{-3}$ in reasonable agreement with the expected value from Eq. (13).

This technique is much less computationally prohibitive than full-PIC simulations and could potentially lead to regimes which are closer to astrophysical parameters. Some improvements such as $\delta f$ method could be implemented to reach better signal-to-noise ratio and distribution function representation (e.g., Denton & Kotschenreuther 1995; Bai et al. 2019) as well as non-periodic boundary conditions (Bai 2022) that are important to properly describe the saturation of the instability. This work is deferred to a forthcoming study.

## 5. Conclusion

To conclude, we investigated in detail the evolution of a resonant electron–positron beam-driven streaming instability that develops in a cold background electron-ion plasma. We considered a pair beam following a Maxwell–Jüttner distribution and drifting with respect to the background plasma at a super-Alfvénic speed with initial value $V_D = 10V_A$. The simulations were performed using 1D3V Particle-in-Cell simulations.

We find the following main results:

1. Maximum linear growth rates at scales $kR_{L,cr} \sim 0.7$ are well recovered. The solutions at smaller scales show polarisation-dependent growth rates, with forward right modes taking over left-handed modes. This effect can be explained by relaxing the condition over the wave frequency in the moving frame to be $\ll \Omega_{ci}$. At high $k$ the right and left-handed modes follow the Whistler/electron-cyclotron and the Alfén/ion-cyclotron branches respectively. The simulations show a reasonable agreement with the analytical solutions. At $k \gg d_i^{-1}$ left-handed modes are damped due to ion-cyclotron resonance and the right-handed modes are cut due to the interaction with the Whistler branch;

2. We have checked that the criterion for magnetic saturation by momentum transfer from the beam to the waves is well reproduced, with the exception of the low beam density case where the simulations have not converged.

Our results underline the importance of kinetic simulations when subtle effects determine the saturation level and properties of the self-generated spectrum.

We also compared PIC simulations with 1D3V MHD-PIC simulations using the same setup. In particular, we find a reasonable agreement between the two methods. MHD-PIC results show a bit more dissipation, leading to a typical saturated magnetic field energy that is a factor ≈2 below the one obtained using PIC simulations (see Fig. 8).

We finally discuss the rescaling of our numerical setup to the context of gamma-ray halos around evolved pulsars. We find that the kinetic numerical setups have beam-to-background density ratios still several orders of magnitude above the one currently expected in the Geminga pulsar halo. Moreover, the non-thermal particle population carrying the bulk of momentum has typical kinetic energies several orders of magnitude above the range covered by PIC simulations and is unreachable for now. MHD-PIC simulations may partly fill the energy gap to describe the gamma-ray halo more properly. This deserves a future investigation.

*Acknowledgements.* The authors are grateful to C. Evoli and B. Olmi for fruitful discussions. C. Evoli also provided us with an estimation of the beam-to-background plasma density ratio, which was helpful for the discussion in Sect. 4.1. I.P. is grateful to Silvio Sergio Cerri and Gabriel Fruit for insightful suggestions that improved our understanding of this topic. We acknowledge financial support by ANR for support to the GAMALO project under reference ANR-19-CE31-0014. The simulations have been performed on CALMIP supercomputing resources at Université de Toulouse III (France) under the allocation P20028.

## References

Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017a, Science, 358, 911

Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017b, ApJ, 843, 40

Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2009, Nature, 458, 607

Aguilar, M., Aisa, D., Alvino, A., et al. 2014, Phys. Rev. Lett., 113, 121102

Aharonian, F., An, Q., Axikegu, Bai, L. X., et al. 2021, Phys. Rev. Lett., 126, 241103

Amato, E., & Blasi, P. 2009, MNRAS, 392, 1591

Bai, X.-N. 2022, ApJ, 928, 112

Bai, X.-N., Caprioli, D., Sironi, L., & Spitkovsky, A. 2015, ApJ, 809, 55

Bai, X.-N., Ostriker, E. C., Plotnikov, I., & Stone, J. M. 2019, ApJ, 876, 60

Bao, L.-Z., Fang, K., Bi, X.-J., & Wang, S.-H. 2022, ApJ, 936, 183

Bowen, T. A., Chandran, B. D. G., Squire, J., et al. 2022, Phys. Rev. Lett., 129, 165101

Caraveo, P. A., Bignami, G. F., De Luca, A., et al. 2003, Science, 301, 1345

Cerri, S. S. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202449492

Cerri, S. S., Califano, F., Jenko, F., Told, D., & Rincon, F. 2016, ApJ, 822, L12

Cerri, S. S., Servidio, S., & Califano, F. 2017, ApJ, 846, L18

De La Torre Luque, P., Fornieri, O., & Linden, T. 2022, Phys. Rev. D, 106, 123033
Denton, R. E., & Kotschenreuther, M. 1995, J. Computat. Phys., 119, 283

Derouillat, J., Beck, A., Pérez, F., et al. 2018, Comput. Phys. Commun., 222, 351

Dieckmann, M. E., Alejo, A., & Sarri, G. 2018, Phys. Plasmas, 25, 062122

Dieckmann, M. E., Folini, D., Hotz, I., et al. 2019, A&A, 621, A142

Dieckmann, M. E., Folini, D., Walder, R., Charlet, A., & Marcowith, A. 2022, Phys. Plasmas, 29, 092103

Di Mauro, M., Manconi, S., & Donato, F. 2019, Phys. Rev. D, 100, 123015

Di Mauro, M., Manconi, S., Negro, M., & Donato, F. 2021, Phys. Rev. D, 104, 103002

Engelbrecht, N. E., & Strauss, R. D. T. 2018, ApJ, 856, 159

Evoli, C., Linden, T., & Morlino, G. 2018, Phys. Rev. D, 98, 063017

Evoli, C., Amato, E., Blasi, P., & Aloisio, R. 2021, Phys. Rev. D, 103, 083010

Fang, K., Bi, X.-J., Yin, P.-F., & Yuan, Q. 2018, ApJ, 863, 30

Fang, K., Bi, X.-J., & Yin, P.-F. 2019, MNRAS, 488, 4074

Foote, E. A., & Kulsrud, R. M. 1979, ApJ, 233, 302

Gaensler, B. M., & Slane, P. O. 2006, ARA&A, 44, 17

Gupta, S., Caprioli, D., & Haggerty, C. C. 2021, ApJ, 923, 208

Haverkorn, M., Brown, J. C., Gaensler, B. M., & McClure-Griffiths, N. M. 2008, ApJ, 680, 362

Henri, P., Sgattoni, A., Briand, C., Amiranoff, F., & Riconda, C. 2019, J. Geophys. Res. (Space Phys.), 124, 1475

H.E.S.S. Collaboration (Aharonian, F., et al.) 2023, A&A, 673, A148

Holcomb, C., & Spitkovsky, A. 2019, ApJ, 882, 3

Jóhannesson, G., Porter, T. A., & Moskalenko, I. V. 2019, ApJ, 879, 91

Kasaba, Y., Matsumoto, H., & Omura, Y. 2001, J. Geophys. Res., 106, 18693

Krafft, C., & Savoini, P. 2021, ApJ, 917, L23

Kulsrud, R. M. 2005, Plasma Physics for Astrophysics (Princeton University Press)

Lee, S.-Y., Yoon, P. H., Lee, E., & Tu, W. 2022, ApJ, 924, 36

López-Coto, R., & Giacinti, G. 2018, MNRAS, 479, 4526

Manconi, S., Di Mauro, M., & Donato, F. 2020, Phys. Rev. D, 102, 023015

Marcowith, A., van Marle, A. J., & Plotnikov, I. 2021, Phys. Plasmas, 28, 080601

Marret, A., Ciardi, A., Smets, R., & Fuchs, J. 2021, MNRAS, 500, 2302

Martin, P., Marcowith, A., & Tibaldo, L. 2022a, A&A, 665, A132

Martin, P., Tibaldo, L., Marcowith, A., & Abdollahi, S. 2022b, A&A, 666, A7

Mignone, A., Bodo, G., Vaidya, B., & Mattia, G. 2018, ApJ, 859, 13

Mukhopadhyay, P., & Linden, T. 2022, Phys. Rev. D, 105, 123008

Plotnikov, I., Ostriker, E. C., & Bai, X.-N. 2021, ApJ, 914, 3

Posselt, B., Pavlov, G. G., Slane, P. O., et al. 2017, ApJ, 835, 66

Profumo, S., Reynoso-Cordova, J., Kaaz, N., & Silverman, M. 2018, Phys. Rev. D, 97, 123008

Recchia, S., Di Mauro, M., Aharonian, F. A., et al. 2021, Phys. Rev. D, 104, 123017

Schroer, B., Evoli, C., & Blasi, P. 2023, Phys. Rev. D, 107, 123020

Strauss, R. D., & le Roux, J. A. 2019, ApJ, 872, 125

Tang, X., & Piran, T. 2019, MNRAS, 484, 3491

Torres, D. F., Cillis, A., Martín, J., & de Oña Wilhelmi, E. 2014, J. High Energy Astrophys., 1, 31

Vainio, R. 2000, ApJS, 131, 519

van Marle, A. J., Casse, F., & Marcowith, A. 2018, MNRAS, 473, 3394

Woodham, L. D., Wicks, R. T., Verscharen, D., & Owen, C. J. 2018, ApJ, 856, 49

Zweibel, E. G., & Everett, J. E. 2010, ApJ, 709, 1412
<figure><img src="image_10.png" /><figcaption>**Fig. A.1.** Same as Figure 1 but for the simulation with electron-positron background $m_i = m_e$; run PIC-17. Overall time evolution of field intensities (panel a) and drift velocities of different populations (panel b).</figcaption></figure>

## Appendix A: Electron-positron background plasma

We want to test whether the asymmetry in the dynamic behaviour of CR positron population with respect to the CR electron population (see Figure 1) and if the spectrum asymmetry is due to the ionic nature of the background plasma. In this appendix, we simplify the simulation setup by considering background plasma being composed of $e^{\pm}$ pairs ($m_i = m_e$) instead of ion-electron case studied in the main text. Indeed, such a setup imposes mass and charge symmetry that minimises space charge effects.

Due to typical scales change and the fact that the background mass is no longer dominated by ions we also changed several other parameters as follows: $m_i = m_e$, $n_{cr}/n_i = 10^{-3}$, $T_{cr} = m_e c^2$; $\Omega_{ce}/\omega_{pe} = 0.01$. The initial drift was unchanged $V_D/c = 0.1$, and background temperature $T_{bg} = 10^{-3} m_e c^2$. The deduced Alfvén speed is $V_A = \Omega_{ci}/\omega_p = \Omega_{ce}/\omega_p = 10^{-2}/\sqrt{2}$, where the division by $\sqrt{2}$ is due to $m_i = m_e$, hence total plasma frequency is $\omega_p = \sqrt{2}\omega_{pe}$. These parameters are also reported in Table 2, simulation PIC-17. This setup ensures that the energy density in CRs (kinetic + rest mass) is not dominant compared to the background.

In Figure A.1 we present the time evolution of wave energy averaged over the simulation box and drift velocities of four populations (background electrons and positrons, and CR electrons and positrons) $^9$. Despite a high level of noise floor measured in the longitudinal electric field $E_x$, we can clearly identify the linear phase of the instability between $t\Omega_{ce} \approx 1000$ and $t\Omega_{ce} = 1600$. The saturation occurs after $t\Omega_{ce} = 2500$. The intensity of self-generated waves settles at $\delta B_{sat}^2/B_0^2 \approx 0.45$. As seen in the bottom panel, the decrease in CR drift is identical for CR electrons and CR positrons. Both drop at the same rate from ~14$V_A$ to $V_A$, with good accuracy. Concerning the background plasma, it gains some amount of momentum corresponding to final drift $\langle V_{bg} \rangle_{\text{pos}} = \langle V_{bg} \rangle_e \approx 0.055V_A$.

Concerning the wave spectrum and anisotropy in CR distribution function, in Figure A.2 they are presented in the same format as in Figure 2 of the main text. Both branches of wave polarisation are identical at any time of the simulation. Their intensity at three different times is presented with blue and red lines in the top row of the figure. The same is true for the anisotropy in the distribution function (bottom row of the figure). Both CR electrons and positrons undergo a reduction of anisotropy at the same rate. It also reflects in identical behaviour of overall drift velocities. The final state presents a roughly flat $f(V_x)$, such as the anisotropy being basically erased. We also note that the line-like feature in the high-$k$ region of the wave spectrum, observed in simulations with $m_i \gg m_e$, is absent here.

To complete the analysis of this simulation, we report the derived linear growth rate of the instability in Figure A.3. The agreement between simulation and analytical calculation is good up to $kR_{L,cr} \approx 10$. As in the Fid simulation, we had to use the full dispersion equation that does not assume $\omega \ll \Omega_{ce}$. The test-particle limit solution deviates clearly from the measured rate, especially around the fastest-growing wavenumber. We note that waves are not growing at $kR_{L,cr} > 10$. We recall that the theoretical growth rates were derived in approximation of cold background plasma. But when background thermal effects are considered, the Doppler-shifted cyclotron resonance with thermal electrons and positrons removes waves with $kR_{L,cr} > 20$ ($kR_{L,cr} > 7$) if the intersection between $\omega_r$ and $\omega = \Omega_{ce} - k\nu_{th,e}$ ($\omega = \Omega_{ce} - 3k\nu_{th,e}$) is considered. In Figure A.3 this effect can be seen as an intersection between the black solid line ($\omega_r$) and black dot-dashed line (dotted line, respectively). This explains the absence of self-generated waves in $kR_{L,cr} > 7$ region observed in Figure A.2. At saturation, it seems that the resonance condition with factor of 3, $\omega = \Omega_{ce} - 3k\nu_{th,e}$, provides better agreement for the measured cut-off location in the simulation, consistent with what was found in the main text for the cut-off of $I_-(k)$, see Sect. 3.2.4.

We conclude that the different dynamical behaviour of CR electrons and CR positrons, the asymmetry in the self-generated wave spectrum and line-like emission feature are absent if the background is $e^{\pm}$ plasma. All these features require the background to be of ionic nature.

## Appendix B: Linear growth rates: mono-energetic and power-law cases

For completeness, we provide below the growth rates in the quasi-linear limit for power-law and mono-energetic distributions. These two distributions can physically occur in realistic astrophysical systems and could be explored in future numerical

<footnote>$^9$ This figure is basically the same as Figure 1 but for the $e^{\pm}$ background simulation.</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. A.2. Wave spectrum $I(k)$ (top panels), magnetic field spatial profile (middle panels) and distribution function of beam electrons and positrons $f(V_x)$ and different times of the simulation. Left column corresponds to early time, middle column to the linear phase and right column to the saturated phase of the instability.</figcaption></figure>

simulations. Expressions given below could be useful as analytic expectations for such study.

## B.1. Mono-energetic distribution

In the case the beam has a mono-energetic distribution we have $g(p) = p_0^{-2}\delta(p - p_0)$. It comes

$$ \chi_{cr} = i \left( 1 - \frac{kV_D}{\omega} \right) \frac{2\pi^2 e^2 n_{cr}}{k\omega p_0}. \qquad (B.1) $$

## B.2. Power-law distribution

In the case the beam has a power-law distribution we have for the dimensionless function $g(p) = (\alpha - 3)p_0^{-3}(p/p_0)^{-\alpha}$ (Evoli et al. 2018), where $p_0 \le p \le p_{max}$. It comes if $\alpha > 2$:

$$ \chi_{cr} \approx i \left( 1 - \frac{kV_D}{\omega} \right) \frac{2\pi^2 e^2 n_{cr}}{k\omega p_0} \frac{(\alpha - 3)}{(\alpha - 2)} \left( \frac{p_{\min}(k)}{p_0} \right)^{2-\alpha}. \qquad (B.2) $$

# Appendix C: Wave decomposition in Fourier space

Here we describe the procedure of construction of the wave intensity in Fourier space, $I(k)$. We follow the same procedure as presented in Bai et al. (2019) (Appendix A), but without assuming that the waves are strictly Alfvénic. The main lines are described below.

The wave intensity $I(k)$ is constructed according to

$$ \int_{k_{\min}}^{k_{\max}} I(k)dk = \left\langle \frac{\delta B^2}{B_0^2} \right\rangle. \qquad (C.1) $$

The brackets $\langle \cdot \rangle$ stand for spatial average. A way to decompose circularly polarised modes is to consider the following combination:

$$ C_{\pm}(x) = \frac{\delta B_y \pm i \delta B_z}{B_0}. \qquad (C.2) $$

A discrete Fourier transform of $C^{\pm}(x)$,

$$ W_{\pm}(k_i) = \frac{1}{N} \sum_{n=0}^{N-1} C^{\pm}(x_n) e^{-ik_i x_n}, \qquad (C.3) $$

gives the wave amplitude and phase. The wavenumber of mode $i$ is defined as $k_i = \frac{2\pi}{L_x}i$ and $L_x = N\Delta$ is the length of the simulation box. The wave intensity is then given by

$$ I_{\pm}(k) = |W_{\pm}(k)|^2 \frac{L_x}{2\pi}. \qquad (C.4) $$

With this construction, one defines the dimensionless power in modes of wavenumber $k$ per logarithmic scale as $F(k) = kI(k)$, since $dk = 2\pi/L_x$. For clarity, we mainly used $I(k)$ in the figures of the main part of the article.
<figure><img src="image_12.png" /></figure>

**Fig. A.3.** Real (black solid line) and imaginary part (dashed lines) of the dispersion relation obtained by solving Equation (16), in the case where the background is pair plasma. The growth rates from the corresponding simulation are plotted using blue and red solid lines for right- and left-handed modes, respectively. The fit was not performed on $kR_{L,cr} > 50$ modes as being strongly affected by thermal noise. The black dot-dashed and dotted lines correspond to Doppler-shifted cyclotron resonance of background thermal plasma: $\omega = \Omega_{ce} - kv_{th,e}$ and $\omega = \Omega_{ce} - 3kv_{th,e}$, respectively.