Nonlinear Feedback of the Electrostatic Instability on the Blazar-induced Pair Beam and GeV Cascade

Relativistic pair beams produced in the cosmic voids by TeV gamma-rays from blazars are expected to produce a detectable GeV-scale cascade that is missing in the observations. The suppression of this secondary cascade implies either the deflection of the pair beam by intergalactic magnetic fields or, alternatively, an energy loss of the beam due to the beam-plasma instability. Here, we study how the beam-plasma instability feeds back on the beam, using a realistic two-dimensional beam distribution. We find that the instability broadens the beam opening angles significantly without any significant energy loss, thus confirming a recent feedback study on a simplified one-dimensional beam distribution. However, narrowing diffusion feedback of the beam particles with Lorentz factors less than 106 might become relevant, even though initially it is negligible. Finally, when considering the continuous creation of TeV pairs, we find that the beam distribution and the wave spectrum reach a new quasi-steady state, in which the scattering of beam particles persists and the beam opening angle may increase by a factor of hundreds. Understanding the implications on the GeV cascade emission requires accounting for inverse-Compton cooling.


INTRODUCTION
Blazars are active galactic nuclei with their relativistic jet pointing toward Earth.Observations of the Fermi-LAT telescope and the imaging atmospheric Cerenkov telescopes (such as VERITAS, MAGIC, and HESS) show bright GeV-TeV γ-ray emission from several blazars.During their propagation through the intergalactic medium (IGM), those very high energy γrays interact with the extragalactic background light (EBL), producing a focused beam of electron-positron pairs, that are anticipated to dissipate their energies via inverse Compton scattering on the cosmic microwave background (CMB) (Gould & Schréder 1967;Blumenthal & Gould 1970).mahmoud.s.a.alawashra@uni-potsdam.de martin.pohl@desy.deAlthough primary γ-rays with energies of a few TeV would initiate an electromagnetic cascade in the GeV energy range, such emissions seem to be absent from the γray spectra of certain blazars (Neronov & Semikoz 2009) and possibly the isotropic γ-ray background (Blanco et al. 2023).One possible explanation for the absence of the GeV cascade emission from the γ-ray spectra of blazars is deflection of the TeV pairs by the intergalactic magnetic fields (IGMF) (Elyiv et al. 2009;Neronov & Semikoz 2009;Neronov & Vovk 2010;Taylor et al. 2011;Takahashi et al. 2011;Vovk et al. 2012;Durrer & Neronov 2013).This deflection results in an extended emission or/and a time delay of the cascade emission.In this case, the observed blazar spectra are used to put lower limits on the strength of the IGMF.
The only alternative solution for the missing GeV cascade emission within the standard model of physics is beam energy loss by collective beam-plasma instabilities that is faster than inverse Compton cooling on the CMB.However, whether the non-linear evolution of those instabilities is efficient in taking a significant fraction of the beam energy is still debated (Broderick et al. 2012;Schlickeiser et al. 2012;Miniati & Elyiv 2013;Schlickeiser et al. 2013;Broderick et al. 2014;Sironi & Giannios 2014;Chang et al. 2014;Supsar & Schlickeiser 2014;Chang et al. 2016;Kempf et al. 2016;Rafighi et al. 2017;Vafin et al. 2018Vafin et al. , 2019;;Alves Batista et al. 2019;Shalaby et al. 2020;Perry & Lyubarsky 2021;Alawashra & Pohl 2022).Vafin et al. (2018) calculated the linear growth rate of the electrostatic instability using a realistic pair distribution generated by the annihilation of highenergy gamma rays with the extragalactic background light.Their results demonstrated that the finite angular spread of the beam plays a decisive role in shaping the unstable electrostatic modes.Specifically, their results revealed that the fastest growth rates occurred for wave vectors that are quasi-parallel to the beam direction, while growth rates at oblique directions are smaller compared to the peak values.
Previous studies of the blazar-induced pair beam instabilities didn't consider the instability feedback on the pair beam particles.Perry & Lyubarsky (2021) studied this feedback for the first time in the context of blazarinduced pair beam electrostatic instability.Their findings imply that the back reaction of the electrostatic unstable waves on the pair beam widens the beam opening angles by around one order of magnitude without any significant energy loss.
In this study, unlike the simplified one-dimensional beam distribution used in Perry & Lyubarsky (2021), we use a two-dimensional realistic beam distribution to explore the influence of the instability feedback on the beam.Specifically, we use the beam profile at a distance of 50 Mpc from the blazar found in Vafin et al. (2018).This treatment enables us to look at the feedback influence on the pairs that have the relevant Lorentz factors for cascade emission in the GeV band.
The instability feedback is described as Fokker-Planck diffusion both in momentum and angular space.This treatment was simplified in the analysis by Perry & Lyubarsky (2021), by evaluating only the initially dominant angular widening diffusion and neglecting the other effects involving the momentum diffusion and angular narrowing (D θp ).Here, we check rigorously this assumption by using the 2D spectrum of the expanded beam under the dominant feedback to analyse the possible impact of the momentum diffusion on the beam energy and whether the beam narrowing is still negligible.
The blazar-induced pair beam-plasma instability significantly outpaces other factors that could change the beam profile, such as inverse Compton cooling and pair production.Whereas previous works have predominantly focused on assessing the instability's impact on a stationary beam profile, we incorporate the continuous production of TeV pairs into the transport equation of the beam, in addition to the diffusion terms.
The structure of this paper is as follows.In section 2, we introduce the pair beam realistic 2D profile that we used in this study.In section 3, we present the quasilinear theory of the beam-plasma system, we introduce the linear growth rate of the electrostatic instability, the time evolution of the electrostatic waves spectrum and the Fokker-Planck diffusion of the beam distribution.Finally, we demonstrate the numerical simulation and the results in section 4 and conclude in section 5.

