Plasma effects on relativistic pair beams from TeV blazars: PIC simulations and analytical predictions
^{1} Institute for Physics and Astronomy, University of Potsdam, 14476 Potsdam, Germany
email: rafighi@unipotsdam.de
^{2} DESY, Platanenallee 6, 15738 Zeuthen, Germany
^{3} Instytut Fizyki Ja¸drowej PAN, ul. Radzikowskiego 152, 31342 Kraków, Poland
Received: 8 May 2017
Accepted: 22 August 2017
Pair beams produced by veryhighenergy radiation from TeV blazars emit gamma rays in the GeV band by inverseCompton scattering of soft photons. The observed GeVband 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 beamplasma instabilities that we consider in the present work. For realistic parameters the pair beams cannot 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 twodimensional (2D) particleincell (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.
Key words: galaxies: active / instabilities / waves / relativistic processes / plasmas / gamma rays: general
© ESO, 2017
1. Introduction
Many blazars, a subclass of active galactic nuclei (AGN), have been detected with gammaray telescopes such as HESS, VERITAS, Fermi and MAGIC as sources of gammarays with the energy E ≥ 100 GeV (de Naurois 2015). These veryhighenergy photons interact with extragalactic background light (EBL) producing ultrarelativistic electronpositron pairs with the typical Lorentz factor 10^{5}< Γ < 10^{7} (Schlickeiser et al. 2012a; Miniati & Elyiv 2013). These created pairs are the 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 gammaray emission with characteristic energy in the GeV band, but, as indicated by FermiLAT data, the GeV gammaray 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 (Neronov & Semikoz 2009; Neronov & Vovk 2010; Taylor et al. 2011) which has led to its constrainment. However, these constraints on the IGMF are valid only under the assumptions that the multiTeV gammaray emission persists on long timescales and that the pairs lose their energy only due to ICS. This latter assumption is crucial and remains debatable.
The importance of collective plasma effects has been pointed out by several authors (Broderick et al. 2012; Schlickeiser et al. 2012b; Miniati & Elyiv 2013). In fact, the pairs can induce electrostatic (twostream, oblique) and electromagnetic (filamentation, Weibel) instabilities (Breizman & Ryutov 1974; Breizman 1990; Bret et al. 2004, 2005, 2010; Bret 2006; Godfrey et al. 1975; Lominadze & Mikhailovskii 1979). In this case, waveparticle interactions can reduce the energy of the pairs by 30−50% (Schlickeiser et al. 2002; Bret et al. 2010). Therefore, the collective plasma effects can also substantially suppress the GeVband gammaray 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 (Sironi & Giannios 2014; Kempf et al. 2016). At the same time, an accurate analytical description of the nonlinear evolution of the plasma system is also problematic. In this work, we combine numerical particleincell (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) nonlinear damping of plasma waves. Here, we are concerned only with the first two subjects, and the nonlinear effects (Lazar & Merches 2003; Liu et al. 2011) will be analyzed in future papers. Such treatment is possible in the linear stage that 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 gammaray 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} (Broderick et al. 2012; Schlickeiser et al. 2012b; Sironi & Giannios 2014; Miniati & Elyiv 2013). Thus, the energy density ratio is ϵ = n_{b} ⟨ Γ ⟩ m_{e}c^{2}/ (nk_{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 Kempf et al. (2016) who conducted simulations for ϵ = 0.1. The parameters of the simulations by Sironi & Giannios (2014) are α = n_{b}/n ≈ 10^{2}, Γ ≈ 10^{2}, and k_{B}T/ (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 is shown below. We note that Kempf and Sironi have studied a beam distribution with ΔΓ ≪ ⟨ Γ ⟩. We 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 socalled reactive regime. In this case, the beam can be mathematically treated as a delta function (Schlickeiser et al. 2012b) and the growth rates of the electrostatic and electromagnetic instabilities are maximal perpendicular to the direction of the beam propagation (Godfrey et al. 1975). As the momentum spread increases, the electromagnetic instability becomes stabilized (Bret et al. 2005), while the maximum growth rate of the electrostatic mode shifts to the direction parallel to the beam propagation (Breizman 1990). This is the socalled kinetic regime. Miniati & Elyiv (2013) have argued that the momentum spread of the realistic pair beam drastically reduces the growth rate of the electrostatic instability. Later, Schlickeiser et al. (2013) disputed this statement. Sironi & Giannios (2014) 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 quasiperpendicular direction to the beam). But Schlickeiser et al. (2013) 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 (Bret et al. 2005). Thus, we can conclude that the electrostatic instability for a blazarinduced 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 Kempf et al. (2016) and Sironi & Giannios (2014) the electrostatic instability has evolved in the reactive regime. An adequate behavior of the instability has not yet been simulated.
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 potentially be competitive with the electrostatic one. We 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. Bret et al. (2005) have shown that the Weibel instability is strongly suppressed by the nonrelativistic perpendicular temperature of the beam. In this work, we investigate the case of a relativistic temperature and demonstrate that the Weibel instability is suppressed in the case of a realistic blazarinduced beam. Additionally, we 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 Sect. 2, we develop a simple analytical model of plasma instabilities. In Sect. 3, we evaluate a condition for the parallel electrostatic instability to be in the kinetic regime. In Sect. 4, we discuss our choice of physical parameters for PIC simulations. Section 6 presents simulation results and their discussion. The final summary is given in Sect. 7.
2. 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, Schlickeiser et al. (2013) 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 wellknown growth rate of the twostream instability for a cold plasma, (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 (Califano et al. 1998) as well as in the oblique direction (Bret et al. 2010). Moreover, the work by Bret et al. (2010) 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 following section should include oblique modes as well. To derive analytical results, the beamplasma system is modeled by a waterbag distribution (Bret et al. 2005; Yoon & Davidson 1987). Then, the distributions of the beam and the plasma, respectively, are (2)(3)where ; 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; and θ(x) is 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 separately consider two cases: (i) p_{⊥ ,b} = p_{∥ ,p} = p_{⊥ ,p} = 0 and (ii) p_{∥ ,b} = 0.
2.1. Case p_{⊥,b} = p_{∥,p} = p_{⊥,p} = 0
We derive the dispersion equation for this case in Appendix A. It reads (4)Taking the limiting case p_{∥ ,b} ≪ p_{0}, Eq. (4) provides the classical textbook result (Breizman 1990) (5)where V_{0} = p_{0}/ (m_{e}Γ). Equation (5) predicts an instability with growth rate (Godfrey et al. 1975): (6)Now, we will show that the solution (6) is only slightly different for a large parallel momentum spread p_{∥ ,b} ≫ p_{0}. Assuming and , we obtain Neglecting unity in each bracket in Eq. (4) results in the solution (9)It is easily seen from Eq. (9) that even for , the difference between the solutions (8) and (9) is only a factor of 0.4. Thus, we can neglect the parallel momentum dispersion of the beam and use p_{∥ ,b} = 0.
2.2. Case p_{∥,b} = 0
The dispersion equation for p_{∥ ,b} = 0 is derived in Appendix B and has the following form (10)where , . In principle, one can analyze Eq. (10) analytically, but it is more useful and easier to treat two limiting cases of the cold background plasma and the cold beam.
2.2.1. Cold background plasma v_{∥,p} = v_{⊥,p} = 0
Neglecting unity in each bracket in Eq. (10), 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 (11)Let us now assume that the beam obeys a relativistic Maxwellian distribution: (12)where μ = Γμ_{R} = Γm_{e}c^{2}/ (k_{B}T_{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 (13)In the simulations by Sironi & Giannios (2014), magneticfield 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 (13) not being fulfilled. But in the simulations by Kempf et al. (2016), the Weibel mode was suppressed, because α = 2 × 10^{6}<α_{W} = 10^{5}. For a realistic blazarinduced beam, the Weibel instability is also suppressed, since k_{B}T_{b} ≈ m_{e}c^{2} and α ≪ 1 / ⟨ Γ ⟩.
2.2.2. Cold beam p_{⊥,b} = 0
Again neglecting unity in Eq. (10), we can approximate it as (14)where (15)(16)(17)The growth rate reads (18)If F> 0, Eq. (18) describes the classical Weibel instability (Bret et al. 2005) with growth rate (19)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 Sironi & Giannios (2014), where there was a growth of the magneticfield fluctuations at later time around ω_{p,e}t ≈ 10^{4}. In the opposite case, F< 0 and assuming that v_{⊥ ,p} is large enough, Eq. (18) reduces to Eq. (6).
Fig. 1 Function α_{kin}(Γ): dashed dotted black line μ_{R} = 2.5; dashed black line μ_{R} = 10; dotted black line μ_{R} = 100. Function α_{ϵ}(Γ): green line k_{B}T_{p}/ (m_{e}c^{2}) = 4 × 10^{3}; red line k_{B}T_{p}/ (m_{e}c^{2}) = 10^{4}; blue line k_{B}T_{p}/ (m_{e}c^{2}) = 10^{5}; brown line k_{B}T_{p}/ (m_{e}c^{2}) = 10^{6}. The black point illustrates the parameters chosen for simulation run 1 that satisfies all criteria. 

Open with DEXTER 
3. Condition for the kinetic regime
The parallel electrostatic instability evolves in the kinetic regime (Breizman 1990), if (20)which can be rewritten as (21)An analytical expression for α_{kin} is derived in Appendix C, and its functional behavior is illustrated in Fig. 1.
For the simulation parameters used by Kempf et al. (2016), μ_{R} = 5 × 10^{3} (T_{b} = 10^{6} K), Γ = 10, and α = 2.5 × 10^{6}, we obtain α_{kin} ≈ 2.8 × 10^{9} and Eq. (21) is not fulfilled. For the work of Sironi & Giannios (2014), Γ = 300, α = 10^{2}, μ> 3, it results in α_{kin} ≈ 7.1 × 10^{9}, and Eq. (21) is not satisfied again. Hence, both Kempf et al. (2016) and Sironi & Giannios (2014) did not simulate the electrostatic instability in the appropriate kinematic regime of pair cascades from AGN.
Fig. 2 Dependence of all three constraints on the beam Lorentz factor. The green line represents α_{kin}, the red line α_{ϵ}, and the blue line α_{W}. μ_{R} = 2.5, k_{B}T_{p}/ (m_{e}c^{2}) = 4 × 10^{4}. The black dot indicates parameter values of run 2. 

Open with DEXTER 
Fig. 3 Same as Fig. 2, but for μ_{R} = 10, k_{B}T_{p}/ (m_{e}c^{2}) = 4 × 10^{3}. The black dot indicates parameter values of run 3. 

Open with DEXTER 
Fig. 4 Same as Fig. 2, but for μ_{R} = 10, k_{B}T_{p}/ (m_{e}c^{2}) = 10^{3}. The black dot indicates parameter values of run 4. 

Open with DEXTER 
4. Choice of parameters for PIC simulations
In the introduction, we specified three criteria for a physically relevant setup for the beamplasma system. First, the energy density ratio must satisfy ϵ = αΓm_{e}c^{2}/ (k_{B}T_{p}) ≪ 1 yielding (22)The behavior of α_{ϵ}(Γ) is shown in Figs. 2–4 as a red line. Second, the electrostatic instability should develop in the kinetic regime at all angles which is determined by Eq. (21), for which we indicate α_{kin} by the green line in the figures. Lastly, the Weibel mode must be stable which requires satisfying Eq. (13). Equation (13) is automatically fulfilled due to α_{kin} ≪ α_{W} for Γ > 1 and μ_{R}> 1. Figure 1 compares the functions α_{kin}(Γ) and α_{ϵ}(Γ). To satisfy Eqs. (21) and (22) 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 Fig. 1 by a black dot.
In addition, we have specified three other setups (runs 2−4) that are listed in Table 1. The goal of these tests is to determine the impact of a violation of one of the criteria on the beamplasma 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. 2–4, respectively.
5. The simulation code
For the simulation purposes we use EPOCH 2D, a multidimensional, fully electromagnetic, relativistic particleincell 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 (Dawson 1983; Birdsall & Langdon 2004).
The relevant equations are (23)and (24)where the current density, J(x,t), is computed using the algorithm of Villasenor (1992). 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 higherorder weights and interpolation schemes (Arber et al. 2015). We note that the twodimensional (2D) model can break down on the nonlinear evolution stage, when threedimensional (3D) mode coupling becomes important (Lazar & Merches 2003; Liu et al. 2011). As the electrostatic mode involves a narrow resonance, its modeling in a PIC simulation requires a very good wavenumber resolution of the numerical grid (Shalaby et al. 2017). 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 xy 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 electricfield noise caused by the chargedensity 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 electricfield 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 one eighth of the skin length in size, . The timestep is chosen to satisfy the CourantFriedrichsLewy (CFL) condition and to resolve the plasma frequency, ω_{pe} = (n_{0}e^{2}/ϵ_{0}m_{e})^{1 / 2}.
Table 1 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 xdirection with Lorentz factor Γ_{b}. For the IGM, EPOCH code generates a nonrelativistic thermal distribution using the method of Box & Muller (1958). However, we implemented the algorithm of Zenitani (2015) 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 and , and electric and magnetic fields are given in units of ω_{p,e}cm_{e}/e and ω_{p,e}m_{e}/e, respectively.
Simulation parameters.
In order to reduce the wellknown PICcode phenomena of self heating and statistical noise, all simulations are performed with a high number of CPs (400 particles per species), a 6thorder field particle pusher, and a triangularshaped cloud (TSC) shape function, with the peak of the triangle located at the position of the pseudoparticle.
6. Discussion of simulation results
6.1. 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, that is, at the parallel wave vector k_{} ≈ ω_{p}/c to the beam.
In Fig. 5 we present the Fourier spectrum of the electric field, and it is evident that an electrostatic mode with E ∥ k 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.
Figures 7 and 8 demonstrate that after , corresponding to about eight 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 FermiLAT. Our run 1 suggests that this effect is not efficient for realistic pair beams induced by gamma rays from AGN.
Fig. 5 Twodimensional Fourier spectrum of E ∥ k at ω_{pe}t = 4222 for run 1. 

Open with DEXTER 
Fig. 6 Time evolution of the energy densities of electric and magnetic field, respectively, in SI units for run 1. 

Open with DEXTER 
Fig. 7 Beam momentum distribution in p_{x} for run 1 at two points in time. 

Open with DEXTER 
Fig. 8 Beam momentum distribution in p_{y} for run 1 at two points in time. 

Open with DEXTER 
Figure 6 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. 7 and 8). This development of the beamplasma interaction is caused by the initial momentum spread of the beam. It was also found by Sironi & Giannios (2014) 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 nonlinear Landau damping depend on the resonant wave energy, and therefore they will stabilize the instability at smaller electric field energies.
6.2. Run 2
Fig. 9 Twodimensional Fourier spectrum of E ∥ k at ω_{pe}t = 5036 for run 2. 

Open with DEXTER 
Fig. 10 Time evolution of the energy densities of electric and magnetic field, respectively, in SI units for run 2. 

Open with DEXTER 
In contrast to run 1, run 2 considers the beam/plasma energydensity ratio, ϵ = 5, greater than 1. The only parameter changed compared to run 1 is the plasma temperature that became smaller by an order of magnitude. 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 of around ≃ 5 × 10^{4}ω_{pe}, and the time evolution of the Fourier spectrum (shown in Fig. 9 at ω_{pe}t = 5036) is consistent with the value.
Figure 10 shows that the electric field energy density saturates at nearly the same level as in run 1. We note that due to a smaller plasma temperature the initial electric noise level in run 2 is smaller by about an order of magnitude 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.
6.3. 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 in an oblique direction (at about 30°) to the beam as is evident from the Fourier spectrum shown in Fig. 11. The growth rate for oblique propagation and the parameters of run 3 (assuming a cold beam Breizman 1990) is (25)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 (see Fig. 3).
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 electricfield energy density remains small as shown in Fig. 12, we observe in Fig. 13 a significant transverse widening of the beam that is not seen in runs 1 and 2. Figure 14 indicates that the width of the perpendicular momentum distribution of the beam increased by a factor of 3.
Fig. 11 Twodimensional Fourier spectrum of E ∥ k at ω_{pe}t = 2448 for run 3. 

Open with DEXTER 
Fig. 12 Time evolution of the energy densities of electric and magnetic field, respectively, in SI units for run 3. 

Open with DEXTER 
Fig. 13 Beam momentum distribution in p_{y} for run 3 at two points in time. 

Open with DEXTER 
Fig. 14 Time evolution of the momentum spread of the beam, p_{rms}, for run 3. 

Open with DEXTER 
6.4. Run 4
Fig. 15 Twodimensional Fourier spectrum of E ∥ k at ω_{pe}t = 401 for run 4. 

Open with DEXTER 
Fig. 16 Time evolution of the energy densities of electric, magnetic field, and kinetic energy of IGM, respectively, in SI units for run 4. 

Open with DEXTER 
Fig. 17 IGM momentum distribution in p_{x} for run 4 at two points in time. 

Open with DEXTER 
Fig. 18 Beam momentum distribution in p_{y} for run 4 at two points in time. 

Open with DEXTER 
Fig. 19 Time evolution of the momentum spread of the beam, p_{rms}, for run 4. 

Open with DEXTER 
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 quasiperpendicular to the beam, as is easily seen in Fig. 15. The numerical growth rate perfectly agrees with the analytical estimation for a cold beam, and is about 2.2 × 10^{2}ω_{pe}, which is larger than in runs 1, 2, and 3. Furthermore, Fig. 16 demonstrates that the electricfield 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. 16) than that of the electric field. Actually, the orange line in Fig. 16 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 show in Figs. 17 and 18. In Fig. 19 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 Sironi & Giannios (2014) who observed a similar beamplasma evolution.
7. Summary
We have revisited the issue of plasma instabilities induced by electronpositron 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 beamdriven instabilities on the pairs.
The largest difficulty is the impossibility to simulate realistic blazarinduced 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 (Kempf et al. 2016); and (ii), the electrostatic mode evolves in the kinetic regime (Miniati & Elyiv 2013). 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 blazarinduced pair beams (Sironi & Giannios 2014; Kempf et al. 2016) considered only some of these requirements on the beamplasma 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 relatively small, and nonlinear effects stabilize the beam very efficiently. However, once the electrostatic instability becomes reactive (runs 3−4), as is the case for the studies of Kempf et al. (2016) and Sironi & Giannios (2014), 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 EBL. 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 Sironi & Giannios (2014).
To summarize, we have improved modeling of plasma instabilities for blazarinduced 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 as magneticfield deflection must be at play.
Acknowledgments
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 NorthGerman Supercomputing Alliance (HLRN) under project bbp00003. M.P. acknowledges support through grant PO 1508/12 of the Deutsche Forschungsgemeinschaft. The work of J.N. is supported by Narodowe Centrum Nauki through research project DEC2013/10/E/ST9/00662.
References
 Arber, T. D., Bennett, K., Brady, C. S., et al. 2015, Plasma Phys. Control. Fusion, 57, 113001 [NASA ADS] [CrossRef] [Google Scholar]
 Birdsall, C., & Langdon, A. 2004, Plasma Physics Via Computer Simulation, Series in Plasma Physics (Taylor & Francis) [Google Scholar]
 Box, G. E. P., & Muller, M. E. 1958, Annu. Math. Stat., 29, 610 [CrossRef] [Google Scholar]
 Breizman, B. N. 1990, RvPP, 15, 61 [Google Scholar]
 Breizman, B. N., & Ryutov, D. D. 1974, Nucl. Fusion, 14, 873 [NASA ADS] [CrossRef] [Google Scholar]
 Bret, A. 2006, Europhys. Lett., 74, 1027 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bret, A., Firpo, M.C., & Deutsch, C. 2004, Phys. Rev. E, 70, 046401 [NASA ADS] [CrossRef] [Google Scholar]
 Bret, A., Firpo, M.C., & Deutsch, C. 2005, Phys. Rev. E, 72, 016403 [NASA ADS] [CrossRef] [Google Scholar]
 Bret, A., Gremillet, L., & Dieckmann, M. E. 2010, Phys. Plasmas, 17, 120501 [NASA ADS] [CrossRef] [Google Scholar]
 Broderick, A. E., Chang, P., & Pfrommer, C. 2012, ApJ, 752, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Califano, F., Prandi, R., Pegoraro, F., & Bulanov, S. V. 1998, Phys. Rev. E, 58, 7837 [NASA ADS] [CrossRef] [Google Scholar]
 Dawson, J. M. 1983, Rev. Mod. Phys., 55, 403 [NASA ADS] [CrossRef] [Google Scholar]
 de Naurois, M. 2015, Proc. 34th International Cosmic Ray Conference, Netherlands [Google Scholar]
 Godfrey, B. B., Shanahan, W. R., & Thode, L. E. 1975, Phys. Fluids, 18, 346 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Kempf, A., Kilian, P., & Spanier, F. 2016, A&A, 585, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lazar, M., & Merches, I. 2003, Phys. Lett. A, 313, 418 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, Y., Liu, S. Q., & Li, X. Q. 2011, Contrib. Plasm. Phys., 51, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Lominadze, D. G., & Mikhailovskii, A. B. 1979, Sov. Phys. JETP, 49, 483 [NASA ADS] [Google Scholar]
 Meierovich, B. E., & Sukhorukov, S. T. 1976, Sov. Phys. JETP, 41, 895 [NASA ADS] [Google Scholar]
 Miniati, F., & Elyiv, A. 2013, ApJ, 770, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Neronov, A., & Semikoz, D. V. 2009, Phys. Rev. D, 80, 123012 [NASA ADS] [CrossRef] [Google Scholar]
 Neronov, A., & Vovk, I. 2010, Science, 328, 73 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Schlickeiser, R. 2004, Phys. Plasmas, 11, 5532 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R., Vainio, R., Böttcher, M., et al. 2002, A&A, 393, 69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schlickeiser, R., Ibscher, D., & Supsar, M. 2012a, ApJ, 758, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R., Ibscher, D., & Supsar, M. 2012b, ApJ, 758, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R., Krakau, S., & Supsar, M. 2013, ApJ, 777, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Shalaby, M., Broderick, A. E., Chang, P., et al. 2017, ApJ, 818, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Sironi, L., & Giannios, D. 2014, ApJ, 787, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, A. M., Vovk, I., & Neronov, A. 2011, A&A, 529, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Villasenor, J., & Buneman, O. 1992, Comput. Phys. Commun., 69, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Watson, K. M., Bludman, S. A., & Rosenbluth, M. N. 1960, Phys. Fluids, 3, 741 [NASA ADS] [CrossRef] [Google Scholar]
 Yoon, P. H., & Davidson, R. C. 1987, Phys. Rev. A, 35, 2718 [NASA ADS] [CrossRef] [Google Scholar]
 Zenitani, S. 2015, Phys. Plasmas, 22, 042116 [NASA ADS] [CrossRef] [Google Scholar]
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 (A.1)(A.2)where δ(x) is the Dirac delta function. The dielectric tensor is given by Breizman (1990) and Schlickeiser (2004): (A.3)Evaluating the dielectric tensor (A.3) for the distribution functions (A.1), (A.2) and for the wave vector k = (0,0,k) yields where Here, we have introduced . The dispersion equation reads: (A.11)Thus, the dispersion equation for electromagnetic fluctuations is (A.12)
Appendix B: Derivation of the dispersion equation forp_{∥,b} = 0 and k_{∥} = 0
For p_{∥ ,b} = 0, the distribution function of the beam reads (B.1)We will assume p_{⊥ ,b} ≪ p_{0}. We will still model background protons with the distribution (A.2), whereas the distribution function of the background electrons is given by Eq. (3). Moreover, we will assume that the background electrons are nonrelativistic, 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 (B.2)(B.3)(B.4)(B.5)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_{0}p_{⊥ ,b}/p_{0}.
Finally, the Weibel instability is described by the equation (B.6)
Appendix C: Approximation for α_{kin} at large values of μ_{R}
For μ_{R} ≫ 1, we can use the series expansion (C.1)Then Eq. (12) can be approximated as Watson et al. (1960), Meierovich & Sukhorukov (1976)(C.2)or (C.3)where , , , . It is easy to find that
All Tables
All Figures
Fig. 1 Function α_{kin}(Γ): dashed dotted black line μ_{R} = 2.5; dashed black line μ_{R} = 10; dotted black line μ_{R} = 100. Function α_{ϵ}(Γ): green line k_{B}T_{p}/ (m_{e}c^{2}) = 4 × 10^{3}; red line k_{B}T_{p}/ (m_{e}c^{2}) = 10^{4}; blue line k_{B}T_{p}/ (m_{e}c^{2}) = 10^{5}; brown line k_{B}T_{p}/ (m_{e}c^{2}) = 10^{6}. The black point illustrates the parameters chosen for simulation run 1 that satisfies all criteria. 

Open with DEXTER  
In the text 
Fig. 2 Dependence of all three constraints on the beam Lorentz factor. The green line represents α_{kin}, the red line α_{ϵ}, and the blue line α_{W}. μ_{R} = 2.5, k_{B}T_{p}/ (m_{e}c^{2}) = 4 × 10^{4}. The black dot indicates parameter values of run 2. 

Open with DEXTER  
In the text 
Fig. 3 Same as Fig. 2, but for μ_{R} = 10, k_{B}T_{p}/ (m_{e}c^{2}) = 4 × 10^{3}. The black dot indicates parameter values of run 3. 

Open with DEXTER  
In the text 
Fig. 4 Same as Fig. 2, but for μ_{R} = 10, k_{B}T_{p}/ (m_{e}c^{2}) = 10^{3}. The black dot indicates parameter values of run 4. 

Open with DEXTER  
In the text 
Fig. 5 Twodimensional Fourier spectrum of E ∥ k at ω_{pe}t = 4222 for run 1. 

Open with DEXTER  
In the text 
Fig. 6 Time evolution of the energy densities of electric and magnetic field, respectively, in SI units for run 1. 

Open with DEXTER  
In the text 
Fig. 7 Beam momentum distribution in p_{x} for run 1 at two points in time. 

Open with DEXTER  
In the text 
Fig. 8 Beam momentum distribution in p_{y} for run 1 at two points in time. 

Open with DEXTER  
In the text 
Fig. 9 Twodimensional Fourier spectrum of E ∥ k at ω_{pe}t = 5036 for run 2. 

Open with DEXTER  
In the text 
Fig. 10 Time evolution of the energy densities of electric and magnetic field, respectively, in SI units for run 2. 

Open with DEXTER  
In the text 
Fig. 11 Twodimensional Fourier spectrum of E ∥ k at ω_{pe}t = 2448 for run 3. 

Open with DEXTER  
In the text 
Fig. 12 Time evolution of the energy densities of electric and magnetic field, respectively, in SI units for run 3. 

Open with DEXTER  
In the text 
Fig. 13 Beam momentum distribution in p_{y} for run 3 at two points in time. 

Open with DEXTER  
In the text 
Fig. 14 Time evolution of the momentum spread of the beam, p_{rms}, for run 3. 

Open with DEXTER  
In the text 
Fig. 15 Twodimensional Fourier spectrum of E ∥ k at ω_{pe}t = 401 for run 4. 

Open with DEXTER  
In the text 
Fig. 16 Time evolution of the energy densities of electric, magnetic field, and kinetic energy of IGM, respectively, in SI units for run 4. 

Open with DEXTER  
In the text 
Fig. 17 IGM momentum distribution in p_{x} for run 4 at two points in time. 

Open with DEXTER  
In the text 
Fig. 18 Beam momentum distribution in p_{y} for run 4 at two points in time. 

Open with DEXTER  
In the text 
Fig. 19 Time evolution of the momentum spread of the beam, p_{rms}, for run 4. 

Open with DEXTER  
In the text 