© The Authors 2024

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 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 upcoming 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 ballistic 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 emission requires a power input that exceeds the available spin-down luminosity in some cases (Bao et al. 2022). Alternatively, 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 number 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 origin of this reduced coherence length is that the pulsar still finds itself in a medium under the influence of the associated supernova remnant turbulence (Fang et al. 2019). Finally, Evoli et al. (2018) investigate the triggering of the resonant streaming instability 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 including 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 distances 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 thousands of years later, at which time the pulsar has lost much of its spin-down power (Gaensler & Slane 2006); for (ii) and (iii), realistic 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 expansion of the remnant¹, 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². 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 highest. 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 simulations. 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 Sec. 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 Sec. 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 composed 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 (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 .$$(1)

The plasma response function (or susceptibility) χ has two components, one for the background χ_bg and one for the CR beam χ_cr. 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_D with respect to the background plasma assumed to be cold. We first calculate the CR contribution χ_cr. Assuming an isotropic distribution in the momentum space g(p) given by f_cr = n_crg(p)/4π we have $$\begin{array}{rl}{\chi }_{\text{cr }}=& \sum _{\alpha }\frac{\pi q_{\alpha }^2}{\omega }{n}_{\text{cr },\alpha }{\int }_{p_0}^{p_{max}}\mathrm{d}p{p}^2\upsilon (p)\frac{\mathrm{d}{g}_{\alpha }(p)}{\mathrm{d}p}\\ & \times {\int }_{-1}^1\frac{(1-{\mu }^2)}{\omega +k\upsilon (p)\mu \pm {\mathrm{\Omega }}_{\alpha }}\mathrm{d}\mu ,\end{array}$$(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_∥ = k in the above expression. Each particle species has a distribution g_α(p) such that ∫ dpp²g_α(p) = 1; different forms of particle distribution are discussed below. Finally, n_cr,α is the CR density. If the calculation is carried in the background plasma frame, a supplementary term (1 – kV_D/ω) 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: $$\begin{array}{rl}{\int }_{-1}^1\frac{(1-{\mu }^2)}{\omega +k\upsilon \mu \pm {\mathrm{\Omega }}_{\alpha }}\mathrm{d}\mu =& \mathcal{P}{\int }_{-1}^1\frac{(1-{\mu }^2)}{\omega +k\upsilon \mu \pm {\mathrm{\Omega }}_{\alpha }}d\mu \\ & -i\pi {\int }_{-1}^1(1-{\mu }^2)\delta (\omega +k\upsilon \mu \pm {\mathrm{\Omega }}_{\alpha })d\mu ,\end{array}$$(3)

where the particle speed is a function of the momentum p and can be written as υ(p). We find: $$\begin{array}{rl}{\chi }_{\mathrm{c}\mathrm{r}}=& \frac{\pi e^2}{k\omega }\sum _{\alpha }{N}_{CR,\alpha }{\int }_{p_0}^{p_{max}}\mathrm{d}p{p}^2\frac{\mathrm{d}{g}_{\alpha }(p)}{\mathrm{d}p}\\ & [(1-{\mu }_{\mathrm{R},\alpha }^2)\log \left|\frac{1-{\mu }_{\mathrm{R},\alpha }^{-1}}{1+{\mu }_{\mathrm{R},\alpha }^{-1}}\right|-2{\mu }_{\mathrm{R},\alpha }-i\pi (1-{\mu }_{\mathrm{R},\alpha }^2)],\end{array}$$(4)

where the wave-particle resonance selects a particular pitch-angle cosine $${\mu }_{\mathrm{R},\alpha }=\mp \frac{{\mathrm{\Omega }}_{\alpha }}{k\upsilon }-\frac{\omega }{k\upsilon }.$$(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⁻p⁺) 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_D we have: $${\chi }_{\mathrm{b}\mathrm{g}}=-\frac{4\pi e^2}{{\omega }^2}(\frac{n_{\mathrm{i}}}{m_{\mathrm{i}}}\frac{\omega +k{V}_{\mathrm{D}}}{\omega +k{V}_{\mathrm{D}}\pm {\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}}+\frac{n_{\mathrm{e}}}{m_{\mathrm{e}}}\frac{\omega +k{V}_{\mathrm{D}}}{\omega +k{V}_{\mathrm{D}}\pm {\mathrm{\Omega }}_{\mathrm{c}\mathrm{e}}}).$$(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_D) ≪ Ω_ci ≪ Ω_ce. Hence, if we have n_e = n_i, the contribution from the background e⁻p⁺plasma is $${\chi }_{\mathrm{b}\mathrm{g}}\simeq \frac{{(\omega +k{V}_{\mathrm{D}})}^2}{{\omega }^2}\frac{c^2}{V_{\mathrm{A}}^2}.$$(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 Sec. 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 ω ≪ Ω_e+. The non-resonant term vanishes and we recover the results by Evoli et al. (2018) (their Eq. (17)), namely: $${\chi }_{\mathrm{c}\mathrm{r}}=-2\frac{i{\pi }^2{e}^2}{k\omega }{n}_{\mathrm{c}\mathrm{r}}{\int }_{p_{min}}^{p_{max}}\mathrm{d}p\frac{\mathrm{d}g(p)}{\mathrm{d}p}(p^2-p_{min}(k{)}^2),$$(8)

where the maximum particle momentum is p_max and $p_{\text{res }}(k)=\frac{m_{\mathrm{e}}\left|{\mathrm{\Omega }}_{\mathrm{c}\mathrm{e}}\right|}{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 ⋍ pc ≫ m_ec². The MJ distribution is then defined as: $$g(p)=\frac{1}{m_{\mathrm{e}}^3{c}^3}\frac{1}{\theta K_2\left(\frac{1}{\theta }\right)}{e}^{-\gamma (p)/\theta },$$(9)

where θ = k_BT_cr/(m_ec²), k_B is the Boltzmann constant, T_cr is the temperature of CR electrons (same for CR positrons), and K₂ is the modified Bessel function of the second kind. The particle Lorentz factor is γ(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{array}{rl}{\chi }_{\mathrm{c}\mathrm{r}}& \simeq i\frac{4{\pi }^2{e}^2}{\omega k}\frac{n_{\mathrm{c}\mathrm{r}}}{m_{\mathrm{e}}c}\frac{\theta }{K_2(1/\theta )}I(k,{\gamma }_{min}),\\ I(k,{\gamma }_{min}(k))& =\exp (-\frac{{\gamma }_{min}}{\theta })(\frac{{\gamma }_{min}(k)}{\theta }+1),\\ {\gamma }_{min}(k)& =\sqrt{1+{\left(\frac{{\mathrm{\Omega }}_{\mathrm{c}\mathrm{e}}}{kc}\right)}^2},\end{array}$$(10)

and γ_min = γ(p_min). We have assumed that γ(p_max) ≫ θ. The dispersion relation then reads as $$k^2{V}_{\mathrm{A}}^2={\omega }^2{\left(\frac{V_{\mathrm{A}}}{c}\right)}^2+{(\omega +k{V}_{\mathrm{D}})}^2+i\pi \frac{\omega {\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}}{k{R}_{\mathrm{L},\mathrm{e}}}\frac{n_{\mathrm{c}\mathrm{r}}}{n_{\mathrm{i}}}\frac{\theta }{K_2(1/\theta )}I(k,{\gamma }_{\text{min }}),$$(11)

where R_L,e = c/Ω_ce = and R_L,cr = k_BT_cr/eB. The instability growth rate is given by the system $$k^2{V}_{\mathrm{A}}^2\simeq {\tilde{\omega }}^2+i\omega \frac{n_{\mathrm{c}\mathrm{r}}}{n_{\mathrm{i}}}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}{I}_2(k),$$(12)

with $\tilde{\omega }=\omega +k{V}_{\mathrm{D}},I_2(k,\theta )=\frac{\pi }{k{R}_{\mathrm{L},\mathrm{e}}}\frac{\theta }{K_2(1/\theta )}{e}^{-\frac{{\gamma }_{min}}{\theta }}(1+\frac{{\gamma }_{min}}{\theta })$. The system is of second order as we assume $\mathfrak{R}(\tilde{\omega })\ll k{V}_{\mathrm{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}_{\mathrm{D}}\frac{n_{\mathrm{c}\mathrm{r}}}{n_{\mathrm{i}}}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}$, the linear growth rate $\mathrm{\Gamma }\equiv \mathfrak{I}(\tilde{\omega })$ reads as $$\mathrm{\Gamma }\simeq \frac{\pi }{2}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}\frac{1}{k{R}_{\mathrm{L},\mathrm{c}\mathrm{r}}}\frac{V_{\mathrm{D}}}{V_{\mathrm{A}}}\frac{n_{\mathrm{c}\mathrm{r}}}{n_{\mathrm{i}}}\frac{{\theta }^2}{K_2(1/\theta )}{e}^{-\frac{{\gamma }_{min}}{\theta }}(1+\frac{{\gamma }_{min}}{\theta }).$$(13)

and the dispersion relation is almost unmodified $\mathfrak{R}(\tilde{\omega })=k{V}_{\mathrm{a}}$. In the relativistic regime γ_min(k) ⋍ θ/kR_L,cr) and θ²/K₂(1/θ) → 0.5 (for θ ≫ 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 (ω + kV_D) ≪ Ω_ci. In this case, processes at frequencies Ω_ci < ω < Ω_ce are now included (e.g., ion-cyclotron resonances, excitation of Whistler waves). Yet, the two following approximations remain valid: Ω_ci ≪ Ω_ce and ω ≪ ω_pe. The former is always satisfied for m_i ≫ 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 }_{\mathrm{b}\mathrm{g}}=\frac{c^2}{V_{\mathrm{A}}^2}\frac{{\tilde{\omega }}^2}{{\omega }^2}\frac{1-{\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}/{\mathrm{\Omega }}_{\mathrm{c}\mathrm{e}}}{[1\pm \tilde{\omega }/{\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}][1\pm \tilde{\omega }/{\mathrm{\Omega }}_{\mathrm{c}\mathrm{e}}]},$$(14)

leads to a third-order dispersion relation $$k^2{V}_{\mathrm{A}}^2={\omega }^2\frac{V_{\mathrm{A}}^2}{c^2}+\frac{{\tilde{\omega }}^2}{[1\pm \tilde{\omega }/{\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}][1\pm \tilde{\omega }/{\mathrm{\Omega }}_{\mathrm{c}\mathrm{e}}]}+\frac{i\omega C}{k{V}_{\mathrm{D}}},$$(15)

where $C\equiv \pi V_{\mathrm{D}}\frac{{\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}}{R_{\mathrm{L},\mathrm{c}\mathrm{r}}}\frac{n_{\mathrm{c}\mathrm{r}}}{n_{\mathrm{p}}}\frac{{\theta }^2}{K_2\left(\frac{1}{\theta }\right)}I(k,{\gamma }_{min})$. After reorganising the terms $$\begin{array}{rl}& {\tilde{\omega }}^3\frac{iC}{k{V}_{\mathrm{D}}}+{\tilde{\omega }}^2({\mathrm{\Omega }}_{\mathrm{c}\mathrm{e}}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}-{\left(k{V}_{\mathrm{A}}\right)}^2-iC\pm \frac{iC{\mathrm{\Omega }}_{\mathrm{c}\mathrm{e}}}{k{V}_{\mathrm{D}}})\\ & +\tilde{\omega }{\mathrm{\Omega }}_{\mathrm{c}\mathrm{e}}(\mp {\left(k{V}_{\mathrm{A}}\right)}^2\mp iC+\frac{iC{\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}}{k{V}_{\mathrm{D}}})\\ & -{\mathrm{\Omega }}_{\mathrm{c}\mathrm{e}}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}[{\left(k{V}_{\mathrm{A}}\right)}^2+iC]=0.\end{array}$$(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 ΔP_cr ⋍ (4/3)[(V_D – V_A)/c]ncr⟨p⟩ with the momentum density acquired by forward propagating Alfvén waves $\mathrm{\Delta }{P}_{\mathrm{w}}\simeq \rho V_{\mathrm{A}}\left(\delta B^2/B_0^2\right)$. Here, n_cr is the CR density and ⟨p⟩ the mean CR momentum, ρ is the local gas (thermal proton) density. We find: $$\frac{\delta B_{\text{sat }}^2}{B_0^2}\simeq \frac{8}{3}\frac{(V_{\mathrm{D}}-V_{\mathrm{A}})}{V_{\mathrm{A}}}\frac{n_{\mathrm{c}\mathrm{r}}}{n_{\mathrm{i}}}\frac{m_{\mathrm{e}}}{m_{\mathrm{i}}}\frac{⟨\gamma v⟩}{c}.$$(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 ∥ B₀, see Sec. 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 ≪ 1. The external magnetic field is parallel to the simulation box B = B₀e_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/ω_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/ω_pe, is resolved with 6 grid points. The time step is Δt = 0.95Δ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³. 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 ≪ 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⁻³m_ec², that is θ_bg = 10⁻³ (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⁻³ and 10⁻². 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⁻².

• The initial distribution function of the beam is drifting Maxwell–Jüttner distribution (see Eq. (9)) with temperature of T_cr (typically θ = 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ω_pe = 10⁶, hat generally corresponds to tΩ_ci = 10³ (for the case of Fid run: m_i/m_e = 100 and ω_pe/Ω_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 \simeq 8\pi n_{\mathrm{i}}{k}_{\mathrm{B}}{T}_{\mathrm{i}}/B_0^2={\left({\upsilon }_{\mathrm{t}\mathrm{h},\mathrm{i}}/V_{\mathrm{A}}\right)}^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, β = 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 (δf / f₀ ≪ 1) and is equivalent to an average drift velocity in x-direction that is small compared to c, V_D ≪ 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, 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Ω_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 ⋍ 2V_A (see orange and magenta lines in the bottom panel). At some point (around tΩ_ci = 150) the wave growth ends, and the instability enters the saturated phase. During the saturated phase (at tΩ_ci > 150) the evolution becomes slow as the isotropization process of remaining particles takes more time. At the final time of the simulation (tΩ_ci = 10³) the level of magnetic fluctuations is $\delta B^2/B_0^2=(\delta B_y^2+\delta B_z^2)/B_0^2\simeq 0.06$, while the drift velocity of beam positrons decreased to ~ 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 ⟨V_i⟩_bg ⋍ ⟨V_e⟩_bg ⋍ 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 δB_y,z/B₀ (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Ω_ci = 7.6), in the middle of the linear phase (tΩ_ci = 76), and the late post-saturation time (tΩ_ci = 186), respectively. The wave intensity in Fourier space for different polarisations I_±(k) was done by transforming real-space magnetic fields δB_x,y and according to $\int I(k)dk=⟨\delta B^2/B_0^2⟩$. Details of the method are presented in Appendix C.

During the initial time of the evolution, at tΩ_ci ~ 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 μ = 0 (i.e., near V_x = 0). In the middle of the linear phase of the instability, at tΩ_ci = 76.0, wave-modes with 0.1 ≤ kR_L,0 ≤ 10 are rapidly growing. The fastest growth occurs around kR_L,0 ≈ 0.7 as expected from the analytical calculation in Eq. (13). Panel e shows that the wave fluctuations reach δB_y,z/B₀ ⋍ 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 ≤ kR_L,0 ≤ 20. The amplitude of spatial fluctuations reached δB/B₀ ⋍ 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 ⋍ 170 or kd_e ⋍ 1.7 (panels d and g). While being energetically subdominant to the bulk of excited wavemodes at kR_L,0 ≤ 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).

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.

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 (Δt_lin) and then perform the linear fit of $|\log \left({\tilde{B}}_{y,z}(k,\mathrm{\Delta }{t}_{\mathrm{lin}})\right)|$. This procedure provides the numerical growth rate of different wavemodes, available in the simulation. The smallest wavemode is k_min = 2π/L and the largest is k_max = π/Δ, where Δ is the spatial grid cell size⁴.

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 ⋍ 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 ~ 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.

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 δBy,z/B0 (middle panels: b, e, and h) and distribution function of beam electrons and positrons f(Vx) and different times of the simulation (bottom panels: c, f, and i). Left column corresponds to early time (tΩci = 7.6), middle column to the linear phase (tΩci = 76) and right column to the saturated phase of the instability (tΩci = 186).

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_{\text{sat }}^2=\delta B_{y,\text{end }}^2+\delta B_{z\text{.end }}^2$.

In Fig. 4, we show the dependence of saturation level of the magnetic field, $\delta B_{\mathrm{s}\mathrm{a}\mathrm{t}}^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)⁵. 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⁶.

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 × 10⁻³ and n_cr/n_i = 10⁻² 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, ${⟨V_i⟩}_{\mathrm{b}\mathrm{g}}/V_{\mathrm{A}}\simeq \delta B_{\mathrm{s}\mathrm{a}\mathrm{t}}^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 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.

Fig. 3 Linear growth rate of the instability Γ 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 ncr/ni = 0.01, middle panel corresponds to ncr/ni = 5 × 10−3, and bottom panel corresponds to ncr/ni = 2 × 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).

Fig. 4 Saturation level of the magnetic fluctuations. Panel a: wave saturation intensity as function of ncr/ni. 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 ncr/ni = 10−3 is a lower limit as the simulation ended before a saturated state was reached. Panel b: same but varying mi/me at fixed ncr/ni = 10−2 (blue symbols) and ncr/ni = 5 × 10−3 (red symbols).

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_{\text{sat }}(k)=\delta B_{\text{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 ~ 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 ≥ 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 ℛ-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 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 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 ⋍ 6, while for m_i/m_e = 1000 it is located at kR_L,cr ⋍ 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⁻² and red circles are from simulations with n_cr/n_i = 5 × 10⁻³. All other parameters are the same as in Fid run. The black dashed line follows k_cutd_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 ω = Ω_ci − 3kυ_th,i (red dotted line in the same panel)⁷. This gives the cut-off at k_cutd_i = 0.47, 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_{\mathrm{c}\mathrm{u}\mathrm{t},-}=\frac{{\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}}{V_{\mathrm{A}}+3{\upsilon }_{\mathrm{t}\mathrm{h},\mathrm{i}}}.$$(18)

Left-handed waves with k > k_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_cut,−. Numerical solutions of the dispersion relation accounting for thermal effects in background plasma – using plasma dispersion function – shows as well that k_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_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 ω = Ω_ce − 3kυ_th,e occurs at kd_e ⋍ 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 ⋍ 0.1 (i.e., kR_L,cr ⋍ 10), see top panel of Fig. 5. This value is sensitively lower than expected k_cut,+. By inspecting in more detail, we found that there are indeed growing waves up to kd_e ~ 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 ≤ 5 × 10⁻³) 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_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.

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 ω (solid lines) and instability growth rate for the two polarizations that is the imaginary part of ω (dashed lines). The light-wave is shown with black solid line. The red star symbol at the top left side locates the intersection with ω = 2ωpe line. Alfvén wave corresponds to low frequency and low-k part of red and blue solid lines. At kRL,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 kdi = 1) for left-handed waves, and electron-cyclotron regime kde = 1 for right-handed waves. The lower red star symbol locates the intersection between the ion-cyclotron branch (red solid line) and Ωce – 3kυth,i (red dotted line).

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 mi/me = 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 mi/me. Blue circles correspond to values derived from simulations with ncr/ni = 10−2 and red circles are from simulations with ncr/ni = 5 × 10−3. The dashed black line follows kcutdi = 0.5.

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 ⋍ 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 ω(k) from the simulation. This was done by storing the EM fields on the spatial grid at every time-step between tΩ_ci = 95 ant tΩ_ci = 96.6. It corresponds to the 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 ω(k).

In Fig. 7, we present the ω(k) dispersion spectrum performed on B_z component of the magnetic field. The light-wave (EM wave) is well-seen with dispersion ${\omega }^2\simeq k^2{c}^2+{\omega }_{\mathrm{p}\mathrm{e}}^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 cutoff frequency at ω_rh ⋍ ω_pe + Ω_ce/2 while the left-handed branch has the cut-off frequency at ω_LH ⋍ ω_pe − Ω_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 ω = 2ω_pe. In terms of the wavenumber, theoretically, it corresponds to $k{d}_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_{\mathrm{L},\mathrm{c}\mathrm{r}}^{-1}$. In this panel, one can identify the dominant contribution of kd_e < 0.1 waves and a much less intense region slightly below $k{d}_{\mathrm{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 $k{d}_{\mathrm{e}}\simeq \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 ω < Ω_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 ω/ω_pe < 0.1 and kd_e < 0.1. It is not resolved in this figure as the time window for diagnostics was too small.

Fig. 7 Wave dispersion analysis. Top panel a: wave dispersion ω(k)/ωpe at the final times of the linear phase. From figure simulation (see Table 2 for parameters). The diagnostic was performed on the δBz component between tΩci = 95 ant tΩ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 ω < Ωce the Whistler branch is clearly seen as well (theoretical dispersion line over plotted with dashed line). Plasma emission at ω = ω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_{\mathrm{e}}^{-1}$. This improves the visibility of the high-k part of the spectrum.

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 ≫ 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_⊥ related phenomena. Yet, it is known that at kinetic scales with k_⊥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 × 10³⁴ erg/s and an initial spin-down timescale t₀ ⋍ 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_{\text{age }}\frac{{(1+t_{\text{age }}/t_0)}^2}{{(1+t/t_0)}^2}.$$(19)

This implies an initial spin-down luminosity L₀ = 2.9 × 10³⁷ erg/s at early times t ≪ t₀, three orders of magnitude above the present-day value. The particle injection spectrum follows a broken power law⁸ $$\begin{array}{rlrl}Q(E,t)& =Q_0(t){\left(\frac{E}{E_{\mathrm{b}}}\right)}^{-a_{\mathrm{L}}}& & \text{ for }\mathrm{E}\le {\mathrm{E}}_{\mathrm{b}}\\ & =Q_0(t){\left(\frac{E}{E_{\mathrm{b}}}\right)}^{-a_{\mathrm{H}}}& & \text{ for }\mathrm{E}>{\mathrm{E}}_{\mathrm{b}},\end{array}$$

where we use E_b ~ 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 η ⋍ 1 of the pulsar spin-down power, such that the normalisation Q₀ of the injection spectrum is obtained by $${\int }_{E_{min}}^{\mathrm{\infty }}\mathrm{d}EEQ(E,t)=\eta L(t),\eta \le 1.$$(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{array}{rl}⟨E⟩& =H(a_{\mathrm{L}},a_{\mathrm{H}},r)\times E_{\mathrm{b}}\\ & \simeq \frac{{(2-a_{\mathrm{L}})}^{-1}+{(a_{\mathrm{H}}-2)}^{-1}-{(2-a_{\mathrm{L}})}^{-1}{r}^{2-a_{\mathrm{L}}}}{{(a_{\mathrm{H}}-1)}^{-1}-{(a_{\mathrm{L}}-1)}^{-1}+{(a_{\mathrm{L}}-1)}^{-1}{r}^{1-a_{\mathrm{L}}}}\times E_{\mathrm{b}},\end{array}$$(21)

where r = E_min/E_b. For the adopted parameters, this yields H ⋍ 0.12 hence ⟨E⟩ ⋍ 24 GeV. If E_min = 100 MeV then H ⋍ 0.03 and H goes as r^0.6 as r → 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 ⋍ 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}\text{ and so }{R}_{\mathrm{p}\mathrm{w}\mathrm{n}}(t)=R_{\mathrm{p}\mathrm{w}\mathrm{n}}\left(t_{\text{age }}\right)\sqrt{L(t)/L_{\text{age }}}$ and S_PWN ∝ L.

The number density of pairs injected by the nebula at time t can therefore be written as $$n_{\mathrm{c}\mathrm{r}}(t)=\frac{\eta L(t)}{⟨E⟩S_{\mathrm{P}\mathrm{W}\mathrm{N}}{V}_{\mathrm{D}}}.$$(22)

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_α emission surrounding the Geminga bow shock to be in the range between n_i ~ 10⁻³ cm⁻³ (Caraveo et al. 2003) and n_i ~ 10⁻² cm⁻³ (Posselt et al. 2017). We use the latter higher bound together with a magnetic field strength of to compute V_A = 22 km s⁻¹.

According to Eq. (22), the particle density immediately beyond the surface of the nebula does not depend on time because S_PWN ∝ L and all other parameters are considered constant. Using present day values for L = 3.26 × 10³⁴ erg s⁻¹, R_pwn = 0.1 pc and other parameters being specified above (i.e., η = 1, V_D = 10V_A = 220 km s⁻¹, ⟨E⟩ = 24 GeV), we get n_cr ~ 3 × 10⁻⁸ cm⁻³. 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₀ 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_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⁻⁶ − 10⁻⁵ (with the caveat that the background density at early times might have been higher than 10⁻² cm⁻³). 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⁻⁶ − 10⁻⁵ and a beam with a bulk energy of E_b ~ 10⁵m_ec². 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_cr/n_i ≤ 10⁻³ 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 energy of E ~ 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 × 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_cr = 10 m_ec² and a relative particle density of n_cr/n_i = 5 × 10⁻³ for comparison with the PIC simulations. The second simulation has a beam temperature of T_cr = 100 m_ec² and a relative particle density of n_cr/n_i = 10⁻³ to demonstrate MHD-PIC’s ability to reach 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_ec² and = 5 × 10⁻³ (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_ec², the resulting magnetic wave energy evolution in our simulation is shown in Fig. 9. The maximum value $\delta B_{\mathrm{s}\mathrm{a}\mathrm{t}}^2/B_0^2\sim 0.05$ is reached at tΩ_ci ≈ 400. The level of magnetic fluctuations at the saturation state is consistent with the expected value from Eq. (17) which gives $\delta B_{\mathrm{s}\mathrm{a}\mathrm{t}}^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Ω_ci ⋍ 90 until tΩ_ci ⋍ 180. From our simulation, we find a growth rate over that period of approx. Γ/Ω_ci ~ 3 × 10⁻³ 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 δ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.

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.

Fig. 9 Magnetic energy evolution of a MHD-PIC simulation with Tb = 100mec2 and relative beam density ncr/ni = 10−3.

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 ~ 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 ≪ Ω_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 ≫ d_i⁻¹ 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.

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^± 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⁻³, T_cr = m_ec²; Ω_ce/ω_pe = 0.01. The initial drift was unchanged V_D/c = 0.1, and background temperature T_bg = 10⁻³m_ec². The deduced Alfvén speed is $V_{\mathrm{A}}={\mathrm{\Omega }}_{\mathrm{c}\mathrm{i}}/{\omega }_{\mathrm{p}}={\mathrm{\Omega }}_{\mathrm{c}\mathrm{e}}/{\omega }_{\mathrm{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 }_{\mathrm{p}\mathrm{e}}$. 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)⁹. 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Ω_ce ⋍ 1000 and tΩ_ce = 1600. The saturation occurs after tΩ_ce = 2500. The intensity of self-generated waves settles at $\delta B_{\mathrm{s}\mathrm{a}\mathrm{t}}^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 ~ 14V_A to V_A, with good accuracy. Concerning the background plasma, it gains some amount of momentum corresponding to final drift ⟨V_bg⟩_pos = ⟨V_bg⟩_e ⋍ 0.055V_A.

Fig. A.1 Same as Figure 1 but for the simulation with electron–positron background mi = me; run PIC-17. Overall time evolution of field intensities (panel a) and drift velocities of different populations (panel b).

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 ≫ m_e, is absent here.

Fig. A.2 Wave spectrum I(k) (top panels), magnetic field spatial profile (middle panels) and distribution function of beam electrons and positrons f(Vx) 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.

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 ⋍ 10. As in the Fid simulation, we had to use the full dispersion equation that does not assume ω ≪ Ω_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 ω_r and ω = Ω_ce − kυ_th,e (ω = Ω_ce − 3kυ_th,e) is considered. In Figure A.3 this effect can be seen as an intersection between the black solid line (ω_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, ω = Ω_ce − 3kυ_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.

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 kRL,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: ω = Ωce − kυth,e and ω = Ωce − 3kυth,e, respectively.

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^± 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 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 }_{\mathrm{c}\mathrm{r}}=i(1-\frac{k{V}_{\mathrm{D}}}{\omega })\frac{2{\pi }^2{e}^2}{k\omega }\frac{n_{\mathrm{c}\mathrm{r}}}{p_0}.$$(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}{\left(p/p_0\right)}^{-\alpha }$ (Evoli et al. 2018), where p₀ ≤ p ≤ p_max. It comes if α > 2: $${\chi }_{\mathrm{c}\mathrm{r}}\simeq i(1-\frac{k{V}_{\mathrm{D}}}{\omega })\frac{2{\pi }^2{e}^2}{k\omega }\frac{n_{\mathrm{c}\mathrm{r}}}{p_0}\frac{(\alpha -3)}{(\alpha -2)}{\left(\frac{p_{min}(k)}{p_0}\right)}^{2-\alpha }.$$(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=⟨\frac{\delta B^2}{B_0^2}⟩.$$(C.1)

The brackets ⟨·⟩ 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 \mathrm{i}\delta B_z}{B_0}.$$(C.2)

A discrete Fourier transform of C^±(x), $$W_{\pm }\left(k_i\right)=\frac{1}{N}\sum _{n=0}^{N-1}{C}^{\pm }\left(x_n\right)e^{-\mathrm{i}{k}_i{x}_n}\text{, }$$(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Δ 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 }.$$(C.4)

With this construction, one defines the dimensionless power in modes of wavenumber k per logarithmic scale as ℱ(k) = kI(k), since dk = 2π/L_x. For clarity, we mainly used I(k) in the figures of the main part of the article.

References

All Tables

All Figures

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.

In the text

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 δBy,z/B0 (middle panels: b, e, and h) and distribution function of beam electrons and positrons f(Vx) and different times of the simulation (bottom panels: c, f, and i). Left column corresponds to early time (tΩci = 7.6), middle column to the linear phase (tΩci = 76) and right column to the saturated phase of the instability (tΩci = 186).

In the text

Fig. 3 Linear growth rate of the instability Γ 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 ncr/ni = 0.01, middle panel corresponds to ncr/ni = 5 × 10−3, and bottom panel corresponds to ncr/ni = 2 × 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).

In the text

Fig. 4 Saturation level of the magnetic fluctuations. Panel a: wave saturation intensity as function of ncr/ni. 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 ncr/ni = 10−3 is a lower limit as the simulation ended before a saturated state was reached. Panel b: same but varying mi/me at fixed ncr/ni = 10−2 (blue symbols) and ncr/ni = 5 × 10−3 (red symbols).

In the text

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 ω (solid lines) and instability growth rate for the two polarizations that is the imaginary part of ω (dashed lines). The light-wave is shown with black solid line. The red star symbol at the top left side locates the intersection with ω = 2ωpe line. Alfvén wave corresponds to low frequency and low-k part of red and blue solid lines. At kRL,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 kdi = 1) for left-handed waves, and electron-cyclotron regime kde = 1 for right-handed waves. The lower red star symbol locates the intersection between the ion-cyclotron branch (red solid line) and Ωce – 3kυth,i (red dotted line).

In the text

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 mi/me = 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 mi/me. Blue circles correspond to values derived from simulations with ncr/ni = 10−2 and red circles are from simulations with ncr/ni = 5 × 10−3. The dashed black line follows kcutdi = 0.5.

In the text

Fig. 7 Wave dispersion analysis. Top panel a: wave dispersion ω(k)/ωpe at the final times of the linear phase. From figure simulation (see Table 2 for parameters). The diagnostic was performed on the δBz component between tΩci = 95 ant tΩ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 ω < Ωce the Whistler branch is clearly seen as well (theoretical dispersion line over plotted with dashed line). Plasma emission at ω = ω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_{\mathrm{e}}^{-1}$. This improves the visibility of the high-k part of the spectrum.

In the text

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.

In the text

	Fig. 9 Magnetic energy evolution of a MHD-PIC simulation with Tb = 100mec2 and relative beam density ncr/ni = 10−3.
In the text

	Fig. A.1 Same as Figure 1 but for the simulation with electron–positron background mi = me; run PIC-17. Overall time evolution of field intensities (panel a) and drift velocities of different populations (panel b).
In the text

	Fig. A.2 Wave spectrum I(k) (top panels), magnetic field spatial profile (middle panels) and distribution function of beam electrons and positrons f(Vx) 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.
In the text

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 kRL,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: ω = Ωce − kυth,e and ω = Ωce − 3kυth,e, respectively.

In the text