BLAZAR-INDUCED PAIR BEAM DISTRIBUTION
The pair-beam distribution function is the crucial quantity that determines the beam-plasma instability growth rate (Vafin et al. 2018).Thus, using the realistic spectrum of blazar-induced pair beams is essential for examining the influence of the beam-plasma instability on the beam and the GeV-scale cascade emission.In this study, we used the realistic beam distribution at a distance of 50 Mpc from the blazar, as reported in Vafin et al. (2018).Here, we introduce this beam spectrum and explain the ingredients used to find it.
The propagation of the beam distribution in the IGM is driven by two primary factors.The first one is the pair's production due to the interaction of the highenergy gamma rays with the EBL, along with their subsequent cooling processes.The second effect is the dispersion of the primary gamma-ray flux with the propagation distance, leading to an inverse proportionality of the beam density with the square of the distance from the blazar.These two fundamental mechanisms collectively shape the evolution of the pair-beam distribution along the propagation distance in the IGM.
The two effects have been combined in Vafin et al. (2018), neglecting the IC cooling, to calculate the accumulated pair spectrum over the IC cooling mean-free path of pairs with Lorentz factor of 10 7 starting at the distance 50 Mpc from the blazar.The neglect of IC cooling is driven by the necessity to investigate beamplasma instabilities that provide the dominant energy loss.They used an intrinsic power-law gamma-ray spectrum with a spectral index of 1.8 and a cut-off step function at the energy of 50 TeV.The normalized initial beam distribution 2πp 3 θ 2 f (p, θ)/n b at distance 50 Mpc from the blazar (Vafin et al. 2018).
We define the normalized beam momentum distribution where n b is the pair-beam density.The density at 50 Mpc from the blazar is estimated as Vafin et al. 2018).Factorizing the distribution as where f cos θ (γ, θ) is an angular differential part, f γ (γ) is a momentum differential part, and β is the normalized speed.The angular part is approximated by a Gaussian (Broderick et al. 2012;Miniati & Elyiv 2013;Vafin et al. 2018) with the angular spread of ∆θ = 1 γ .The momentum part, f γ (γ), is given by eq.B14 in the appendix B, where we replaced the sharp step-function cut-off used in Vafin et al. (2018) by an exponential cut-off at the Lorentz factors higher than 6 × 10 6 using a part of a logarithmic Gaussian as shown in Fig. 12.
The initial normalized realistic beam spectrum is shown in Fig. 1, and the main energy bulk of the pair is located at Lorentz factors of a few 10 6 .We also see that the pairs are concentrated in a narrow band around the production angles of θ ∼ γ −1 .We will use this distribution to find the linear growth rate of the instability in section 3.1 and as the initial condition for the Fokker-Planck simulation of the instability feedback in section 4.

QUASILINEAR THEORY OF THE BEAM-PLASMA SYSTEM
The beam-plasma instabilities manifest in both electrostatic and electromagnetic modes, including the twostream instability (k × δE = 0 where δE is the perturbed electric field), the transverse Weibel, and filamentation modes (k • δE = 0) (Bret et al. 2010).The electrostatic modes dominate the wave spectrum for the blazar-induced TeV beams, whereas Weibel-type modes are suppressed (Bret et al. 2005;Rafighi et al. 2017).Consequently, we consider only the electrostatic oblique modes, which is sufficient to recover the essential physics (Chang et al. 2016).
In section 3.1, we present the linear growth rate of the electrostatic instability.In section 3.2, we introduce the balance equation for waves.Lastly, in section 3.3, we introduce the Fokker-Planck diffusion equation describing the instability feedback on the beam distribution.

