Plasma effects on relativistic pair beams from TeV blazars: PIC simulations and analytical predictions
Iman Rafighi Sergei Vafin Jacek Niemiec Martin Pohl















  





























    B⃗E⃗v⃗x⃗k⃗k⃗^'r⃗k_∥k_⊥J_n(z)J_n^'(z)2J_n^2(z)V⃗f_a^(0)v_∥v_⊥Ωp_∥p_⊥p⃗p_⊥p_∥łαωω _RΩ _aΩ _eω _p,aω ^2_p,aω ^2_p,eΘ _pΘ _eΘu_aβ _∥ , aβ _∥w^2αδ̣ϵ _a,sδ_a,sϵ _nsrϵ _psrβ _c∂ _σaperiodic 
            [ u_a√(1-E^-2)μ _aarctan
                𝕀∋[ κP_e, critR_⊥R_∥PIC simulations and analytical predictionsInstitute for Physics and Astronomy, University of Potsdam, D-14476 Potsdam, Germany
                             
                                 DESY, Platanenallee 6, D-15738 Zeuthen, Germany
                                 
                    Instytut Fizyki Ja̧drowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland
                                 
                    
                     
                      Pair beams produced by very high-energy radiation from TeV blazars emit gamma rays in the GeV band by inverse-Compton scattering of soft photons. The observed GeV-band signal is smaller than that expected from the full electromagnetic cascade. This means that the pair beams must be affected by other physical processes reducing their energy flux. One possible loss mechanism involves beam-plasma instabilities that we consider in the present work. For realistic parameters the pair beams can not be simulated by modern computers. Instead, we use a simple analytical model to find a range of the beam parameters that (i) provides a physical picture similar to that of realistic pair beams and (ii) at the same time can be handled by available computational resources. Afterwards, we performed corresponding 2D PIC simulations. We confirm that the beams experience only small changes in the relevant parameter regime, and other processes such as deflection in magnetic field must be at play.   Plasma effects on relativistic pair beams from TeV blazars: 
                        I. Rafighi,1rafighi@uni-potsdam.de, S. Vafin1, M. Pohl1,2, and J. Niemiec3
                        Accepted August 22, 2017 
                    ==============================================================================
                    
                    
                    
                    
                    
                    § INTRODUCTION
                    
                    
                    Many blazars, a sub-class of active galactic nuclei, have been detected with gamma-ray telescopes such as HESS, VERITAS, Fermi and MAGIC as sources of gamma-rays with the energy E≥ 100 GeV <cit.>. These very-high energy photons interact with extragalactic background light (EBL) producing ultra-relativistic electron-positron pairs with the typical Lorentz factor 10^5<Γ<10^7<cit.>. The created pairs are subject of many investigations, as they can be affected by several physical processes: (i) inverse Compton scattering (ICS), (ii) deflection by the intergalactic magnetic field (IGMF), or (iii) collective plasma effects. The ICS would result in gamma-ray emission with characteristic energy in the GeV band. But, as indicated by Fermi-LAT data, the GeV gamma-ray emission is suppressed meaning that the ICS is not the fastest of the three processes. The effect of deflection by the IGMF has been well investigated <cit.> which led to constraints on the IGMF. However, these constraints are valid only under the assumptions that the multi-TeV gamma-ray emission persists on long timescales and that the pairs lose their energy only due to ICS. The last assumption is very crucial and remains debatable. 
                    
                    The importance of collective plasma effects has been pointed out by several authors <cit.>. In fact, the pairs can induce electrostatic (two-stream, oblique) and electromagnetic (filamentation, Weibel) instabilities <cit.>. In this case, wave-particle interactions can reduce the energy of the pairs by 30-50 %<cit.>. Therefore, the collective plasma effects can also substantially suppress the GeV-band gamma-ray emission affecting as well the IGMF constraints. 
                    
                    The pair beams constitute an extremely small fraction of the plasma density in the intergalactic medium (IGM), α=n_b/n≈ 10^-16-10^-18. This circumstance prohibits direct computer simulations of the beams due to insufficient computational power, and substantial adjustments in parameter values have been made in published simulation studies <cit.>. 
                    At the same time, an accurate analytical description of the non-linear evolution of the plasma system is also problematic. In this work, we combine numerical PIC simulations with a simple analytical model to determine physical parameters of the beam and plasma, so that (i) the problem can be treated with reasonable computational power, and (ii) the physical picture is adequate to realistic pair beams. The physical picture is determined by several aspects: (i) the ratio of the energy densities of the beam and background plasma, (ii) instabilities and their growth rates, and (iii) non-linear damping of plasma waves. Here, we are concerned only with the first two subjects, and the non-linear effects <cit.> will be analyzed in future papers. 
                    Such treatment is possible in the linear stage which we are interested in here.
                    
                    
                    The created pairs are subject to the ICS and a full electromagnetic cascade that modifies their parameters. However, the goal of the present paper is to explore the potential dominance of plasma effects on the beam evolution. Therefore, we consider a pair beam created only by the initial TeV gamma-ray emission neglecting ICS. In this case, the typical parameters of the created beams depend on the distance from a blazar, and they are ⟨Γ⟩=10^5, Γ=10^3-10^8, n_b=10^-25-10^-19 cm^-3, Δθ≈1/⟨Γ⟩≈10^-5 (Δθ is the angular spread), whereas typical parameters of the IGM are T=10^4-10^7 K, n=10^-7 cm^-3<cit.>. Thus, the energy density ratio is ϵ=n_b⟨Γ⟩ m_e c^2/ (n k_B T)≈  10^-10-10^-1 (k_B is the Boltzmann constant, m_e is the electron mass) indicating that the pair beam cannot considerably heat the IGM plasma. This point was realized by <cit.> who conducted simulations for ϵ=0.1. The parameters of the simulations by <cit.> are α=n_b/n≈ 10^-2, Γ≈10^2, and k_BT/(m_e c^2)≈ 10^-8, providing with ϵ≈ 10^8, a parameter regime that is not relevant for realistic pair beams. Moreover, such a high energy-density ratio causes anisotropic plasma heating that can eventually drive the Weibel instability as it will be shown below. Note that Kempf and Sironi have studied a beam distribution with ΔΓ≪⟨Γ⟩. We will also investigate this case in the present work, whereas a realistic distribution with ΔΓ≫⟨Γ⟩ will be studied in a separate paper. 
                    
                    The pair beam can induce two unstable modes: electrostatic and electromagnetic. The growth rate of these instabilities sensitively depends on the momentum spread of the beam. If the momentum spread is small enough, then the instabilities evolve in the so-called reactive regime. In this case, the beam can be mathematically treated as a delta function <cit.> and the growth rates of the electrostatic and electromagnetic instabilities are maximal perpendicular to the direction of the beam propagation <cit.>. As the momentum spread increases, the electromagnetic instability becomes stabilized <cit.>, while the maximum growth rate of the electrostatic mode shifts to the direction parallel to the beam propagation <cit.>. This is the so-called kinetic regime. <cit.> have argued that the momentum spread of the realistic pair beam drastically reduces the growth rate of the electrostatic instability. Later, <cit.> have disputed this statement. <cit.> have demonstrated that the maximum growth rate occurs in the direction almost parallel to the beam (contrary to the reactive regime, when the maximum growth rate occurs in quasi-perpendicular direction to the beam). But <cit.> have assumed the parallel direction of the wave vector from the very beginning. In this case, the electrostatic growth rate, indeed, only weakly depends on the beam temperature <cit.>. Thus, we can conclude that the electrostatic instability for a blazar-induced beam evolves in the kinetic regime at all angles with the maximum growth rate parallel to the beam propagation. It can be shown (see below) that for the beam parameters used by <cit.> and <cit.> the electrostatic instability has evolved in the reactive regime. Thus, an adequate behaviour of the instability has not been simulated before.  
                    
                    So far we have discussed only the electrostatic instability. Usually, the electromagnetic (Weibel) instability can be neglected due to its smaller growth rate, but that is not always the case. If we compare the growth rate of the parallel electrostatic instability with the maximum Weibel growth rate γ_W≈ (V_b/c)(α/Γ)^1/2, then γ_W/γ_react,∥≈α^1/6Γ^1/2=0.3-1 for α=10^-18-10^-15. Thus, the Weibel instability can be potentially competitive with the electrostatic one. Note that we have used γ_W assuming that the beam does not have any momentum spread, and the situation can be different for a beam with a finite temperature. <cit.> have shown that the Weibel instability is strongly suppressed by the non-relativistic perpendicular temperature of the beam. In this work, we will investigate the case of a relativistic temperature and demonstrate that the Weibel instability is suppressed in the case of a realistic blazar-induced beam. Additionally, we will demonstrate that for other conditions (relevant for PIC simulations) this mode can grow. 
                    
                    Consequently, three criteria for a physically relevant simulation setup can be specified: (i) the energy density ratio, ϵ, must be much smaller than unity, (ii) the beam temperature must be high enough, so that the parallel electrostatic instability evolves in the kinetic regime at all angles, and (iii) the Weibel instability must be suppressed. The goal of the current work is to find parameters satisfying all these requirements and to model them using PIC simulations. 
                    
                    In Sec. <ref>, we develop a simple analytical model of plasma instabilities. In Sec. <ref>, we evaluate a condition for the parallel electrostatic instability to be in the kinetic regime. In Sec. <ref>, we discuss our choice of physical parameters for PIC simulations. Sec. <ref> presents simulation results and their discussion. The final summary is given in Sec. <ref>.
                    
                    
                    
                    
                    
                    § ANALYTICAL MODEL
                    
                    
                    We already noted that the electrostatic instability evolves in the kinetic regime and has its maximum growth rate in the direction almost parallel to the beam propagation. At the same time, <cit.> demonstrated that the growth rate of the parallel electrostatic instability very weakly depends on the momentum spread of the beam. Therefore, we can use the well-known growth rate of the two-stream instability for a cold plasma,
                    
                    γ_TS= 3^1/2/2^4/3ω_pα^1/3Γ^-1.
                    
                    
                    Thus, we need to investigate only the electromagnetic Weibel instability. It should be noted that the most unstable wave vector of the Weibel mode can be in transverse direction  to the beam <cit.> as well as in the oblique direction <cit.>. Moreover, the work by <cit.> shows that for dilute beams the maximum growth rates of Weibel mode in the transverse and oblique directions can differ by a factor 2. Therefore, to make a rough estimation, we will study the Weibel instability only for wave vectors perpendicular to the beam. The PIC simulations described in the next section should include oblique modes as well. To derive analytical results, the beam-plasma system is modeled by a waterbag distribution <cit.>. Then, the distributions of the beam and the plasma, respectively, are
                    
                    
                        f_b( p)= n_b 4p_⊥,b^2 (p_∥,b^+-p_∥,b^-)ł[ θł(p_z+p_⊥,b)̊-θł(p_z-p_⊥,b)̊]̊
                        ×ł[ θł(p_y+p_⊥,b)̊-θł(p_y-p_⊥,b)̊]̊ł[ θł(p_x-p_∥,b^-)̊-θł(p_x-p_∥,b^+)̊]̊ ,
                    
                        f_p( p)= n 8p_⊥,p^2 p_∥,pł[ θł(p_z+p_⊥,p)̊-θł(p_z-p_⊥,p)̊]̊
                        ×ł[ θł(p_y+p_⊥,p)̊-θł(p_y-p_⊥,p)̊]̊ł[ θł(p_x+p_∥,p)̊-θł(p_x-p_∥,p)̊]̊ ,
                    
                    where p_∥,b^±=p_0± p_∥,b; p_0 is the beam drift momentum; p_∥,b and p_⊥,b, respectively, the parallel and perpendicular momentum spreads of the beam; p_∥,p and p_⊥,p, respectively, the parallel and perpendicular momentum spreads of the background plasma; θ(x), the Heaviside step function. The beam and background plasma are assumed to be homogeneous with number densities, accordingly, n_b and n. It is useful to consider separately two cases: 
                    (i) p_⊥,b=p_∥,p=p_⊥,p=0 and (ii) p_∥,b=0.
                    
                    
                    
                     §.§ Case p_⊥,b=p_∥,p=p_⊥,p=0
                    
                    
                    We derive the dispersion equation for this case in Appendix A. It reads
                    
                    
                        ł[ 1-ω_p^2ω -ω_b^2ω^2U_1 ]̊ł{  1- ł(kcω)̊^2 - ω_p^2ω^2- .̊
                        ł. -ω_b^2ω^2ł[ ł(kcω)̊^2 U_1   +  ł( 1- ł(kcω)̊^2)̊U_2  ]̊}̊- 
                         - ł(ω_b^2ω^2kcωU_3 )̊^2=0.
                    
                    
                    Taking the limiting case p_∥,b≪ p_0, Eq. (<ref>) provides the classical text book result <cit.>
                        ł[ 1-ω_p^2ω -ω_b^2Γω^2]̊ł{  1- ł(kcω)̊^2 - ω_p^2ω^2-   ω_b^2Γ^3ω^2- .̊
                        ł. - ω_b^2Γω^2ł(kV_oω)̊^2     }̊-  ł(ω_b^2Γω^2kcω)̊^2=0,
                    
                    where V_0=p_0/(m_eΓ). Eq. (<ref>) predicts an instability with growth rate <cit.>:
                    
                    γ=ω_pβ_0 ł(αΓ)̊^1/2kcł(ω_p^2+ (kc)^2 )̊^1/2.
                    
                    
                    Now, we will show that the solution (<ref>) is only slightly different for a large parallel momentum spread p_∥,b≫ p_0. Assuming p_∥,b^+≫ mc and p_∥,b^-≪ mc, we obtain
                    
                    
                    U_1≈ U_3≈m_ec p_∥,b^+lnp_∥,b^+ m_ec,
                    
                    U_2≈m_ec p_∥,b^+.
                    
                    
                    Neglecting unity in each bracket in Eq. (<ref>) results in the solution
                    
                    γ= ω_p ł[ αm_e c p_∥,b^+lnp_∥,b^+ m_ec]̊^1/2kcł(ω_p^2+ (kc)^2 )̊^1/2.
                    
                    
                    It is well seen from Eq. (<ref>) that even for p_∥,b^+=10^2 p_0, the difference between the solutions (<ref>) and (<ref>) is only a factor of 0.4. Thus, we can neglect the parallel momentum dispersion of the beam and use p_∥,b=0.
                    
                    
                    
                     §.§ Case p_∥,b=0
                    
                    
                    The dispersion equation for p_∥,b=0 is derived in Appendix B and has the following form 
                    
                    
                        ł[ 1 - ω_p,p^2ω^2- ω_p,e^2 ω^2- (kv_⊥,p)^2 - ω_b,⊥^2ω^2 - (ku)^2]̊×
                        ×ł[ 1 - (kc)^2+ω_p^2 + ω_b,∥^2 ω^2 -  (kv_∥,p)^2 ω_p,e^2 3ω^2( ω^2 - (kv_⊥,p)^2 )    - .̊
                        ł. 
                        - (kV_0)^2ω_b,⊥^2 ω^2 (ω^2 - (ku)^2)]̊ - ł[ kV_0ω_b,⊥^2 ω(ω^2- (ku)^2 )]̊^2=0,
                    
                    where ω_b,⊥^2=ω_b^2/Γ, ω_b,∥^2=ω_b^2/Γ^3. In principle, one can analyze Eq. (<ref>) analytically, but it is more useful and easier to treat two limiting cases of the cold background plasma and the cold beam.
                    
                    
                    
                      §.§.§ Cold background plasma v_∥,p=v_⊥,p=0
                    
                    
                    Neglecting unity in each bracket in Eq. (<ref>), we obtain that for p_⊥,b≥ p_0(α/Γ)^1/2 the solution is purely real (no instability can arise), whereas for p_⊥,b< p_0(α/Γ)^1/2 the Weibel mode is unstable for k<(ω_p/c)[(α/Γ)(p_0/p_⊥,b)^2-1]^1/2 with growth rate
                    
                    γ= ł( ω_p^2β_0^2αΓ(kc)^2  (kc)^2+ω_p^2  - (ku)^2 )̊^1/2 .
                    
                    
                    Let us now assume that the beam obeys a relativistic Maxwellian distribution: 
                    
                    f( p)= μ 4π (m_ec)^3Γ^2 K_2(μ/Γ)e^-μł[ ł( 1+p^2(m_ec)^2)̊^1/2- β_0 p_x m_ec]̊,
                    
                    where μ=Γμ_R =Γ m_ec^2/(k_BT_b), β_0=V_0/c. Here, T_b is the temperature of the beam in its rest frame. Then we can evaluate Δ p_⊥ (see Appendix C) and write the condition for the Weibel mode stability as
                    α≤α_W= 2-π/2 Γμ_R=(2-π/2) k_BT_b/Γ m_ec^2.
                    
                    
                    In the simulations by <cit.>, magnetic-field fluctuations grew at early times due to the Weibel instability driven by the beam, because α=10^-2 and α_W<5×10^-4 led to condition (<ref>) not being fulfilled. But in the simulations by <cit.>, the Weibel mode was suppressed, because α=2×10^-6<α_W=10^-5. For a realistic blazar-induced beam, the Weibel instability is also suppressed, since k_BT_b≈ m_ec^2 and α≪ 1/<Γ>.
                    
                    
                    
                      §.§.§ Cold beam p_⊥,b=0
                    
                    
                    Again neglecting unity in Eq. (<ref>), we can approximate it as
                    
                    
                    Eω^4 + Fω^2 +G = 0.
                    
                    where
                    
                    
                    E= (kc)^2+ω_p^2 , 
                    
                    F= ł(kV_0ω_b,⊥)̊^2+ 13(kv_∥,pω_p,e)^2  - (kv_⊥,p)^2((kc)^2+ω_p^2) ,
                    
                    G= - ł(kV_0)̊^2(kv_⊥,p)^2ω_b,⊥^2.
                    
                    
                    The growth rate reads
                    
                    γ= ł[F+ (F^2-4EG)^1/2 2E  ]̊^1/2≈ł[ F+|F|2E- G|F|]̊^1/2. 
                    
                    
                    If F>0, Eq. (<ref>) describes the classical Weibel instability <cit.> with growth rate
                    
                    γ≈ł( F E)̊^1/2≈ł( (kv_∥,pω_p,e)^2/3  - (kv_⊥,p)^2((kc)^2+ω_p^2)  (kc)^2+ω_p^2  )̊^1/2.
                    
                    
                    Due to v_⊥,p≠0, the instability is stabilized at large wave vectors, but at small k the plasma is unstable for (v_∥,p/v_⊥,p)^2>3. These conditions were fulfilled in the simulations by <cit.>, where there was a growth of the magnetic-field fluctuations at later time around ω_p,et≈ 10^4. In the opposite case F<0 and assuming that v_⊥,p is large enough, Eq. (<ref>) reduces to Eq. (<ref>). 
                    
                    
                    
                    § CONDITION FOR THE KINETIC REGIME
                    
                    
                    The parallel electrostatic instability evolves in the kinetic regime <cit.>, if
                    
                    ł| v_∥,b c |̊≫α^1/3Γ^-1
                    which can be re-written as
                    α≪α_kin= ł(Γł| v_∥,b/c|̊)̊^3.
                    
                    
                    An analytical expression for α_kin is derived in Appendix C, and its functional behavior is illustrated in Fig. <ref>. 
                    
                    For the simulation parameters used by <cit.>, μ_R=5×10^3 (T_b=10^6 K), Γ=10, and α=2.5×10^-6, we obtain α_kin≈2.8× 10^-9 and Eq. (<ref>) is not fulfilled. 
                    For the work of <cit.>, Γ=300, α=10^-2, μ>3, it results in α_kin≈ 7.1×10^-9, and Eq. (<ref>) is not satisfied again. Hence, both <cit.> and <cit.> did not simulate the electrostatic instability in the appropriate kinematic regime of pair cascades from AGN.
                    
                    
                    
                    
                    § CHOICE OF PARAMETERS FOR PIC SIMULATIONS
                    
                    
                    In the introduction, we have specified three criteria for a physically relevant setup for the beam-plasma system. First, the energy density ratio must satisfy ϵ=αΓ m_ec^2/(k_B T_p)≪1 yielding 
                    α≪α_ϵ(Γ)= k_BT_pΓ m_ec^2≈10^-6-10^-3Γ.
                    
                    The behavior of α_ϵ(Γ) is shown in Figures <ref>, <ref>, and <ref> as red line.
                    Second, the electrostatic instability should develop in the kinetic regime at all angles which is determined by Eq. (<ref>), for which we indicate α_kin by the green line in the figures.
                    Lastly, the Weibel mode must be stable which requires satisfying Eq. (<ref>). Eq. (<ref>) is automatically fulfilled due to α_kin≪α_W for Γ>1 and μ_R>1. Fig. <ref> compares the functions α_kin(Γ) and α_ϵ(Γ). To satisfy Eqs. (<ref>) and (<ref>) for given values of Γ and μ, the value of α must be below both curves α_kin(Γ) and α_ϵ(Γ). We defined a simulation setup, henceforth referred to as run 1, that would satisfy all criteria. The main parameter values are Γ=5 and α=2×10^-4, and it is indicated in Figure <ref> by a black dot. 
                    
                    In addition, we have specified three other setups (runs 2-4) that are listed in 
                    Table <ref>. The goal of these tests is to determine the impact of a violation of one of the criteria on the beam-plasma evolution. For run 2, the energy density ratio ϵ=2.5 is higher than unity, and one might expect a strong heating of the background plasma and subsequently the development of other instabilities. Run 3 considers the evolution of the electrostatic instability in the reactive regime (α>α_kin), and beam energy losses are expected to be larger. Finally, all the conditions are violated for run 4. The values of (α;Γ) for runs 2-4 are demonstrated by the black dots in Figs. <ref>-<ref>, respectively.
                    
                    
                    
                    § THE SIMULATION CODE
                    
                    For the simulation purposes we use EPOCH 2D, a multi-dimensional, fully electromagnetic, relativistic particle-in-cell code developed by the Collaborative Computational Plasma
                    Physics (CCPP) consortium and funded
                    by the Engineering and Physical Sciences Research Council
                    (EPSRC). PIC codes solve Maxwell’s equations on a numerical (Eulerian) grid and follow charged computational particles (CP) as they move under the influence of the electromagnetic field and provide charge and current density <cit.>.
                    
                    The relevant equations are
                    μ_0ϵ_0  ∂ E(x,t)/∂ t =∇× B(x,t) -μ_0 J(x,t)
                    
                    and 
                    ∂ B(x,t)/∂ t =-∇× E(x,t)  ,
                    
                    where the current density, J(x,t), is computed using the algorithm of <cit.>.
                    Collisionless plasma is set up with a Maxwellian velocity distribution. For each CP the field pusher solves the relativistic equation of motion with a numerical approximation of Lorentz force equation. EPOCH is  a refined version of the basic explicit PIC algorithm with higher-order weights and interpolation schemes <cit.>. 
                    Note that the 2D model can break down on the non-linear evolution stage, when 3D mode coupling becomes important <cit.>. As the electrostatic mode involves a narrow resonance, its modeling in a PIC simulation requires a very good wavenumber resolution of the numerical grid <cit.>. This implies a large number of grid points in any direction which we can establish only in 2D. Waves of arbitrary orientation will be included, albeit with only one linear polarization, as is nonlinear wave coupling, provided it does not build on the polarization out of the simulation plane. In the current study we are mainly interested in the linear growth of the instabilities, and so we accept these limitations.
                    
                    
                    The simulation resolves the x–y plane with periodic boundary conditions. The simulation volume is filled with a beam of electrons and positrons and the background plasma of protons and electrons with real mass ratio. We performed a series of tests to verify the stability of the simulation against numerical artifacts. Of particular interest is avoiding artificial plasma heating arising from electric-field noise caused by the charge-density granularity in a particle simulation. We found that using 400 particles per cell and species is required to keep the plasma temperature as desired and the electric-field noise at a level significantly below the intensity of the electrostatic mode. The desired density ratio, α=n_b/n, is established with numerical weights. The simulation box contains 1024 ×1024 cells, each 1/8 of the skin length in size,  λ_e=c/ω_pe=8 Δ_x. The timestep is chosen to satisfy the CFL condition and to resolve the plasma frequency, ω_pe=(n_0e^2/ϵ_0m_e)^1/2.  
                    
                    Table <ref> lists the temperature of the IGM plasma, T_p, and of the beam in its rest frame, T_b. The IGM particles are initially at rest, while the beam is moving in x-direction with Lorentz factor Γ_b. For the IGM, EPOCH code generates a non relativistic thermal distribution using the method of <cit.>. However, we implemented the algorithm of <cit.> to set up the relativistic Maxwellian distribution for the beam. 
                    For the graphical presentation 
                    we use the following normalization: distance and time
                    are normalized to c/ω_p,e and ω_p,e^-1, and electric and magnetic fields are given in units of ω_p,ecm_e/e and ω_p,em_e/e, respectively. 
                    
                    
                    
                    In order to reduce the well-known PIC-code phenomena of self heating and statistical noise, all simulations are performed with a high number of CPs (400 particles per species), a 6th-order field particle pusher, and a triangular-shaped cloud (TSC) shape function, with the peak of the triangle located at the position of the pseudoparticle. 
                    
                    
                    
                    § DISCUSSION OF SIMULATION RESULTS
                    
                    
                     §.§ Run 1
                    
                    
                    As mentioned above, for run 1 all relevant criteria for the beam are fulfilled. First of all, the beam/plasma energy density ratio ϵ=0.5 for run 1 is smaller than unity. Moreover, the beam is stable with respect to the Weibel instability, while the electrostatic mode grows as expected in the kinetic regime, i.e. at the parallel wave vector k_||≈ω_p/c to the beam.
                    
                    In Fig. <ref> we present the Fourier spectrum of the electric field, and it is evident that an electrostatic mode with 𝐄∥𝐤 dominates with peak intensity for wave vectors roughly aligned with the beam direction. The linear growth rate of the electric field is about γ≃ 4×10^-4ω_pe. The theoretically calculated maximum growth rate for parallel wave vectors is 5×10^-4ω_pe which approximately agrees with that derived numerically. 
                    
                    Figs. <ref> and <ref> demonstrate that after 16,237  ω_pe^-1, corresponding to about 8 growth times, the instability has saturated with negligible energy loss and heating of the beam. The latter is of interest because a widening of the lateral beam distribution would impose a temporal smearing of the ICS signal that would reduce the expected flux seen with Fermi-LAT. Our run 1 suggests that this effect is not efficient for realistic pair beams induced by gamma rays from AGN.
                    
                    
                    
                    Fig. <ref> illustrates the time evolution of the electric and magnetic field energy density. The electric field energy saturates after ∼ 7 growth times. It is clear that the beam transferred only a tiny fraction (∼ 10^-4%) of its initial kinetic energy into the electromagnetic fields. Accordingly, the change of the beam distribution is also very small (see Figs. <ref>-<ref>). This development of the beam-plasma interaction is caused by the initial momentum spread of the beam. It was also found by <cit.> that the beam momentum distribution does not relax to the plateau form when Δ p_⊥,b/m_e c∼1. The physical reason is that the electrostatic growth rate simply becomes much smaller than in the reactive regime. At the same time, the damping rates of the modulation instability and non-linear Landau damping depend on the resonant wave energy, and therefore they will stabilize the instability at smaller electric field energies. 
                    
                    
                    
                     §.§ Run 2
                    
                    
                    In contrast to run 1, run 2 considers the beam/plasma energy-density ratio, ϵ=5, greater than 1. The only parameter changed compared to run 1 is the plasma temperature that became by an order of magnitude smaller. Due to the fact that the beam parameters remained the same, the Weibel mode is still stable. The electrostatic instability also evolves in kinetic regime with a growth rate around ≃ 5×10^-4ω_pe, and the time evolution of the Fourier spectrum (shown in Fig. <ref> at ω_pet=5036) is consistent with the value.
                    
                    Fig. <ref> shows that the electric field energy density saturates at nearly the same level as in run 1. Note that due to a smaller plasma temperature the initial electric noise level in run 2 is about by an order of magnitude smaller compared to run 1. Although the peak intensity of the electrostatic modes is now observed at a 10^∘ angle to the beam direction, the distribution function again did not evolve appreciably, in particular not to a plateau distribution, and the beam experienced only a tiny energy loss or widening.
                    
                    
                    
                     §.§ Run 3
                    
                    
                    With run 3, we explore the reactive regime of the electrostatic mode in contrast to runs 1 and 2, where the instability was kinetic. To do this, we have reduced the temperature of the beam and increased its gamma factor. Now, the electrostatic instability grows at an oblique direction (at about 30^∘) to the beam as is evident from the Fourier spectrum shown in Fig. <ref>. The growth rate for oblique propagation and the parameters of run 3 (assuming a cold beam <cit.>) is
                    γ_TS= 3^1/2/2^4/3ω_peł(αΓ)̊^1/3ł(k_||^2 k^2γ^2  + k_⊥^2 k^2)̊^1/3≃ 9×10^-3ω_pe ,
                    
                    where the last equality applies for the parameters of run 3. The numerically determined growth rate is smaller than that by a factor 2-3. This difference may be explained by the fact that run 3 operates not very far from the condition α=α_kinł(Γ)̊ (see Fig. <ref>).
                    
                    The instability growth rate of run 3 is larger by an order of magnitude compared to runs 1 and 2. Therefore, we can expect a more substantial modification of the beam. Although the electric-field energy density remains small as shown in Fig. <ref>, we observe in Fig. <ref> a significant transverse widening of the beam that is not seen in runs 1 and 2. Fig. <ref> indicates that the width of the perpendicular momentum distribution of the beam increased by a factor of 3. 
                     
                    
                    
                    
                     §.§ Run 4
                    
                    
                    Finally, run 4 considers a situation in which all three constraints on the beam parameters are violated. Unlike runs 1 and 2, the fastest electrostatic mode develops for wave vectors that are quasi-perpendicular to the beam, as is well seen in Fig. <ref>. The numerical growth rate perfectly agrees with the analytical estimation for a cold beam, and it is about 2.2×10^-2ω_pe which is larger than in the runs 1, 2, and 3. Furthermore, Fig. <ref> demonstrates that the electric-field energy density assumes a considerably higher value than in three other runs on account of a higher growth rate. At the same time, the Weibel mode is destabilized resulting in a strong growth of magnetic field to a field strength even larger (see the red line in Fig. <ref>) than that of the electric field. Actually, the orange line in Fig. <ref> indicates that the dominant energy transfer is that to IGM electrons (∼ 0.5%), while the magnetic field receives only ∼ 10^-5%.
                    
                    This affects the momentum distribution of both the beam and IGM, as we present in Figs. <ref>-<ref>. In Fig. <ref> we also see a remarkable increase in the momentum spread of the beam and the IGM. This run 4 is similar to the simulations by <cit.> who observed a similar beam-plasma evolution.
                    
                    
                    
                    § SUMMARY
                    
                    
                    We have revisited the issue of plasma instabilities induced by electron-positron beams in the fully ionized intergalactic medium. This problem is related to pair beams produced by TeV radiation of blazars. The main objective of our study is to clarify the feedback of the beam-driven instabilities on the pairs. 
                    
                    The largest difficulty is the impossibility to simulate realistic blazar-induced beams, even with modern computational resources. Therefore, parameters must be found that permit numerical modeling with similar physical properties. Two important criteria of the realistic pair beams have been noticed before: (i), the beam/IGM energy density ratio is much smaller than unity <cit.>, and (ii), the electrostatic mode evolves in the kinetic regime <cit.>. However, the simple estimation presented in the introduction shows that the Weibel mode can potentially compete with the kinetic electrostatic instability. To clarify this point, we have used a simple analytical model and demonstrated that the Weibel mode is actually stable for realistic parameters. This adds a third criterion for the pair beams. 
                    
                    Previous PIC studies of the blazar-induced pair beams <cit.> considered only some of these requirements on the beam-plasma system. In contrast, we have performed a simulation (run 1), for which all of them are taken into account. Then, we have compared this case with three other simulations (runs 2-4), for which some criteria were violated. 
                    The results of run 1 indicate that the pair beam does not experience any significant modification. The electrostatic growth rate turns out to be quite small, and non-linear effects stabilize the beam very efficiently. However, once the electrostatic instability becomes reactive (runs 3-4), as is the case for the studies of <cit.> and <cit.>, the beam momentum distribution widens drastically in the transverse direction. A significant widening of the beam could in principle account for the observed low flux of cascade gamma rays in the GeV band on account of temporal smearing, but that requires a widening by a factor ≫ 10. In any case beam widening is only observed if the instability develops in the reactive regime, and that is not relevant for realistic pair beams arising from interactions of AGN gamma rays with extragalactic background light. Also, if the beam/IGM energy density ratio is high, then the beam effectively heats the IGM (run 4), as was seen in the simulations by <cit.>. 
                    
                    To summarize, we have improved modeling of plasma instabilities for blazar-induced pair beams by including three relevant criteria for the beam. Our results suggest that such instabilities play a negligible role and cannot suppress the flux of cascade gamma rays in the GeV band. Thus, other suppression mechanisms of the energy flux from TeV blazars such magnetic-field deflection must be at play. 
                    
                    
                    
                    The numerical simulations were performed with the EPOCH code that was in part funded by the UK EPSRC grants EP/G054950/1, EP/G056803/1, EP/G055165/1 and EP/ M022463/1. The numerical work was conducted on resources provided by The North-German Supercomputing Alliance (HLRN) under project bbp00003. M.P. acknowledges support through grant PO 1508/1-2 of the Deutsche Forschungsgemeinschaft. The work of J.N. is supported by Narodowe Centrum Nauki through research project DEC-2013/10/E/ST9/00662.
                    
                    
                    § APPENDIX A: DERIVATION OF THE DISPERSION EQUATION FOR P_⊥,B=P_∥,P=P_⊥,P=0 AND K_∥=0
                    
                    
                    In the case p_⊥,b=p_∥,p=p_⊥,p=0, the beam and plasma distributions, respectively, reads
                    
                    
                        f_b( p)= n_b (p_∥,b^+-p_∥,b^-)δ(p_z) δ(p_y) ×
                        ×ł[ θł(p_x-p_∥,b^-)̊-θł(p_x-p_∥,b^+)̊]̊ ,
                    
                    f_p( p)= n   δ(p_x) δ(p_y)δ(p_z) ,
                    
                    where δ(x) is the Dirac delta function. The dielectric tensor is given by <cit.> and <cit.>:
                    
                    
                        ϵ_i,j=δ_i,j+ ∑_a=p,b4π e^2ω^2×
                        ×∫ d^3p ł( v_i  f_a( p) p_j  - v_iv_jk_l  k v-ω f_a( p)  p_l)̊.
                    
                    
                    Evaluating the dielectric tensor (<ref>) for the distribution functions (<ref>)-(<ref>) and for the wave vector k=(0,0,k) yields
                    
                    ϵ_zy=ϵ_yz=ϵ_yx=ϵ_xy=0,
                    ϵ_zz=ϵ_yy=1-ω_p^2ω -ω_b^2ω^2U_1 ,
                    ϵ_xx=1- ω_p^2ω^2-ω_b^2ω^2ł[ ł(kcω)̊^2 U_1   +  ł( 1- ł(kcω)̊^2)̊U_2  ]̊,
                    ϵ_xz=ϵ_zx= - ω_b^2ω^2kcωU_3,
                    
                    
                    where
                    
                    U_1= m_ec p_∥,b^+-p_∥,b^-lnł| p_∥,b^++ [ (p_∥,b^+ )^2 +m_e^2c^2]^1/2 p_∥,b^-+ [ (p_∥,b^-)^2 +m_e^2c^2]^1/2|̊,
                    
                    U_2= m_ec p_∥,b^+-p_∥,b^-ł( p_∥,b^+ [ (p_∥,b^+)^2 +m_e^2c^2]^1/2 - p_∥,b^- [ (p_∥,b^-)^2 +m_e^2c^2]^1/2)̊,
                    
                    U_3=m_ec 2(p_∥,b^+-p_∥,b^-)lnł| (p_∥,b^+)^2 +m_e^2c^2  (p_∥,b^-)^2 +m_e^2c^2 |̊.
                    
                    
                    Here, we have introduced ω_b^2=8π n_be^2/m_e. The dispersion equation reads:
                    
                    
                        (Λ_i,j)=ł(ϵ_i,j+ k_ik_jc^2ω^2- ł(kcω)̊^2δ_i,j)̊= 
                        
                        =Λ_yy(Λ_zzΛ_xx-Λ_zx^2)=0.
                    
                    
                    Thus, the dispersion equation for electromagnetic fluctuations is
                    
                    
                        ł[ 1-ω_p^2ω -ω_b^2ω^2U_1 ]̊ł{  1- ł(kcω)̊^2 - ω_p^2ω^2- .̊
                        ł. -ω_b^2ω^2ł[ ł(kcω)̊^2 U_1   +  ł( 1- ł(kcω)̊^2)̊U_2  ]̊}̊- 
                         - ł(ω_b^2ω^2kcωU_3 )̊^2=0.
                    
                    
                    § APPENDIX B: DERIVATION OF THE DISPERSION EQUATION FOR P_∥,B=0 AND K_∥=0
                    
                    
                    For p_∥,b=0, the distribution function of the beam reads
                    
                    
                        f_b( p)= n_b 4p_⊥,b^2 ł[ θł(p_z+p_⊥,b)̊-θł(p_z-p_⊥,b)̊]̊×
                        ×ł[ θł(p_y+p_⊥,b)̊-θł(p_y-p_⊥,b)̊]̊δ(p_x-p_0).
                    
                    We will assume p_⊥,b≪ p_0. We will still model background protons with the distribution (<ref>), whereas the distribution function of the background electrons is given by Eq. (<ref>). Moreover, we will assume that the background electrons are non-relativistic, p_∥,p=m_e v_∥,p and p_⊥,p=m_e v_⊥,p. Now, it is easy to find that again ϵ_zy=ϵ_yz=ϵ_yx=ϵ_xy=0, but 
                    
                    ϵ_zz= 1 - ω_p,p^2ω^2- ω_p,e^2 ω^2- (kv_⊥,p)^2 - ω_b^2/Γω^2 - (ku)^2,
                    
                        ϵ_yy= 1- ω_p^2+ω_b^2/Γω^2 - (kv_⊥,p)^2 ω_p,e^2 3ω^2( ω^2 - (kv_⊥,p)^2 )    - 
                         
                        -  (ku)^2ω_b^2/Γ 3ω^2 (ω^2 - (ku)^2 ) ,
                    
                        ϵ_xx= 1 - ω_p^2 + ω_b^2/Γ^3 ω^2 -  (kv_∥,p)^2 ω_p,e^2 3ω^2( ω^2 - (kv_⊥,p)^2 )    - 
                         
                        - (kV_0)^2ω_b^2/Γω^2 (ω^2 - (ku)^2),
                    ϵ_xz=ϵ_zx= - kV_0ω_b^2/Γω(ω^2- (ku)^2 ),
                    
                    where ω_p,e=(4π ne^2/m_e)^1/2, ω_p,p=(4π ne^2/m_p)^1/2 (m_p is the proton mass), u=V_0p_⊥,b/p_0.
                    
                    Finally, the Weibel instability is described by the equation
                    
                    
                        ł[ 1 - ω_p,p^2ω^2- ω_p,e^2 ω^2- (kv_⊥,p)^2 - ω_b^2/Γω^2 - (ku)^2]̊×
                        ×ł[ 1 - (kc)^2+ω_p^2 + ω_b^2/Γ^3 ω^2 -  (kv_∥,p)^2 ω_p,e^2 3ω^2( ω^2 - (kv_⊥,p)^2 )    - .̊
                        ł. 
                        - (kV_0)^2ω_b^2/Γω^2 (ω^2 - (ku)^2)]̊ - ł[ kV_0ω_b^2/Γω(ω^2- (ku)^2 )]̊^2=0.
                    
                    
                    § APPENDIX C: APPROXIMATION FOR Α_KIN AT LARGE VALUES OF Μ_R
                    
                    For μ_R≫1, we can use the series expansion 
                    
                    c(p^2+m_e^2c^2)^1/2- V_0p_x ≈mc^2Γ+ (p_x-p_0)^2 2m_eΓ^3 + p_z^2+p_y^2 2m_eΓ .
                    
                    Then Eq. (<ref>) can be approximated as <cit.>
                    f_b( p)= n_bπ^3/2 p_⊥,b^2 p_∥,be^-p_z^2+p_y^2 p_⊥,b^2- (p_x-p_0)^2 p_∥,b^2
                    or
                    
                    f_b( v)= n_bπ^3/2 v_⊥,b^2 v_∥,be^-v_z^2+v_y^2 v_⊥,b^2- (v_x-v_0)^2 v_∥,b^2,
                    
                    where p_⊥,b^2=2 m_ek_BT_b, p_∥,b^2=2Γ^2 m_ek_BT_b, v_⊥,b^2=2k_BT_b/(m_eΓ^2), v_∥,b^2=2k_BT_b/(m_eΓ^4). It is easy to find that
                    
                     p_⊥,b=(2-π/2)^1/2mcμ_R^1/2
                    v_∥, b=ł( ⟨ł(v_x-⟨ v_x⟩)̊^2 ⟩)̊^1/2≈ł( k_B T_b m_eΓ^4)̊^1/2
                    and 
                    α_kin= 1Γ^3 μ_R^3/2.
                    [Arber et al.(2015)]Arber15 Arber, T. D.,  Bennett, K., Brady, C. S. , Lawrence-Douglas, A. ,
                    Ramsay, M. G., Sircombe, N. J. ,  Gillies, P.,  Evans, R. G., Schmitz, H. ,Bell, A. R., and Ridgers, C. P.,  2015, Plasma Phys. Controlled Fusion 57, 113001
                    [Birdsall & Langdon(2004)]Birdsall04Birdsall, C., & Langdon, A. 2004, Plasma  Physics  Via  Computer  Simulation, Series in Plasma Physics (Taylor & Francis, 2004).
                    [Box & Muller(1958)]BoxM58Box, G. E. P., & Muller , M. E. 1958, A
                    note on the generation random normal deviates, Annals of Mathematical Statistics, 29: 610-611
                    [Breizman & Ryutov(1974)]Breizman74Breizman, B. N., & Ryutov, D. D. 1974, Nucl. Fus., 14, 873
                    [Breizman(1990)]Breizman90Breizman, B. N. 1990, RvPP, 15, 61
                    [Bret et al.(2004)]Bret04Bret, A., Firpo, M.-C., & Deutsch, C. 2004, Phys. Rev. E, 70, 046401 
                    [Bret et al.(2005)]Bret05Bret, A., Firpo, M.-C., & Deutsch, C. 2005, Phys. Rev. E, 72, 016403
                    [Bret(2006)]Bret06Bret, A. 2006, EPL (Europhysics Letters), 74, 1027
                    [Bret et al.(2010)]Bret10Bret, A., Gremillet, L., & Dieckmann M. E. 2010, Phys. Plasmas, 17, 120501
                    [Broderick et al.(2012)]Broderick12Broderick, A. E., Chang, P., & Pfrommer, C. 2012, ApJ, 752, 22
                    [Califano et al.(1998)]Califano98Califano, F., Prandi, R., Pegoraro, F., & Bulanov, S. V. 1998, PRE, 58, 7837
                    [Dawson(1983)]Dawson83Dawson, J. M. 1983, Rev. Mod. Phys., 55, 403 
                    [de Naurois(2015)]Nauroisde Naurois, M. 2015, Proceedings of the 34^th International Cosmic Ray Conference, Netherlands
                    [Godfrey et al.(1975)]Godfrey75Godfrey, B. B., Shanahan, W. R., & Thode, L. E. 1975, Phys. Fluids, 18, 346
                    [Kempf et al.(2016)]Kempf16Kempf, A., Kilian, P., & Spanier, F. 2016, A&A, 585, A132
                    [Lazar & Merches(2003)]Lazar03Lazar, M., & Merches, I. 2003, Phys. Lett. A, 313, 418 
                    [Liu et al.(2011)]Liu11Liu, Y., Liu, S. Q., & Li, X. Q. 2011, Contrib. Plasma Phys., 51, 51 
                    [Lominadze & Mikhailovskii(1979)]Lominadze79Lominadze, D. G., & Mikhailovskii, A. B. 1979, Sov. Phys. JETP, 49, 483
                    [Meierovich & Sukhorukov(1976)]Meierovich76Meierovich, B. E., & Sukhorukov, S. T. 1976, Sov. Phys.-JETP, 41, 895
                    [Miniati & Elyiv(2013)]Miniati13Miniati, F., & Elyiv, A. 2013, ApJ, 770,54
                    [Neronov & Semikoz(2009)]Neronov09Neronov, A., & Semikoz, D. V. 2009, Phys. Rev. D, 80, 123012
                    [Neronov & Vovk(2010)]Neronov10Neronov, A., & Vovk, I. 2010, Science, 328, 73
                    [Schlickeiser et al.(2002)]RS02Schlickeiser, R., Vainio, R., Böttcher, M., Lerche, I., Pohl, M., & Schuster, C. 2002, A&A, 393, 69
                    [Schlickeiser(2004)]RS04Schlickeiser, R. 2004, Phys. Plasmas, 11, 5532
                    [Schlickeiser et al.(2012a)]RS12aSchlickeiser, R., Ibscher, D., & Supsar, M. 2012, ApJ, 758, 101
                    [Schlickeiser et al.(2012b)]RS12Schlickeiser, R., Ibscher, D., & Supsar, M. 2012, ApJ, 758, 102
                    [Schlickeiser et al.(2013)]RS13Schlickeiser, R., Krakau, S., & Supsar, M. 2013, ApJ, 777, 49
                    [Shalaby et al.(2017)]2017arXiv170400014S Shalaby, M., Broderick, A. E., Chang, P., et al. 2017, arXiv:1704.00014 
                    [Sironi & Giannios(2014)]Sironi14Sironi, L., & Giannios, D. 2014, ApJ, 787, 49
                    [Taylor et al.(2011)]Taylor11Taylor. A. M., Vovk, I., & Neronov, A. 2011, A&A, 529, A144
                    [Villasenor(1992)]Villasenor92Villasenor, J.,& Buneman, O. 1992, Comput. Phys. Commun. 69, 306
                    [Watson et al.(1960)]Watson60Watson, K. M., Bludman, S. A., and Rosenbluth, M. N. 1960, Phys. Fluids, 3, 741
                    [Yoon & Davidson(1987)]Yoon87Yoon, P. H., & Davidson, R. C. 1987, Phys. Rev. A, 35, 2718 
                    [Zenitani(2015)]Zenitani15Zenitani, S. 2015, Phys. Plasmas, 22, 042116 ] ]
