EDP Sciences
Free Access
Issue
A&A
Volume 607, November 2017
Article Number A112
Number of page(s) 11
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201731127
Published online 23 November 2017

© ESO, 2017

1. Introduction

Many blazars, a sub-class of active galactic nuclei (AGN), 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 (de Naurois 2015). These very-high-energy photons interact with extragalactic background light (EBL) producing ultra-relativistic electron-positron pairs with the typical Lorentz factor 105< Γ < 107 (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 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 (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 multi-TeV gamma-ray 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 (two-stream, 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, wave-particle interactions can reduce the energy of the pairs by 3050% (Schlickeiser et al. 2002; Bret et al. 2010). 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), α = nb/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 non-linear evolution of the plasma system is also problematic. In this work, we combine numerical particle-in-cell (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 (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 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 ⟨ Γ ⟩ = 105, Γ = 103−108, nb = 10-25−10-19 cm-3, Δθ ≈ 1 / ⟨ Γ ⟩ ≈ 10-5 (Δθ is the angular spread), whereas typical parameters of the IGM are T = 104−107 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 ϵ = nb ⟨ Γ ⟩ mec2/ (nkBT) ≈ 10-10−10-1 (kB is the Boltzmann constant, me 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 α = nb/n ≈ 10-2, Γ ≈ 102, and kBT/ (mec2) ≈ 10-8, providing with ϵ ≈ 108, 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 so-called 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 so-called 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 quasi-perpendicular 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 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 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 ≈ (Vb/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 non-relativistic 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 blazar-induced 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 well-known growth rate of the two-stream 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 beam-plasma 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 ; p0 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, nb 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∥ ,bp0, Eq. (4) provides the classical textbook result (Breizman 1990) (5)where V0 = p0/ (meΓ). 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∥ ,bp0. 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⊥ ,bp0(α/ Γ)1 / 2 the solution is purely real (no instability can arise), whereas for p⊥ ,b<p0(α/ Γ)1 / 2 the Weibel mode is unstable for k< (ωp/c) [ (α/ Γ)(p0/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 = Γmec2/ (kBTb), β0 = V0/c. Here, Tb 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), 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 (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 blazar-induced beam, the Weibel instability is also suppressed, since kBTbmec2 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 magnetic-field fluctuations at later time around ωp,et ≈ 104. In the opposite case, F< 0 and assuming that v⊥ ,p is large enough, Eq. (18) reduces to Eq. (6).

thumbnail Fig. 1

Function αkin(Γ): dashed dotted black line μR = 2.5; dashed black line μR = 10; dotted black line μR = 100. Function αϵ(Γ): green line kBTp/ (mec2) = 4 × 10-3; red line kBTp/ (mec2) = 10-4; blue line kBTp/ (mec2) = 10-5; brown line kBTp/ (mec2) = 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 re-written 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 × 103 (Tb = 106 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.

thumbnail 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, kBTp/ (mec2) = 4 × 10-4. The black dot indicates parameter values of run 2.

Open with DEXTER

thumbnail Fig. 3

Same as Fig. 2, but for μR = 10, kBTp/ (mec2) = 4 × 10-3. The black dot indicates parameter values of run 3.

Open with DEXTER

thumbnail Fig. 4

Same as Fig. 2, but for μR = 10, kBTp/ (mec2) = 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 beam-plasma system. First, the energy density ratio must satisfy ϵ = αΓmec2/ (kBTp) ≪ 1 yielding (22)The behavior of αϵ(Γ) is shown in Figs. 24 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 24) 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 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 24 are demonstrated by the black dots in Figs. 24, respectively.

5. 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 (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 higher-order weights and interpolation schemes (Arber et al. 2015). We note that the two-dimensional (2D) model can break down on the non-linear evolution stage, when three-dimensional (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 non-linear 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, α = nb/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 Courant-Friedrichs-Lewy (CFL) condition and to resolve the plasma frequency, ωpe = (n0e2/ϵ0me)1 / 2.

Table 1 lists the temperature of the IGM plasma, Tp, and of the beam in its rest frame, Tb. 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 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,ecme/e and ωp,eme/e, respectively.

Table 1

Simulation parameters.

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.

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 Ek 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 Fermi-LAT. Our run 1 suggests that this effect is not efficient for realistic pair beams induced by gamma rays from AGN.

thumbnail Fig. 5

Two-dimensional Fourier spectrum of Ek at ωpet = 4222 for run 1.

Open with DEXTER

thumbnail Fig. 6

Time evolution of the energy densities of electric and magnetic field, respectively, in SI units for run 1.

Open with DEXTER

thumbnail Fig. 7

Beam momentum distribution in px for run 1 at two points in time.

Open with DEXTER

thumbnail Fig. 8

Beam momentum distribution in py 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 beam-plasma 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/mec ~ 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.

6.2. Run 2

thumbnail Fig. 9

Two-dimensional Fourier spectrum of Ek at ωpet = 5036 for run 2.

Open with DEXTER

thumbnail 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 energy-density 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 ωpet = 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 23. 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 electric-field 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.

thumbnail Fig. 11

Two-dimensional Fourier spectrum of Ek at ωpet = 2448 for run 3.

Open with DEXTER

thumbnail Fig. 12

Time evolution of the energy densities of electric and magnetic field, respectively, in SI units for run 3.

Open with DEXTER

thumbnail Fig. 13

Beam momentum distribution in py for run 3 at two points in time.

Open with DEXTER

thumbnail Fig. 14

Time evolution of the momentum spread of the beam, prms, for run 3.

Open with DEXTER

6.4. Run 4

thumbnail Fig. 15

Two-dimensional Fourier spectrum of Ek at ωpet = 401 for run 4.

Open with DEXTER

thumbnail 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

thumbnail Fig. 17

IGM momentum distribution in px for run 4 at two points in time.

Open with DEXTER

thumbnail Fig. 18

Beam momentum distribution in py for run 4 at two points in time.

Open with DEXTER

thumbnail Fig. 19

Time evolution of the momentum spread of the beam, prms, 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 quasi-perpendicular 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 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. 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 beam-plasma evolution.

7. 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 (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 blazar-induced pair beams (Sironi & Giannios 2014; Kempf et al. 2016) 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 24), 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 non-linear effects stabilize the beam very efficiently. However, once the electrostatic instability becomes reactive (runs 34), 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 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 as magnetic-field 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 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.

References

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⊥ ,bp0. 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 non-relativistic, p∥ ,p = mev∥ ,p and p⊥ ,p = mev⊥ ,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πne2/me)1 / 2, ωp,p = (4πne2/mp)1 / 2 (mp is the proton mass), u = V0p⊥ ,b/p0.

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

(C.4)

(C.5)and (C.6)

All Tables

Table 1

Simulation parameters.

All Figures

thumbnail Fig. 1

Function αkin(Γ): dashed dotted black line μR = 2.5; dashed black line μR = 10; dotted black line μR = 100. Function αϵ(Γ): green line kBTp/ (mec2) = 4 × 10-3; red line kBTp/ (mec2) = 10-4; blue line kBTp/ (mec2) = 10-5; brown line kBTp/ (mec2) = 10-6. The black point illustrates the parameters chosen for simulation run 1 that satisfies all criteria.

Open with DEXTER
In the text
thumbnail 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, kBTp/ (mec2) = 4 × 10-4. The black dot indicates parameter values of run 2.

Open with DEXTER
In the text
thumbnail Fig. 3

Same as Fig. 2, but for μR = 10, kBTp/ (mec2) = 4 × 10-3. The black dot indicates parameter values of run 3.

Open with DEXTER
In the text
thumbnail Fig. 4

Same as Fig. 2, but for μR = 10, kBTp/ (mec2) = 10-3. The black dot indicates parameter values of run 4.

Open with DEXTER
In the text
thumbnail Fig. 5

Two-dimensional Fourier spectrum of Ek at ωpet = 4222 for run 1.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 7

Beam momentum distribution in px for run 1 at two points in time.

Open with DEXTER
In the text
thumbnail Fig. 8

Beam momentum distribution in py for run 1 at two points in time.

Open with DEXTER
In the text
thumbnail Fig. 9

Two-dimensional Fourier spectrum of Ek at ωpet = 5036 for run 2.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 11

Two-dimensional Fourier spectrum of Ek at ωpet = 2448 for run 3.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 13

Beam momentum distribution in py for run 3 at two points in time.

Open with DEXTER
In the text
thumbnail Fig. 14

Time evolution of the momentum spread of the beam, prms, for run 3.

Open with DEXTER
In the text
thumbnail Fig. 15

Two-dimensional Fourier spectrum of Ek at ωpet = 401 for run 4.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 17

IGM momentum distribution in px for run 4 at two points in time.

Open with DEXTER
In the text
thumbnail Fig. 18

Beam momentum distribution in py for run 4 at two points in time.

Open with DEXTER
In the text
thumbnail Fig. 19

Time evolution of the momentum spread of the beam, prms, for run 4.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.