Electrostatic linear growth rate
In the kinetic regime that is applicable for blazarinduced pair beams (Miniati & Elyiv 2013), the linear growth rate of an unstable wave with a wave vector k can be found by (Breizman 1990) where ω p = (4πn e e 2 /m e ) 1/2 is the plasma frequency of the intergalactic background plasma with density n e .
Here, we neglected the influence of intergalactic medium temperature, T e , on the plasma frequency, ω, expressed as ω 2 = ω 2 p + 3k 2 v 2 th,e ≈ ω 2 p , where v th,e is the IGM electros' thermal speed.For an intergalactic temperature of 10 4 K, the plasma frequency experiences a negligible shift of around 2.6 × 10 −6 (ω = ω p (1 + 2.6 × 10 −6 )), insignificantly affecting linear growth rate calculations.However, it's important to mention that intergalactic medium temperature may significantly impact the nonlinear evolution of waves, as discussed in studies of nonlinear landau damping by (Miniati & Elyiv 2013;Chang et al. 2014;Vafin et al. 2019).
Exploiting the cylindrical symmetry around the beam propagation axis (z-axis in our case), we fixed the wave vector of the electrostatic waves to k = (k ⊥ , 0, k || ), where k ⊥ and k || are the perpendicular and the parallel components to the beam propagation direction respectively.After integrating over the azimuthal angle of the beam we get the following where the boundaries are Here θ ′ is the wave vector angle with the beam propagation direction (z-axis), θ is the angle between momentum and the beam axis, and v b ≃ c(1 − 1 2γ 2 ) is the particle speed.
In Fig. 2, we present the linear growth rate (eq.5 -6), using the beam distribution introduced in section 2, the density of IGM electron as n e = 10 −7 (1 + z) 3 cm −3 , and a redshift z = 0.15.Previous treatments in the literature used the approximation (v b = c) when evaluating eq.5 (Miniati & Elyiv 2013;Vafin et al. 2018;Perry & Lyubarsky 2021).However, we found that in the regimes of wave numbers with (ck ωp − 1 , the difference between the particle speed and the speed of light becomes relevant.We have taken this difference into account in our calculations of Fig. 2. In Fig. 2, we see that the growth rate is maximal and constant in the range of perpendicular wave numbers, 10 −6 < ck ⊥ /ω p < 1, with a sharp drop at the oblique angles ck ⊥ /ω p ∼ 1, whereas for the parallel modes, k ⊥ c/ω p < 10 −6 , it is smaller by around a factor of 3. Note that the growth rate of the parallel modes is sensitive to the beam distribution, for a simplified monoenergetic Gaussian the growth of the parallel modes is larger than the quasi-parallel ones (Perry & Lyubarsky 2021), where it's smaller for the distribution we used and for a Maxwell-Jüttner distribution (Chang et al. 2016).
The turnover of the linear growth rate spectral shape around the wave numbers of (ck ωp − 1 is due to the change in the corresponding resonant beam angles.In the regime of (ck ⊥ /ω p ) 2 >> ck || ωp − 1 , the resonant beam angles are constrained by the minimum angles of θ 1 = (ck || /ω p − 1)/(ck ⊥ /ω p ) and large θ 2 .In the regime of (ck ⊥ /ω p ) 2 << ck || ωp − 1 , the resonant angles boundaries approach each other to a ck ⊥ /ω p in- . The maximum growth rate, ω i,max ∼ 6.7 × 10 −8 s −1 , is much faster than the IC cooling rate of the beam, However, the instability-induced energy-loss rate significantly depends on the nonlinear evolution of the instability (Miniati & Elyiv 2013;Schlickeiser et al. 2013;Chang et al. 2014;Vafin et al. 2019).
In this study we focus on the instability feedback, therefore we will consider only the linear regime of the instability and neglect the restrictions on the growth of the waves due to non-linear interactions.In the next section, we briefly introduce the linear evolution equation of unstable waves and levels of the unstable wave's energy density where the nonlinear processes become relevant.

Evolution of the wave spectrum
The quasi-linear evolution of the wave spectrum for homogeneous plasma is governed by the following equation where W (k) is the spectral energy density of the electric field oscillations, ω i (k) is the linear growth rate as defined in section 3.1, and ω c is the collisional damping rate (Tigik et al. 2019), Here g = (n e λ 3 D ) −1 is the plasma parameter, λ D = 6.9 cm Te/K ne/cm −3 is the Debye length, n e = 10 −7 (1 + z) 3 cm −3 is the density of IGM electrons, and T e = 10 4 K is their temperature.We start integrating eq.( 7) at the very low thermal fluctuations level.
The collisional damping rate given by eq.8 is approximately 20 times smaller than the approximation employed in other studies (i.e.Miniati & Elyiv (2013); Vafin et al. (2019); Perry & Lyubarsky (2021)), which did not account for the microscopic wave-particle interactions.Those interactions were included under the generalized weak turbulence theory in Yoon et al. (2016), deriving an accurate general kinetic formulation of the collisional damping rate of the electrostatic plasma waves, that was used in Tigik et al. (2019) to find the collisional damping rate (eq.8).
The total electric field energy density is calculated by Accounting for the energy equipartition between kinetic electrostatic fluctuations, the energy loss rate of the beam due to the growth of the electrostatic waves at time t is given by (Vafin et al. 2018) where the total beam energy density is defined as In the previous section, we found that the modes with the maximum growth, 10 −6 < ck ⊥ /ω p < 1, grow at the same rate (Fig. 2), maintaining a similar spectral amplitude.However, the energy density in those modes is proportional to their wave number volume element, 2πk ⊥ ∆k ⊥ ∆k || (eq.9).Therefore, we can focus on the quasi-parallel and oblique modes, 10 −3 < ck ⊥ /ω p < 1, since they dominate the energy density of the unstable mode spectrum.We can also neglect inhomogeneity of the background plasma since it is relevant only for the strictly parallel modes (Perry & Lyubarsky 2021;Shalaby et al. 2020).
The amplitude of the unstable modes grows exponentially until their wave intensity is high enough to trigger nonlinear processes.One of the main non-linear processes is the modulation instability that moves wave energy from resonant to non-resonant modes.This process operates when the total electric field energy density hits the threshold of (kλ D ) 2 n e T e (Miniati & Elyiv 2013), resulting in the saturation of the resonant unstable mode at around 10 −3 of the total beam energy we consider here.
Another non-linear process is non-linear Landau damping, where the non-linear scattering of the unstable waves on the background plasma ions results in severe damping of the resonant modes.Non-linear Landau damping becomes effective when the total electric field energy density reaches around 10 −2 of that of the beam (Chang et al. 2014;Vafin et al. 2019).
The impact of these non-linear interactions is still uncertain (Schlickeiser et al. 2012;Miniati & Elyiv 2013;Chang et al. 2014;Vafin et al. 2019).The numerical noise in simulations (such as PIC) is too high, and the numerical growth rate is too small, for a reliable assessment, on account of the very small beam density.Upscaling of the beam density and downscaling the beam Lorentz factor is possible, but the results of those simulations are difficult to scale back to the realistic parameters (Sironi & Giannios 2014;Rafighi et al. 2017).
In this work, we focus on the nonlinear feedback of the instability on the beam and so we consider only the linear phase of the instability growth.We discuss in section 4 that under the instability feedback on the beam, the total electric field energy density stays always below the non-linear thresholds for the beam density we consider.In the next section, we look at the Fokker-Planck diffusion equation that describes the feedback of the electrostatic waves on the beam during the quasilinear regime.

Fokker-Planck diffusion equation for the pair beam
The quasilinear regime is applicable when the total wave energy density is much smaller than that of the plasma.In this regime, the feedback of the electrostatic unstable waves on the beam is governed by the following Fokker-Planck diffusion equation (Brejzman & Ryutov 1974) where the diffusion coefficients are defined by the following resonance integrals (Rudakov 1971) where the electric charge, e, is given in cgs units.The pair-beam distribution function, f , is given in spherical coordinates (p, θ, φ), and so is the wave-vector k (k, θ ′ , φ ′ ).The angles θ and θ ′ are defined with respect to the beam propagation direction (z −axis).Due to the azimuthal symmetry of the pair-beam distribution function, we can set φ = 0 and integrate over φ ′ , yielding (see appendix A) where and v b = c(1 − 1 2γ 2 ) is the particle speed for Lorentz factor γ. The boundaries of the cos θ ′ integration are fixed by the resonance condition The integrands are largest at the peak of the wave spectrum, therefore a proper numerical resolution of the spectrum is necessary when calculating the diffusion coefficients.We have changed the integration variables in appendix A arriving at the coordinates (k ⊥ , θ R ) with ωp − 1 /(ck ⊥ /ω p ), for which the peak of the unstable modes is numerically well resolved and the diffusion coefficients are well defined by eq.A11.We also found that the Lorentz factor, γ dependence of the diffusion coefficients is negligible compared to the beam angle, θ.
In the next section, we describe the numerical setup of our study of the instability feedback and present our results.

NUMERICAL RESULTS
We have calculated the rate of change for every term on the right-hand side of the Fokker-Planck equation (eq.12), using the diffusion coefficients of a wave spectrum generated by the growth rate presented in section 3.1.We found that the diffusion term D θθ exceeds the other terms by orders of magnitude in the phase-space region containing the bulk of the beam particles.Evaluating the maximum rate across the entire parameter space for every term, we found the following ratio between the different terms: θθ : θp : pθ : pp ≈ 1 : 10 −3 : 10 −5 : 10 −8 .
Given this result, we initially neglect all the subdominant terms and in section 4.1 consider only the diffusion  term D θθ .We will check the validity of this approximation in section 4.2 as we analyse the effect of the subdominant terms as the θθ diffusion modifies the beam.We also analyse the dependence of our results on the beam parameters in section 4.3.Finally, we add the continuous pair production to our simulation setup in section 4.4.

Simulation of the θθ angular diffusion feedback
Having established that D θθ initially dominates over the other diffusion terms by orders of magnitudes, we perform here a numerical simulation of instability feedback including only this term.We introduce the simulation setup in section 4.1.1 and present the results in section 4.1.2.

Simulation setup
The first numerical simulation of the beam-plasma system only includes the first term on the right-hand side of eq.12 coupled time-dependently with the waves' spectral evolution equation (eq.7).The linear growth rate of the instability (eq.5) and the diffusion coefficients (eq.14) involve integration over the beam distribution function and the wave spectrum, respectively.We solve eq.17 using the Crank-Nicolson scheme along with the FTCS scheme for the wave equation, eq.7.We used a dynamical time step of ω −1 i,max as the default time step with an upper limit set by the fastest rate of change of the distribution.We tested this by using time steps that are 10 times smaller.In order to properly resolve the narrow wave spectrum we use a logarithmic grid in the coordinates (ck ⊥ /ω p , θ R ) where θ R = (ck || /ω p − 1)/(ck ⊥ /ω p ).We verified convergence in our grid resolution for both the wave spectrum and the beam distribution.The initial beam distribution is as described in section 2, and the initial wave energy density corresponds to the fluctuation level (Vafin et al. 2019).

Results
We found that the instability feedback severely increased the beam's angular spread.This broadening strongly depends on the Lorentz factor of the beam particles.In Fig. 3, we show the angular spread for different beam Lorentz factors.We see that particles with larger Lorentz factors get scattered earlier since those particles are in resonance with faster-growing wave modes, and so the scattering feedback affects them earlier.
The angular spreading of the beam immediately shifts the resonant wave numbers.In Fig. 4, we see the reduction of the growth rate for the parallel wave numbers during the simulation time for a fixed perpendicular wave number of 0.1.This reduction starts at the fastest growing modes as they quickly scatter their resonant particles, and with time it extends to slower growing modes at higher parallel wave numbers.We found  that the initial profile of the linear growth with respect to the perpendicular wave numbers, as shown in Fig. 2, doesn't change during the time evolution.
The resulting wave spectrum of the time-dependent linear growth rate is shown in Fig. 5.In the beginning, the fastest-growing modes form a spectral peak.Once the beam widens, the slower modes at higher parallel wave numbers start forming a second peak until the wave's intensity is sufficient to kick the beam particles to higher angles.The process keeps repeating until the linear growth rate becomes less than or comparable to the collisional damping rate (presented by the dashed black line in Fig. 4).By the time we stop the simulation at 5 × 10 12 seconds, all modes are collisionally damped.
In Fig. 6, we show the diffusion coefficient, D θθ , at various times.The variation of the diffusion coefficient with beam angle, θ, closely resembles that of the wave spectrum with parallel wave numbers, as represented in Fig. 5.This is due to the resonance relation between the beam angle and the parallel wave numbers (θ = (ck || /ω p − 1)/(ck ⊥ /ω p ) in the regime ωp − 1 ).Peaks in the wave spectrum (Fig. 5) and diffusion coefficient (Fig. 6) arise from the evolving linear growth rate (Fig. 4).The initial peak forms at wave numbers with the highest growth rates, causing resonant particles to diffuse, smoothly shifting the peak.A second peak emerges as particles diffuse further, resonating with higher wave numbers.Eventually, collisional damping leads to the decay of both the first and second peaks.A third peak follows in the same mechanism, decaying at a lower amplitude due to the fall of the growth rate below the collisional damping shortly.
The appearance of the last peak in the diffusion coefficient profile, manifesting at angular values around 10 −5 radians, causes the observed surge in the beam angular spread after 7 × 10 11 seconds (Fig. 3).The time scale of this upturn in angular spread is governed by the interplay between the evolving linear growth rate and the collisional damping rate.We also see that around this time pairs with different Lorentz factors react to the same resonant unstable modes, because in the regime (ck ⊥ /ω p ) 2 >> ck || ωp − 1 waves with a certain parallel wave number are resonant with the beam particles at a certain angle, whatever their momentum.
We see in Fig. 3 that by the time the instability has saturated, the angular spread of pairs with Lorentz factor 10 6 has increased by around two orders of magnitudes, much more than the factor of ten reported by Perry & Lyubarsky (2021).The main reasons for this higher spread are the smaller collisional damping rate and the higher beam density we used.
We found that during the entire simulation time, the wave energy density never exceeded 10 −3 of the beam energy density.This level of the wave intensity is lower than that needed for efficient operation of nonlinear Landau damping and the Modulation instability (Vafin et al. 2019;Chang et al. 2014;Miniati & Elyiv 2013).Therefore, the effect of these non-linear processes on the instability development might be minimal compared to that of the diffusive feedback on the beam.We also calculated the total energy transferred from the beam to the waves by integrating the energy loss rate of the beam given in eq.10 over time.The result is given by the black dashed line in Fig. 7.We see that the beam lost less than 1% of its total initial energy by the time the instability development was saturated by the widening feedback.Those results suggest that the feedback widening severely limits the energy transfer from the beam to the waves.We explore whether this situation changes as we use different beam densities in section 4.3.
Up to here, we only included the initially dominant term D θθ of the right-hand side of eq.12.In the next section, we analyse the feedback of the other subdominant terms as the dominant θθ diffusion widens the pair beam.

2D analysis of the diffusion equation
We analyse here the effect of the subdominant terms as the beam widens, using the time-dependent beam distribution that we numerically derived and discussed in the previous section.
For the momentum diffusion of the beam (third and fourth terms on the RHS of eq.12), we can calculate the energy loss or gain rate of the beam by inserting the corresponding time derivative of the beam distribution eq.12 in the total rate of change of the beam energy.After integrating by parts we get the following relation for pθ diffusion  eq.20).The diffusion θp dominates over the term (θθ) in the orange and red areas with values higher than zero, while it contributes less than 10% in the dark blue area.The drop in the ratio just before the rim at 10 12 seconds is due to the increase in the widening as a result of wave growth outside the initial resonance region.After 10 12 seconds, the collisional damping effectively damps the waves, and the impact of both terms declines.Looking at the overall sign of eq.18, we see that the diffusion pθ involves a global energy loss of the beam since D pθ and the angular derivative are always negative.We also found that the dominant feedback of the diffusion pp is an energy gain of the beam, as D pp is always positive and the beam distribution function declines for γ ≳ 10 5 (see Fig. 12).
Integrating eq.18 and eq.19 over the simulation time and dividing by the total beam initial energy, we see in Fig. 7 the accumulated fraction of the beam energy lost and gained.We observe that pθ diffusion could eliminate only around 0.1% of the beam total energy by the end of the simulation whereas pp diffusion increases it by a negligible fraction.Therefore, it is evident that the cumulative effect of the diffusive momentum flux on the beam is insignificant compared to the scattering.Now, we proceed to the analysis of the second term on the RHS of eq.12, θp diffusion.This diffusion involves angular flux as the θθ diffusion, but it can result in both the narrowing and widening of the beam depending on the beam momentum gradient.Pairs with negative momentum gradient, γ > 10 5 , experience narrowing whereas the ones with positive momentum gradient, γ < 10 5 , experience a widening.In Fig. 8, we have compared the normalized angular integral of the absolute rate of change of the θθ diffusion for a certain beam Lorentz factor with that one of the θp diffusion (21) It is noticeable in Fig. 8 that the ratio of I θp and I θθ increases gradually until it drops after 7 × 10 11 s.The reason for the increase is that the diffusive flux of θθ decreases as the beam profile flattens, while the diffusive flux of θp remains relatively constant as the momentum gradients are not impacted by the beam broadening.The drop after 7 × 10 11 seconds is due to the increase of diffusion θθ by the accumulated wave density outside the initial resonance region that we discussed in the previous section.
In Fig. 8 we see that the θp diffusion becomes dominant for Lorentz factors less than 10 6 at times much earlier than their inverse Compton cooling time (≳ 10 13 s).For these particles including this diffusion is necessary.However, there is a minimal impact of the θp diffusion on the pairs that are capable of giving IC emission in the detectable GeV band (Lorentz factors of 10 6 or slightly higher).This indicates that θp diffusion might not impact the GeV-scale cascade emission as strongly as the θθ diffusion does.

Parameters dependence
In the simulation discussed in section 4.1, we used a fiducial pair beam density at a distance of 50 Mpc from the blazar, 3 × 10 −22 cm −3 (Vafin et al. 2018).However, the beam density changes under different conditions, such as varying the distance from the source, changing the source's luminosity, or using different EBL models in the calculations.Here we vary the beam density using the same setup as in section 4.1 and investigate its impact on our results.
In Fig. 9, we see the fraction of the beam energy lost by the instability for different beam densities.As the beam density is increased, the instability develops earlier and takes more energy from the beam.However, the beam lost only 2% even for a very high beam density, 8 × 10 −21 cm −3 .Therefore, the fundamental physical behaviour of the system remains consistent, the beam experiences expansion with a negligible energy loss of its initial energy as the instability is saturated by the beam expansion.
We noticed that the wave intensity during the simulations with beam density higher than 10 −21 cm −3 has exceeded the threshold for the non-linear modulation instability, which we didn't include in our calculations.Those processes will impose further restrictions on the growth of the unstable modes.
We observed that the angular spread increased by a factor of 1.5 when the beam density was inflated by a factor of three.This scaling can be attributed to the fact that the linear growth rate is linearly proportional to the pair beam density and inversely proportional to the square of the beam angular spread, ω i ∝ n b ∆θ 2 .Therefore, increasing the beam density by a factor C requires an increase in the angular spread by a factor of approximately √ C to maintain the reduction of the linear growth rate to the collisional damping rate at the time when the instability has saturated.
In the remainder of this section, we will discuss the influence of the cut-off energy in the intrinsic gammaray spectrum on the results of section 4.1.Vafin et al. (2018) used an intrinsic power-law gamma-ray spectrum with a step function cut-off at the energy of 50 TeV.However, in the end, they used the attenuated gammaray spectrum at a distance of 50 Mpc to calculate the accumulated pair beam spectrum over a certain path length.At a distance of 50 Mpc from the blazar, the majority of gamma rays with energies higher than 10 TeV have already been absorbed.The mean free path for a gamma-ray with energy E γ for pair production with the EBL photons is given by where ξ = 4.5 and ξ = 0 for redshifts of z ≤ 1 and z > 1, respectively (Kneiske, T. M. et al. 2004;Neronov & Semikoz 2009).Therefore, any cut-off energy above the 10-TeV threshold will have only a minimal impact.

Simulation with injection
In section 4.1, we found that the instability growth is severely reduced by beam broadening to the point that it cannot be isolated from the production and the cooling rates of the beam.In this section, we include the pair creation rate in the evolution equation of the beam.
The beam distribution found in Vafin et al. (2018) was calculated as the accumulation of pairs over the path length of 7.7 × 10 12 light-seconds, using a constant production rate Q ee .We added this production rate to the beam evolution equation along with the dominant θθ diffusion term, Using the same simulation setup as described in section 4.1, we numerically solved the coupled system of the evolution equations (eq.23 and eq.7).For times much less than 10 13 seconds we found essentially the same behaviour of the system as without injection.After 10 13 s, a new quasi-steady state of the beam distribution and the waves spectrum emerges.
The creation of highly focused pairs with beam angles of the order γ −1 increases the linear growth rate at wave numbers in resonance with these particles.Ultimately, this leads to a quasi-equilibrium of the wave spectrum and the beam distribution.On the wave side, the linear growth rate and the collisional damping rate balance across the resonant wave numbers, resulting in a steady-state wave spectrum as shown in Fig. 10.On the beam side, the diffusive scattering compensates the pair production, keeping the beam expanding as shown in the angular profile of pairs with a Lorentz factor of 10 6 in Fig. 11.This ongoing expansion of the beam extends the unstable modes to higher parallel wave numbers as shown in Fig. 10 We have stopped the simulation after 5 × 10 13 s, which corresponds to the IC cooling time of pairs with a Lorentz factor of 10 6 .By this time the pairs have experienced a diffusive deflection up to angles of around 4 × 10 −4 radians.This deflection results in an arrival time delay of the secondary GeV-band photons emitted  by those pairs (Neronov & Semikoz 2009).The arrival time delay of secondary gamma rays emitted by pairs that have undergone a deflection by an angle of ∆θ from the primary gamma-ray propagation direction is given by the following formula where D c is the distance between the emitting pairs and the blazar and D b is the distance between the blazar and the Earth.Given that for our simulation setup D c = 50 Mpc and D b = 720 Mpc for the fiducial z = 0.15, the formula reduces to ∆t delay = ∆θ 2 × 7.6 × 10 7 years.Hence, the deflection of pairs with Lorentz factor of 10 6 by 4 × 10 −4 radians implies a time delay of around 10 years for the GeV-scale cascade emission produced at the distance 50 Mpc from the source.Calculating the deflection at different distances from the source is needed to find the impact on the observed cascade emission.This is beyond the scope of this paper and will be covered in future works.
As of the previous simulation in section 4.1, we also found here that the wave energy density never exceeded 10 −3 of the total beam energy density, keeping the wave intensity at levels lower than what is needed for the nonlinear processes to operate efficiently, again justifying their neglect.This paper's calculations did not consider the IC scattering of the beam particles.For particles with Lorentz factors γ < mec 2 4ϵCMB ∼ 2.5 × 10 8 , IC scattering occurs in the Thomson regime, leading to momentum loss without significant angular changes.While we simulate the beam until the IC cooling time of pairs with γ = 10 6 , we anticipate a substantial decrease in the beam steady state for Lorentz factors above 10 6 .The steady-state wave spectrum is expected to undergo less significant changes however as it is influenced by the resonance condition that is set by the particles' angle irrespective of their energy.Nevertheless, the inclusion of IC cooling is crucial for a comprehensive understanding of the physical implications on the arrival time distribution of GeV-scale cascade emissions, and it is part of our future research plans.

CONCLUSIONS
We have explored the feedback of the electrostatic beam-plasma instability on blazar-induced pair beams.The feedback of the beam-plasma instability is described by a Fokker-Planck diffusion equation with diffusion coefficients that are dependent on the resonance condition between the unstable waves and the beam particles.This feedback is crucial for understanding the propagation of the blazar-induced pair beam in the IGM, in particular the question whether or not the instability is capable of draining the beam energy faster than inverse Compton cooling.Such insights hold significance in unravelling the underlying reasons for the absence of secondary GeV-scale emissions in several distant blazar spectra (Neronov & Semikoz 2009;Broderick et al. 2012).
We solved the Fokker-Planck diffusion equation for the beam distribution function, coupled with the linear evolution equation of the plasma-wave spectrum.As the initial condition for the beam, we used the realistic two-dimensional beam distribution computed by Vafin et al. (2018) for a distance of 50 Mpc from the blazar.Initially, the dominant feedback is angular broadening of the beam, stemming from the scattering of the beam particles by the excited waves.As the instability widens the beam, the instability growth rate is severely reduced, leading at the end to a negligible energy transfer from the beam to the plasma waves.These findings align with a recent study on instability feedback (Perry & Lyubarsky 2021).
Using the 2D time-dependent beam profile evolving by the predominant angular diffusion, we found that momentum diffusion does not have any significant impact on the beam.However, we found that another angular diffusion term, which is initially negligible, might become relevant and may narrow the beam particles with Lorentz factors between 10 5 and 10 6 .Therefore, including this term in the feedback calculations is necessary for a comprehensive understanding of the instability impact on those pairs.However, the GeV-scale cascade is emitted by pairs with Lorentz factors of 10 6 or slightly higher, and so the impact of this term on the GeV-scale secondary cascade might be limited.
In our analysis, we neglected non-linear wave interactions in the evolution of the wave spectrum.For beam density lower than 10 −21 cm −3 , we found that the cumulative energy density of the electric field fluctuations remains below the critical thresholds required to trigger the significant impacts of the non-linear processes, such as non-linear Landau damping or the modulation instability.However, for higher beam densities the wave energy density exceeded the threshold for the non-linear modulation instability.Those non-linear processes would impose further restrictions on the growth of the unstable modes.
Lastly, we have included the continuous TeV pairs production in the Fokker-Planck diffusion equation.Unlike the previous simulation discussed in section 4.1, in this particular configuration, the unstable modes do not decay after the beam has expanded but saturate at a finite amplitude.The wave spectrum reaches a quasiequilibrium across the wave numbers resonant with the beam injection angles.The beam particles experience persistent scattering under the diffusive feedback of this steady-state wave spectrum.Then, beam particles with Lorentz factors of 10 6 scatter up to angles of around 4 × 10 −4 radians within their IC cooling time.This results in a time delay of around 10 years in the arrival of the secondary GeV-scale cascade, assuming pairs at a distance of 50 Mpc from a blazar that is 720 Mpc away from Earth.We expect that this estimate depends on the beam density that varies along the propagation distance and with source luminosity.
In the end, calculating the broadening at more points along the beam propagation is needed to understand the impact of the instability broadening on the GeV-scale cascade emission.Also, it's essential to include the inverse Compton cooling in the beam distribution evolution equation to understand the long-term time evolution of the beam-wave system.
After transforming the delta function we get Table 1.The parameters for the approximation in eq.B14.

Figure 2 .
Figure 2. Normalized linear growth rate using the realistic beam distribution function.White areas denote stable modes.

Figure 3 .
Figure 3.The angular spread for different Lorentz factors of the beam as a function of time during the angular diffusion feedback simulation presented in section 4.1.The turn-up around the time of 7 × 10 11 seconds is due to the growth of the wave spectrum's third peak as seen in Fig.5.

Figure 4 .
Figure 4. Evolution of the linear growth rate of the instability for a fixed perpendicular wave number during the angular diffusion feedback simulation presented in section 4.1.The black dashed line represents the collisional damping rate.Legend values are common logarithms of time in seconds.Throughout the simulation, we observed that the linear growth rate has maintained its initial profile with perpendicular wave numbers as in Fig.2.

Figure 5 .
Figure 5.The time evolution of the wave spectrum for fixed perpendicular wave number for the angular diffusion feedback simulation presented in section 4.1.The formation of the peaks here is due to the spectral change of the linear growth rate with time as shown in Fig.4.Legend values are common logarithms of time in seconds.

Figure 6 .
Figure 6.The D θθ for γ = 10 6 as a function of time in the angular diffusion feedback simulation presented in section 4.1.Legend values are common logarithms of time in seconds.

Figure 7 .
Figure 7.The accumulated change in the beam energy during the angular diffusion feedback simulation presented in section 4.1.The black dashed line (∆W ) represents the beam energy fraction going into unstable wave growth.The dashed cyan line (∆ pθ ) and the dashed red line (∆pp) represent the fraction of the beam energy loss and gain due to the momentum diffusion by the pθ and the pp terms, respectively.

Figure 8 .
Figure8.The logarithm of the ratio of I θp (eq.21) and I θθ (eq.20).The diffusion θp dominates over the term (θθ) in the orange and red areas with values higher than zero, while it contributes less than 10% in the dark blue area.The drop in the ratio just before the rim at 10 12 seconds is due to the increase in the widening as a result of wave growth outside the initial resonance region.After 10 12 seconds, the collisional damping effectively damps the waves, and the impact of both terms declines.

Figure 9 .
Figure 9.The accumulated fraction of the beam energy lost as a function of time due to the wave growth during the angular diffusion simulations feedback with different values of the beam density.All the values are in units of 10 −22 cm −3 .

Figure 10 .
Figure 10.The wave spectrum as a function of time during the injection simulation presented in section 4.4.Legend values are common logarithms of time in seconds.We see the steady-state wave spectrum emerging after 10 13 s.

Figure 11 .
Figure 11.The angular profile of pairs with Lorentz factor of 10 6 during the injection simulation presented in section 4.4.Legend values are common logarithms of time in seconds.We see that around the IC cooling time of 5 × 10 13 s for those pairs, they have been deflected by around 4 × 10 −4 radians which yields a time delay of around 10 years for the GeV cascade.