Relativistic Pair Beams from TeV Blazars: A Source of Reprocessed GeV Emission rather than IGM Heating
Abstract
The interaction of TeV photons from blazars with the extragalactic background light produces a relativistic beam of electronpositron pairs streaming through the intergalactic medium (IGM). The fate of the beam energy is uncertain. By means of two and threedimensional particleincell simulations, we study the nonlinear evolution of dilute ultrarelativistic pair beams propagating through the IGM. We explore a wide range of beam Lorentz factors and beamtoplasma density ratios , so that our results can be extrapolated to the extreme parameters of blazarinduced beams ( and , for the most powerful blazars). For cold beams, we show that the oblique instability governs the early stages of evolution, but its exponential growth terminates – due to selfheating of the beam in the transverse direction – when only a negligible fraction of the beam energy has been transferred to the IGM plasma. Further relaxation of the beam proceeds through quasilongitudinal modes, until the momentum dispersion in the direction of propagation saturates at . This corresponds to a fraction of the beam energy being ultimately transferred to the IGM plasma, irrespective of or . If the initial dispersion in beam momentum satisfies (as typically expected for blazarinduced beams), the fraction of beam energy deposited into the IGM is much smaller than . It follows that at least of the beam energy is still available to power the GeV emission produced by inverse Compton upscattering of the Cosmic Microwave Background by the beam pairs.
Subject headings:
gamma rays: general – instabilities – intergalactic medium – plasmas – radiation mechanisms: nonthermal1. Introduction
With the current generation of erenkov telescopes, hundreds of TeV sources have been discovered. By far, the extragalactic TeV sky is dominated by blazars: jets from galactic centers beaming their emission towards our line of sight. The TeV photons from distant blazars cannot travel cosmological distances, since they interact with the extragalactic background light (EBL), producing electronpositron pairs. Studies of the attenuated GeVTeV light from distant blazars can therefore provide contraints on the strength of the EBL (e.g., Aharonian et al., 2006; Abdo et al., 2010).
The produced electronpositron pairs form a relativistic beam moving in the direction of the incident TeV photons. It is usually assumed that the energy of the pair beam is lost via inverse Compton (IC) scattering off the Cosmic Microwave Background (CMB). As a result, the TeV radiation will be reprocessed into the GeV band (Neronov & Semikoz, 2009). While cooling, the pairs gyrate around the IGM magnetic fields. Depending of the field strength and length scale, the GeV emission may form an extended source, show characteristic delays with respect to the TeV flux, or be strongly suppressed. These effects make combined GeV–TeV studies a useful probe of the IGM fields (e.g., Neronov & Vovk, 2010; Tavecchio et al., 2010; Dermer et al., 2011; Dolag et al., 2011; Taylor et al., 2011; Takahashi et al., 2012; Vovk et al., 2012).
Recently, it has been proposed that the destiny of the blazarinduced beams may be different. As they stream through the IGM plasma, the electronpositron pairs are expected to trigger collective plasma instabilities (as opposed to binary Coulomb collisions, that are negligible, as discussed by Miniati & Elyiv, 2013). For the parameters relevant to blazarinduced beams (i.e., dilute and ultrarelativistic), the fastest growing mode is the electrostatic oblique instability (e.g., Fainberg et al., 1970; Bret et al., 2010b), whose linear growth rate can exceed the IC cooling rate by several orders of magnitude (Broderick et al., 2012). Assuming that the instability keeps growing at the linear rate until all the beam energy is deposited into the IGM, the beam energy loss will be dominated by collective beamplasma instabilities, rather than IC cooling. In this case, the blazar TeV emission would not be reprocessed down to multiGeV energies, thus invalidating the IGM field estimates based on the GeV–TeV flux (Broderick et al., 2012, 2013). In addition, as a result of the beam relaxation, a substantial amount of energy would be deposited into the IGM. This “volumetric heating” can have dramatic consequences for the thermal history of the IGM (Chang et al., 2012; Pfrommer et al., 2012; Puchwein et al., 2012).
While plasma instabilities could, in principle, be fast enough to thermalize the pair beam, their nonlinear stages are far more complicated than what linear dispersion analysis predicts (e.g., the studies by Miniati & Elyiv 2013 and Schlickeiser et al. 2012b, 2013 reached opposite conclusions regarding the ultimate fate of blazarinduced beams). The nature of the fastest growing instability can change as the beamplasma system evolves, due to the suppression of temperaturesensitive modes as the beam heats up (e.g., Bret et al., 2008). Also, beamplasma instabilities can saturate at very small amplitudes, in particular for the extremely dilute beams produced by TeV blazars (e.g., Thode & Sudan, 1975; Thode, 1976).
In this work, we use firstprinciples particleincell (PIC) simulations in two and three dimensions to study the nonlinear stages and saturation of the instabilities generated as the blazarinduced pair beams propagate through the IGM. The nonlinear effects of the beamplasma interaction are hard to capture with analytical tools, and they require fullykinetic simulations. We explore a wide range of beam Lorentz factors and beamtoplasma density ratios , so that our results can be extrapolated to the extreme parameters of blazarinduced beams ( and , for the most powerful blazars). We find that, for ultrarelativistic dilute beams that start with a negligible thermal spread, electrostatic beamplasma instabilities can deposit of the beam energy into the background electrons. However, if the beam is born with a significant momentum dispersion (as expected for blazarinduced beams), the fraction of energy going into IGM heating is much smaller. We conclude that at least of the beam energy is still available to power the GeV emission produced by IC upscattering of the CMB. This lends support to the IGM magnetic field estimates that employ the combined GeV–TeV signature of distant blazars.
The paper is organized as follows. In §2 we derive the typical parameters of blazarinduced pair beams in the IGM. In §3 we describe the setup of our PIC simulations, whose results are presented in §4. In particular, in §4.1 we focus on one representative choice of beam parameters (in the regime of dilute ultrarelativistic cold beams) and we describe the complete evolution of the beamplasma unstable system, from the early exponential phase up to the nonlinear stages. In §4.2, we discuss the dependence of our findings on the beam Lorentz factor, the beamtoplasma density contrast and the beam temperature. For the convenience of readers uninterested in the kinetic details of beamplasma instabilities, in §4.3 we summarize our results in application to blazarinduced beams. Finally, in §5 we assess the implications of our findings for the thermal history of the IGM and the detection of reprocessed GeV emission from powerful TeV blazars.
2. Physical Parameters of BlazarDriven Beams
In this section, we summarize the physical parameters of blazarinduced beams, including the density contrast to the IGM, the beam Lorentz factor and velocity spread. We present orderofmagnitude estimates, and we refer to Schlickeiser et al. (2012a) and Miniati & Elyiv (2013) for a more detailed analysis of the beam distribution function, and its dependence on the spectrum of the EBL and of the blazar TeV emission.
Blazar photons of energy travel a distance of
before they interact with the EBL and produce electronpositron pairs (Neronov & Semikoz, 2009). Here, accounts for uncertainties in the intensity of the EBL, with models predicting
for (e.g., Aharonian, 2001).
For a powerful blazar with TeV isotropic equivalent luminosity (e.g., Ghisellini et al., 2010), the number density of the beam pairs is set by the balance of the pair production rate with the energy loss rate (dominated by IC cooling), which gives
(1) 
where the factor of accounts for the rapid energy loss of the pairs due to IC.
(2) 
The density contrast and the beam Lorentz factor are the two crucial parameters determining the plasma physics of the beamIGM interaction. For blazarinduced beams, we expect that they should vary in the range and , respectively.
As discussed by Broderick et al. (2012), collective beamplasma effects can be relevant only if many beam pairs are present within a sphere of radius equal to the wavelength of the most unstable mode. As we show below, the scale of the fastest growing modes is , where is the plasma skin depth of the IGM electrons (with number density ). A sphere of skindepth radius contains beam particles. For , we find that collective phenomena always play a role in the evolution of blazarinduced beams.
Another important parameter is the dispersion in beam momentum at birth. Since the pair creation cross section peaks slightly above the threshold energy, the pairs are born moderately warm (with a comoving temperature of ). Moreover, since the EBL and the blazar TeV spectra are broad, the beam energy distribution will extend over a wide range of Lorentz factors, as discussed by Miniati & Elyiv (2013). In §4.2.2, we explore the role of thermal effects on the nonlinear evolution of blazarinduced beams.
3. Simulation Setup
We investigate blazardriven plasma instabilities in the IGM by means of fullykinetic PIC simulations. We employ the threedimensional (3D) electromagnetic PIC code TRISTANMP (Spitkovsky, 2005), which is a parallel version of the publicly available code TRISTAN (Buneman, 1993), that was optimized for handling ultrarelativistic flows (see, e.g., Sironi & Spitkovsky 2011 and Sironi et al. 2013 for studies of ultrarelativistic collisionless shocks using TRISTANMP). We initialize a relativistic dilute pair beam that propagates along through an unmagnetized electronproton plasma (with the realistic mass ratio ). The simulations are performed in the frame of the background plasma, i.e., of the IGM. No background magnetic field is assumed, so the electric and magnetic fields generated by beamplasma instabilities will grow from noise.
To follow the beamplasma evolution to longer times with fixed computational resources, we mainly utilize 2D computational domains in the plane. In §4.1 we compare 2D and 3D runs, and we show that 2D simulations can capture most of the relevant 3D physics. In the case of 2D simulations with the beam lying in the simulation plane, only the inplane components of the velocity, current and electric field, and only the outofplane component of the magnetic field are present. The simulation box is periodic in all directions. By choosing a periodic domain, we simulate the bulk of the beamplasma system, rather than the “head” of the pair beam.
The background plasma consists of cold electrons and protons, with initial electron temperature . We have tested that higher temperatures of the background electrons do not significant change the development of the relevant instabilities, as long as the electron temperature is nonrelativistic, in agreement with Bret et al. (2005). The background protons are allowed to move, but we obtain similar results when the protons are treated as a static chargeneutralizing background.
The beam consists of electronpositron pairs propagating with Lorentz factor along the direction. To our knowledge, our PIC simulations are the first to address the evolution of an electronpositron beam. All of the previous studies have focused on the case of an electron beam propagating through an electronproton plasma, with the background electrons moving opposite to the beam to compensate for the beam current (e.g., Dieckmann et al., 2006b; Bret et al., 2008; Kong et al., 2009). Any instability triggered by the relative drift between the background electrons and protons (e.g., the Buneman (1958) instability) will then be absent in our setup, where the pair beam carries no net current.
The beamtoplasma density ratio and the beam Lorentz factor expected for blazarinduced pairs streaming through the IGM (see §2) cannot be directly studied with PIC simulations. Yet, by performing dedicated experiments with a broad range of and (in the regime and of ultrarelativistic dilute beams), we can extrapolate the relevant physics to the extreme parameters expected in the IGM. We vary the beam Lorentz factor from up to , and the density contrast from down to .
The results presented below have been extensively tested for convergence. We typically employ 50 particles per computational cell for the background plasma (25 electrons and 25 protons), and the same number for the beam particles (if each carries a weight ). However, we have tested that our results are the same when using up to 256 particles per cell, for both the beam and the plasma. We resolve the skin depth of the background electrons with 8 computational cells, but we have tested that our results do not change when using or cells per skin depth.
The number of beam particles is kept constant during the evolution of the beamplasma system, since the photonphoton interactions that would introduce fresh electronpositron pairs are extremely rare on the timescales covered by our simulations. Also, we neglect IC cooling of the beam pairs, since it is irrelevant over the timespan of our runs. This implies that the total energy in our periodic beamplasma system should be constant over time. However, explicit PIC codes do not conserve energy to machine precision. We track the energy conservation in our runs, and we find that at late times it is still better than . This makes our estimates of the amount of beam energy transferred to the plasma electrons (of order ) extremely robust, for the beam parameters explored in this work.
Finally, we remark that in all PIC codes a numerical heating instability arises when cold relativistic plasma propagates for large distances over the numerical grid (Dieckmann et al., 2006a). Since the numerical speed of light on the grid is smaller than the correct value at large wavenumbers, ultrarelativistic particles will emit numerical erenkov radiation. This might artificially slow down the beam, even in the absence of physical beamplasma instabilities. We have assessed that the results reported below arise from a physical instability (as opposed to the numerical erenkov mode), by comparing our beamplasma simulations with the artificial case of a beam that propagates through the grid in the absence of any background plasma. The beam evolution in the two cases is dramatically different, which provides further confirmation that the beam energy loss that we discuss below arises from the physical interaction of the beam with the background plasma, rather than from the numerical erenkov instability.
4. Results
In this section, we explore the linear and nonlinear evolution of ultrarelativistic dilute pair beams by means of 2D and 3D PIC simulations. In §4.1, we describe the different stages of evolution of the beamplasma system, for a representative choice of beam parameters in the regime of ultrarelativistic dilute beams (, and negligible beam thermal spread at initialization). In §4.2, we investigate the dependence of our results – in particular, of the fraction of beam energy transferred to the background electrons – on the beamtoplasma density contrast, the beam Lorentz factor and the beam temperature at birth. The reader that is not interested in the kinetic details of the beamplasma interaction might proceed to §4.3, where we extrapolate the findings of our PIC simulations to the extreme parameters of blazarinduced beams.
4.1. The Linear and NonLinear Evolution of UltraRelativistic Dilute Cold Beams
In this section, we follow the evolution of a cold beam with and . We start with the analysis of the linear phase, and then we investigate the nonlinear relaxation. We find that the exponential phase of the oblique mode (which is the fastest growing instability for dilute ultrarelativistic beams) terminates due to selfheating of the beam in the direction transverse to the beam motion. At the end of the oblique phase, only a minor fraction of the beam energy has been deposited into the background electrons. Further evolution of the beam is governed by quasilongitudinal modes, which operate on a timescale that is much longer (at least two orders of magnitude) than the oblique growth. At the end of the quasilongitudinal phase, the dispersion of beam momentum in the longitudinal direction saturates at , which corresponds to a fraction of beam energy transferred to the plasma.
The Oblique Exponential Phase
The evolution of the beamplasma system during the oblique phase is presented in panels (a)(d) of Figure 1. The beam is set up with a small thermal spread (), so that the oblique instability initially proceeds in the reactive regime, i.e., all the beam particles are in resonance with each harmonic of the packet of unstable modes, and the instability is the strongest. This is opposed to the kinetic regime, in which the beam velocity spread is considerable. Here, a number of unstable modes with a broad spectrum in phase velocity will be excited, with only a small number of beam particles being in resonance with each mode. This results in a slower growth, as compared to the reactive regime.
In Figure 1, we confirm that the oblique instability is the fastest growing mode for ultrarelativistic dilute beams. The reactive phase of the instability governs the evolution of the system for , where is marked as a dashdotted vertical blue line in panels (a) and (b). Here, is the plasma frequency of the background electrons. The fastest mode grows on a scale , where is the electron skin depth, and its wavevector is inclined at with respect to the beam propagation. This is apparent in the 2D structure of the longitudinal electric field in Figure 1(c), as well as in the 2D plots of the transverse electric and magnetic fields (not shown).
The growth rate of the oblique instability in the reactive regime is (e.g., Fainberg et al., 1970)
(3) 
where the different dependence on and is related to the fact that for relativistic beams the transverse inertia is much smaller than the longitudinal inertia (by a factor of ), so that the modes transverse to the beam are the easiest to be excited (for an intuitive physical description, see Nakar et al. 2011).
(4) 
which nicely agrees with our results. In fact, in Figure 1(a) we show that the fraction of beam kinetic energy deposited into the background electrons (orange line), into the longitudinal (red) and transverse (blue) electric fields, and into the transverse magnetic field (green) all grow at the rate predicted by Equation (4) for (dotted orange line in Figure 1(a)).
The oblique mode is quasielectrostatic, i.e., roughly (e.g., Bret et al., 2010b). Since the angle between the wavevector and the beam is , in agreement with analytical expectations (e.g., Bret et al., 2010b), it follows that the electric field component tranverse to the beam is slightly larger than the longitudinal component, i.e., implies that . This explains the small difference between the red and the blue lines in Figure 1(a). Also, since the mode is quasielectrostatic, the magnetic component will be subdominant relative to the electric fields. In agreement with the analytical considerations of Lemoine & Pelletier (2010), we find that . Given that for ultrarelativistic dilute beams, it follows that .
For electrostatic modes in the linear phase, it is expected that the amount of kinetic energy lost by the beam should be equally distributed between electric fields and plasma heating (e.g., Thode, 1976). This explains why the energy of the background electrons in the exponential phase (orange line in Figure 1(a) at ) is comparable to the energy in electric fields ( and , respectively red and blue lines in Figure 1(a)). We have verified that the fraction of beam energy transferred to the background protons is negligible as compared to the plasma electrons, so the evolution of the protons will be ignored hereafter.
Since during the exponential phase of the oblique mode, both the plasma and the beam are heated quasiisotropically, so that the momentum spreads in the longitudinal and transverse directions are nearly identical (compare the red and blue lines for in Figure 1(b); dashed lines refer to the plasma, solid lines to the beam). The transverse momentum spread of the beam can be related to the transverse electric field via the Lorentz force
(5) 
where we have assumed that the characteristic timescale is set by the oblique growth rate (i.e., ) and that the magnetic force is negligible compared to the electric force (in fact, ).
Since the oblique mode is heating up the beam in the transverse direction (solid blue line in Figure 1(b)), the exponential growth at the reactive rate will necessarily terminate, when the assumption of a cold beam required by the reactive approximation becomes invalid. It is well known that the system will transition from the reactive phase to the kinetic phase when the beam velocity dispersion reaches (e.g., Fainberg et al., 1970; Bret et al., 2010b)
(6) 
namely when the beam, due to its velocity spread, can move across one wavelength of the most unstable mode during the growth time of the reactive instability. In this case, most of the beam particles will lose resonance with the unstable mode, and the instability will transition from the reactive to the kinetic regime. For ultrarelativistic beams with isotropic momentum dispersions (in fact, Figure 1(b) shows that during the oblique reactive phase), the transverse velocity spread is much larger than the longitudinal spread . Since , Equation (6) above reduces to
(7) 
where is the expected transverse dispersion in beam momentum at the end of the reactive oblique phase. The threshold in momentum dispersion can also be recast as a limit in beam temperature (e.g., Bret et al., 2010a). We confirm that the reactive phase of the oblique mode terminates at , when the beam transverse momentum reaches the threshold in Equation (7) (which is shown as a horizontal dashdotted blue line in Figure 1(b)).
We can now derive the expected fraction of beam kinetic energy transferred to the plasma electrons and to the electromagnetic fields at the end of the oblique reactive phase. By setting in Equation (5), we find that the fraction of beam energy converted into transverse electric fields at is
(8) 
Since the oblique mode is quasielectrostatic (i.e., ), it follows that . Moreover, for electrostatic modes, the fraction of the beam kinetic energy converted into plasma heating is comparable to the energy in electric fields (e.g., Thode, 1976), so . Finally, since , the magnetic energy fraction will be . We have extensively verified that the expected scalings of the efficiency parameters and with respect to are in agreement with the results of our simulations, across the whole range of beam Lorentz factors and density contrasts we have explored (see the various curves in Figure 1(a) at , and also §4.2.1).
For , the evolution of the oblique mode will proceed in the kinetic (rather than reactive) regime. The kinetic oblique mode is indeed resposible for the peak in the electric field energy observed at in Figure 1(a) (red and blue lines), which produces a moderate increase in the fraction of beam energy transferred to the background electrons (orange line in Figure 1(a) at ). In this phase, the 2D structure of the longitudinal electric field in Figure 1(d) shows that the wavevector of the kinetic oblique mode is oriented at relative to the beam propagation (as compared to the angle observed during the reactive phase, see Figure 1(c)). A similar pattern is shown in the 3D plot of Figure 2(b). As expected, the increase in the transverse momentum dispersion has suppressed the modes having , which are most sensitive to transverse temperature effects (e.g., Bret et al., 2010b).
For a beam with initial transverse velocity dispersion , the growth rate of the kinetic oblique mode for is (e.g., Breǐzman & Ryutov, 1971)
(9) 
At the end of the reactive oblique phase, Equation (7) prescribes that , so that the growth rate of the kinetic oblique mode will be , i.e., it will have the same scalings with and as the reactive oblique mode (here, we have neglected factors of order unity).
We have explicitly verified that the peak in electric fields at is due to the kinetic oblique mode, by performing a dedicated simulation in which at (i.e., shortly after the end of the reactive stage) we reset by hand the electromagnetic fields and the plasma temperature to their initial values (i.e., no seed fields and ), yet we retain the beam momentum distribution that results selfconsistently from the reactive oblique phase. In this setup, we find that the fastest growing mode has the same 2D pattern as in Figure 1(d) and its growth rate scales as . This confirms that the peak in electric fields at (Figure 1(a)) is indeed associated to the kinetic oblique instability.
The growth of the kinetic oblique mode terminates due to selfheating of the beam, in analogy to the reactive oblique phase. In particular, the growth in the kinetic oblique phase cannot be sustained beyond the point where, due to the selfexcited electric fields, the beam momentum dispersion in the transverse direction exceeds the initial value . At this point, the expression for the growth rate in Equation (9) becomes clearly invalid. This happens when
(10) 
which leads to the same scaling as in Equation (8), if we take , as appropriate for the beam velocity dispersion at the end of the reactive oblique phase. In summary, apart from factors of order unity, the electric fields at the end of the kinetic oblique phase will saturate at a level similar to the reactive oblique stage (compare the two peaks of electric energy in Figure 1(a) at and ). It follows that the kinetic oblique instability will increase the fraction of beam kinetic energy transferred to the plasma electrons only by a factor of order unity, as compared to the reactive oblique mode (see the orange line in Figure 1(a), for ).
The Longitudinal Relaxation Phase
After the saturation of the oblique instability, further evolution of the beamplasma system proceeds via quasilongitudinal modes, as shown in the 2D plot of the longitudinal electric field in Figure 1(g), as well as in the 3D pattern of Figure 2(c). Being quasilongitudinal, these modes are insensitive to thermal spreads in the direction perpendicular to the beam, so they can grow even after the end of the kinetic oblique phase.
The quasilongitudinal oscillations shown in Figure 1(g) are a characteristic signature of the quasilinear relaxation of the beam (see, e.g., Grognard, 1975; Lesch & Schlickeiser, 1987; Schlickeiser et al., 2002; Pavan et al., 2011). In the quasilinear relaxation, the beam generates longitudinal Langmuir waves (see the peak in at in Figure 1(e)), which scatter the beam particles and heat the background plasma (see the increase in the electron thermal energy shown by the orange line at in Figure 1(e)). Since (compare the red and blue curves in Figure 1(e) at ), the background electrons will be heated preferentially in the direction of motion of the beam (see the increase in at in Figure 1(f)). For the same reason, the quasilinear relaxation is accompanied by a substantial increase in the beam momentum spread along the direction of propagation (red solid line in Figure 1(f), showing the growth of at ). The spread in the parallel beam momentum is associated to the formation of phase space holes, that result from the trapping of beam particles by the longitudinal electric oscillations (e.g., O’Neil et al., 1971; Thode & Sudan, 1975).
The quasilinear relaxation occurs on a timescale much longer than the exponential oblique phase. Within the range of beam Lorentz factors and density contrasts probed by our simulations, we find that the characteristic relaxation time is is at least two orders of magnitude longer than the exponential growth time of the oblique instability , in agreement with previous 1D simulations (Grognard, 1975; Pavan et al., 2011). This emphasizes the importance of evolving our PIC simulations to sufficiently long times to capture the physics of the quasilinear relaxation (see §4.2 for further details).
The quasilinear modes broaden the beam momentum spectrum in the longitudinal direction up to the point where (see the red solid line in Figure 1(f), saturating at ). This in agreement with the socalled Penrose’s criterion, stating that the beam will be stable to electrostatic modes only when the longitudinal dispersion in momentum approaches the initial beam Lorentz factor (e.g., Buschauer & Benford, 1977). From , it follows that at the end of the relaxation phase, a fraction of the beam energy has been transferred to the background electrons (see the orange line in Figure 1(e) at ). In §4.2.1 we demonstrate that, irrespective of the beam Lorentz factor or the beamtoplasma density contrast, a generic byproduct of the relaxation of cold ultrarelativistic dilute beams is the conversion of of their energy into plasma heating.
The quasilinear relaxation significantly affects the shape of the beam and plasma longitudinal momentum spectrum, as shown in Figure 3. As a result of the quasilongitudinal relaxation, the plasma distribution at late times () develops a pronounced highenergy tail (at ), that bridges the main thermal peak of the plasma electrons (at ) with the beam particles (that populate the isolated highenergy bump at in Figure 3). This highenergy component in the background electrons, which contains a significant amount of energy, is present only in the forward direction (i.e., along the beam propagation). In the backward direction, the spectrum of plasma electrons (red dotted line in Figure 3, at ) is compatible with a Maxwellian.
During the quasilinear relaxation, the beam spectrum evolves from a quasimonoenergetic distribution into a broad bump (from the black to the red curve at in Figure 3). In agreement with Figure 1(f) (red solid line), most of the evolution occurs at , whereas the beam spectrum at longer times is remarkably steady (we have followed the system up to , finding no further signs of evolution).
We point out that the beam momentum spectrum at late times does not approach the socalled “plateau” distribution (indicated as a black dotted line in Figure 3), which is believed to be the ultimate outcome of the beam relaxation in 1D (e.g., Grognard, 1975; Schlickeiser et al., 2002). As opposed to earlier 1D claims, in our 2D and 3D simulations we find that the beam longitudinal relaxation leads to a momentum spectrum that is harder than the plateau distribution, yet the beamplasma system appears stable.
We argue that the transverse dispersion in beam momentum, which could not be properly captured in previous 1D studies, prevents the longitudinal beam spectrum from relaxing to the plateau distribution (see Appendix B for further details). Our claim is supported by the following experiment. At the end of the kinetic oblique phase (), we artificially set the beam transverse dispersion to be (for comparison, the selfconsistent evolution in Figure 1(b) yields at the end of the kinetic oblique phase, see the solid blue line). Also, in the subsequent evolution, we inhibit any growth in the transverse beam momentum. In this setup, in which any transverse dispersion effects are artificially neglected, the beam relaxation should proceed as in 1D. So, it is not surprising that, in agreement with previous 1D studies, at late times the beam relaxes to the plateau distribution (red dashed line in Figure 3). In other words, we find that the relaxation to the plateau distribution is not a general result of the multidimensional evolution of ultrarelativistic beams, but it only occurs when the beam transverse dispersion stays sufficiently small, so the system is quasi1D.
As we further discuss in Appendix B, the beam relaxation produces a plateau distribution only if , due to the following argument. Complete stabilization of the beamplasma system is achieved when the longitudinal velocity spread of the ultrarelativistic beam reaches (more precisely, when the velocity spread is comparable to the beam speed). The spread in longitudinal velocity includes contributions from both the longitudinal and the transverse momentum dispersions:
(11) 
where the second term on the right hand side is absent in the case of 1D relaxation. It follows that the transverse momentum spread can appreciably modify the relaxation process only if , where we have used that at the end of the relaxation phase. Then, the fact that in the case studied in Figure 1(b) the transverse beam dispersion after the oblique phase is explains why the beam relaxation cannot lead to a plateau distribution. In Appendix B, we provide further evidence that the transition to a plateau distribution requires .
The Magnetic Field Growth
In the previous subsections, we have primarily focused on the electrostatic character of the growing modes, which determines the coupling efficiency between the beam energy and the plasma thermal energy. Here, we comment on the generation of magnetic fields associated with the evolution of the beamplasma system.
As shown in Figure 1(a), the growth of the quasielectrostatic oblique mode (both in the reactive and in the kinetic regime) is accompanied by a minor magnetic component. In §4.1.1, we have estimated that the fraction of beam kinetic energy transferred to the magnetic fields at the end of the oblique phase is , apart from factors of order unity. Since for ultrarelativistic dilute beams, the magnetic fields generated by the oblique instability are generally unimportant.
At the end of the relaxation phase, the beam and the plasma are highly anisotropic, with the longitudinal momentum spread much larger than the transverse one (see Figure 1(f) at ). As a result, the system is prone to the Weibel instability (e.g., Weibel, 1959; Yoon & Davidson, 1987; Silva et al., 2002), which generates the transverse magnetic field pattern shown in Figure 1(h).
The Weibel instability is predominantly magnetic, so it does not mediate any significant exchange of energy from the beam to the plasma electrons. Yet, it might be a promising source for the generation of magnetic fields, as discussed by Schlickeiser et al. (2012b). However, the evolution of the magnetic filaments shown in Figure 1(h) can only be captured with largescale 3D simulations (e.g., Bret et al., 2008; Kong et al., 2009). A detailed 3D investigation of the strength and scale of the magnetic fields resulting from ultrarelativistic dilute pair beams from TeV blazars will be presented elsewhere.
4.2. Dependence on the Beam Parameters
In this section, we discuss the dependence of the beamplasma evolution on the beam parameters. In §4.2.1, we consider the case of cold beams, and we show that the quasilongitudinal relaxation leads to a beam momentum spread along the direction of motion , regardless of the beam Lorentz factor or the beamtoplasma density contrast . In turn, this implies that a fraction of the beam energy is converted into heat of the background plasma, irrespective of or .
In 4.2.2, we discuss the effect of the initial beam thermal spread on the efficiency of the beamtoplasma energy transfer. We find that if the initial dispersion in longitudinal momentum satisfies (as typically expected for blazarinduced beams), the fraction of beam energy deposited into the background plasma is much smaller than .
Cold Beams
In Figure 4 we show how the evolution of the beamplasma system depends on the beam Lorentz factor (that we vary from up to ) and on the beamtoplasma density contrast (from down to ). Most of the previous studies have focused on moderately relativistic electron beams () with (e.g., Gremillet et al., 2007; Bret et al., 2008; Kong et al., 2009). Here, we extend our investigation to the case of ultrarelativistic dilute electronpositron beams, as appropriate for blazarinduced beams.
For each choice of and , we follow the beamplasma system from the oblique phase until the quasilinear relaxation (typically, up to ). We initialize a cold beam with thermal spread , so that the oblique instability initially proceeds in the reactive regime, for the range of and covered by our simulations. In the reactive phase, we confirm that the fastest growing mode has a wavevector oriented at to the beam direction of propagation. The growth rate is in excellent agreement with Equation (4). As described in §4.1.1, the exponential phase of the reactive oblique mode terminates due to selfheating of the beam in the transverse direction. At the end of the reactive phase, we find that the fractions of beam kinetic energy transferred to the plasma electrons, to the electric fields and to the magnetic fields scale respectively as , and , as in §4.1.1.
The reactive oblique phase is followed by the kinetic oblique phase. We find that the characteristic wave pattern in the kinetic oblique phase (see Figure 1(d) in 2D and Figure 2(b) in 3D) appears in the evolution of all the beamplasma systems that we present in Figure 4. In the temporal evolution of the fraction of beam energy converted into plasma heating, the kinetic oblique mode is responsible for the additional increase that is seen in the most relativistic cases () after the end of the reactive oblique phase (see the insets in the top row of Figure 4). Regardless of or , we find that the kinetic oblique phase deposits only a fraction of the beam kinetic energy into the background electrons, i.e., comparable to the reactive oblique phase (see §4.1.1).
The longterm evolution of the system is controlled by the quasilongitudinal relaxation, which operates on a timescale , where is the characteristic efolding time of the oblique mode. In Figure 4, the quasilinear relaxation governs the growth in the plasma thermal energy (top row) and in the beam parallel momentum spread (middle row) occurring at . In the regime , the quasilinear relaxation terminates when the beam momentum spread in the longitudinal direction reaches , regardless of or (middle row in Figure 4). Correspondingly, the fraction of beam energy transferred to the plasma saturates at (top row in Figure 4).
Mildly relativistic beams with moderate density contrasts deviate from such simple scalings, for the following reason. The quasilinear relaxation will not operate if the beam dispersion at the end of the oblique phase is already . According to Equation (7), the beam spread at the end of the oblique phase is , so that the quasilinear relaxation will be suppressed if , i.e., for mildly relativistic beams with moderate .
A similar argument explains why a plateau distribution in the longitudinal momentum spectrum (bottom row in Figure 4) is established only for relatively small . As we have argued in §4.1.2, a transverse spread prevents the beam relaxation to the plateau distribution (shown as a dotted black line in the bottom row of Figure 4). At the end of the oblique phase, , so that the beam will relax to the plateau distribution only if . Clearly, this constraint is hardest to satisfy for highly relativistic beams, which explains why, at fixed , the momentum spectrum of more relativistic beams shows stronger deviations from the plateau distribution.
Hot Beams
In the previous subsection, we have assumed that the beam is born with a negligible thermal spread. Here, we discuss how the results presented above for cold beams will be modified by temperature effects. We describe separately the role of thermal spreads in the development of the oblique instability, and in the relaxation phase.
As we have anticipated in §4.1.1, the oblique instability will proceed in the kinetic (rather than reactive) regime if the initial beam dispersion in the transverse direction is .
(12) 
which results in a fraction of the beam energy transferred to the electric fields at the end of the kinetic phase ( is defined in Equation (9)). Since the kinetic oblique mode is an electrostatic instability, the fraction of beam energy converted into heat will also be . Moreover, due to the fact that – they are comparable only if , i.e., at the boundary between reactive and kinetic regimes – we expect the kinetic oblique mode to be less efficient in heating the plasma electrons, as compared to the reactive phase. For beams with and , we have indeed verified with PIC simulations (not presented here) that the exponential growth of the kinetic oblique mode terminates at smaller for larger values of , in good agreement with the expected scaling (at fixed and ).
Regardless of the character of the oblique mode (reactive or kinetic), the longterm evolution of the beamplasma system is controlled by the quasilinear relaxation. In Figure 5, we show how the relaxation phase is affected by a finite beam temperature . For the set of beam parameters employed in Figure 5 ( and ), the oblique phase is expected to occur in the reactive regime, as long as the beam comoving temperature is nonrelativistic. In fact, for all the choices of presented in Figure 5 (in the plot, is in units of ), the early increase in the heating efficiency proceeds at the reactive rate (compare the curves in the inset of Figure 5(a) with the dotted red line, that scales with the oblique growth rate). Also, the kinetic oblique instability, which is responsible for the further growth in at (see the inset in Figure 5(a)), does not show any dependence on temperature, in the range explored in Figure 5.
The beam temperature has profound effects on the quasilinear relaxation phase, for the beam parameters employed in Figure 5. For cold beams (, yellow and red lines in Figure 5), the relaxation phase does not depend on the beam temperature. In agreement with the results presented in §4.2.1, the longitudinal spread in the beam momentum increases during the relaxation stage until , as shown in Figure 5(b). This corresponds to a fraction of the beam kinetic energy being converted into plasma heating (Figure 5(a)). Similar conclusions hold for moderate beam temperatures (, green line), whereas the quasilinear relaxation is suppressed if the beam temperature at birth is such that the initial beam spread (cyan, blue and black lines in Figure 5(b)). In this case, the quasilinear phase does not mediate any further increase in the heating efficiency , beyond the early oblique phase (see the black line in Figure 5(a)).
The longitudinal momentum spectrum in Figure 5(c) clarifies why the quasilinear relaxation is suppressed for hot beams. As we have discussed in §4.1.2, the relaxation is mediated by Langmuir waves excited by the beam. The beam particles that are slightly faster than the wave tend to transfer energy to the wave, and thus excite it, while those that are slightly slower than the wave tend to receive energy from it, and thus damp the wave. It follows that a mode is unstable if the number of beam particles moving slightly faster than the wave exceeds that of those moving slightly slower. More precisely, the instability will be stronger for a harder slope of the longitudinal momentum spectrum at momenta (this is the relevant momentum scale, since the phase velocity of the most unstable mode is slightly smaller than the beam speed). For a sharply peaked beam (i.e., with ), the slope will be extremely hard, and the excitation of the Langmuir modes that mediate the relaxation process will be most efficient. As a result of the quasilinear relaxation, the beam particles will be scattered down to lower energies, and when the beam momentum spread exceeds the slope of the momentum spectrum below the peak becomes too shallow (see the red and yellow curves in Figure 5(c) at ), and the quasilinear relaxation terminates. If the beam temperature at birth is too hot (i.e., such that ), the initial slope of the momentum spectrum at is already too shallow to trigger efficient excitation of Langmuir waves, and the quasirelaxation process is suppressed.
In Figure 5(c), we further support this argument by showing that, in the latest stages of evolution, the shape of the beam momentum spectrum below (in Figure 5(c), at ), is nearly independent of the beam temperature at birth. For initially cold beams (yellow and red lines), the quasilinear relaxation increases the beam momentum spread over time (see Figure 5(b)), until the beam spectrum relaxes to the broad bump shown in Figure 5(c). At this point, the beam is stable, although it has not relaxed to the socalled plateau distribution (as we explain in §4.1.2, the relaxation to the plateau distribution requires , which is not satisfied here). Below , the initial momentum spectrum of hot beams resembles the final momentum distribution of initially cold beams. This explains why hot beams having at birth will not experience the quasilinear relaxation phase.
So far, we have discussed the effect of thermal spreads on the beam relaxation, assuming that the beam spectrum at birth is a drifting Maxwellian. However, the conclusions derived above hold for more complicated beam distributions. In particular, we have performed a set of PIC simulations, assuming that the beam spectrum at birth is a power law for . For blazar environments, such a beam spectrum is expected as a result of intense IC cooling off the CMB.
If the beam powerlaw distribution has a negligible spread in the transverse direction (which might not be the case for blazarinduced beams, see §4.3), then the oblique instability will proceed in the reactive regime. Apart from factors of order unity, we find that the exponential growth rate is , similar to the case of monoenergetic beams. As we have discussed above, the effectiveness of the quasilinear relaxation – which ultimately determines whether the heating efficiency can reach – is determined by the spread in the beam longitudinal spectrum below the peak, i.e., at . If the beam spectrum has a sharp lowenergy cutoff at , then the quasilinear relaxation does operate, and the beam deposits of its energy into the background electrons. However, if the lowenergy end of the beam distribution is broader, the quasilinear relaxation will be inhibited, and the amount of beam energy transferred to the plasma will be much smaller, in agreement with the results in Figure 5. We find that, if the beam longitudinal spectrum below can be modeled as a power law , the quasilinear relaxation is suppressed if , with the case corresponding to the plateau distribution.
4.3. Implications for BlazarInduced Beams
We now analyze the implications of our findings for the evolution of blazarinduced beams in the IGM. As we have discussed in §3, the Lorentz factors and density contrasts of blazardriven beams ( and ) cannot be directly studied with PIC simulations. However, by performing a number of experiments with a broad range of and , we have been able to assess how the relaxation of ultrarelativistic dilute beams depends on the beam Lorentz factor and the beamtoplasma density ratio. Our results can then be confidently extrapolated to the extreme parameters of blazarinduced beams.
We have demonstrated that the oblique instability, which governs the earliest stages of evolution of ultrarelativistic dilute beams, proceeds in the kinetic regime if the initial velocity dispersion in the direction transverse to the beam is , where is the growth rate of the reactive oblique mode in units of the plasma frequency rads of the IGM electrons. Since the beam pairs are born with a mildly relativistic thermal spread in the center of mass of the photonphoton interaction (see §2), the transverse velocity spread in the IGM frame will be . If follows that the oblique instability will proceed in the kinetic regime if , which is marginally satisfied for blazarinduced beams (). The oblique instability will grow at the kinetic rate in Equation (9), which for , as appropriate for blazarinduced beams, reduces to . For the parameters of blazarinduced beams, this is a factor of
(13) 
larger than the IC cooling rate , where is the IC cooling length computed in §2. In agreement with Broderick et al. (2012) and Schlickeiser et al. (2012b), we find that the kinetic oblique instability has ample time to grow, before the beam loses energy to IC emission. Furthermore, the typical lifetime of blazar activity of is sufficient for the instability to operate.
Due to selfheating of the beam in the transverse direction (see §4.1 and 4.2.2), the exponential phase of the kinetic oblique instability terminates when only a minor fraction of the beam energy has been transferred to the IGM electrons. As we have argued in §4.2.2, the heating efficiency of blazarinduced beams at the end of the kinetic oblique phase is only
(14) 
It follows that, even though the oblique instability can grow faster than the IC cooling time, its efficiency in heating the plasma electrons is extremely poor.
As we have emphasized in §4.1.2, a larger amount of beam energy (up to ) can be deposited into the plasma electrons by the quasilinear relaxation phase. We find that the quasilinear relaxation occurs on a timescale much longer than the growth time of the oblique instability, . Yet, since the oblique growth rate is much faster than the IC cooling rate, see Equation (13), the quasilinear relaxation should have enough time to operate before the beam energy is lost to IC emission. Here, we are conservatively neglecting the possibility that might be much longer than for more extreme beam parameters, as a result of nonlinear plasma processes that reduce the strength of the electric fields available for the quasilinear relaxation (see Schlickeiser et al., 2012b).
In §4.2.2, we have demonstrated that the quasilinear relaxation can occur only if the beam momentum spectrum along the longitudinal direction is sufficiently narrow. More precisely, we find that if the beam distribution peaks at (in the case of a powerlaw tail, this would be the lowenergy cutoff), the quasilinear relaxation can operate only if the parallel momentum spread below the peak satisfies . This constraint is hard to fulfill by blazarinduced beams. Since the pair creation cross section peaks slightly above the threshold energy, the pairs are born moderately warm, with a comoving temperature . This corresponds to a longitudinal momentum spread at birth (measured in the IGM frame) of , which is already prohibitive for the development of the quasilinear relaxation. Moreover, the momentum dispersion of blazarinduced beams might be even larger, since the spectrum of both the EBL and the blazar TeV emission are usually modeled as broad power laws (Miniati & Elyiv, 2013).
In summary, a solid assessment of the shape of the beam momentum distribution below the peak is essential to predict the amount of beam energy deposited into the IGM. If the beam spectrum at birth were to have (which is unlikely to be the case, but see Schlickeiser et al. 2012b), then the quasilinear relaxation would transfer of the beam energy to the IGM. Even in this optimistic case, we remark that the heating efficiency would reach at most , so of the energy would still remain in the beam. At the end of the quasilinear relaxation phase, the beam spectrum will be harder than the socalled plateau distribution, since the transverse dispersion in beam velocity does not meet the requirement for relaxation to the plateau spectrum. In the more realistic case , the quasilinear relaxation will be inhibited, resulting in a lower efficiency of IGM heating – with a firm lower limit being the heating fraction at the end of the oblique phase, see Equation (14).
Comparison with Earlier Studies
We now compare our numerical work with earlier analytical studies of the relaxation of blazarinduced beams in the IGM. As we have emphasized in §3, our PIC simulations are the first to address the evolution of dilute ultrarelativistic electronpositron beams, as appropriate for blazarinduced beams in the IGM. Most of the previous PIC studies (e.g., Dieckmann et al., 2006b; Gremillet et al., 2007; Bret et al., 2008; Kong et al., 2009) have focused on mildly relativistic electron beams.
In agreement with Broderick et al. (2012) and Schlickeiser et al. (2012a), we find that the fastest growing instability for blazarinduced beams propagating through the IGM is the oblique mode. If the pair beam were to be initially cold, the oblique instability would evolve at the reactive rate in Equation (4), with wavevector oriented at relative to the beam (see Schlickeiser et al. 2012a). However, the initial transverse spread in beam momentum is large enough such that the oblique mode evolves at the kinetic rate in Equation (9), as argued by Broderick et al. (2012) and Miniati & Elyiv (2013). We remark that the transverse beam spread does affect the growth of the oblique mode, but it will not impact the evolution of longitudinal waves, as found by Schlickeiser et al. (2013). Yet, since the early evolution of blazarinduced beams is controlled by the oblique (rather than longitudinal) mode, transverse thermal effects are indeed important for the beamplasma interaction at early times.
The evolution of blazarinduced beams at late times is governed by nonlinear plasma processes, which are extremely hard to capture with analytical tools. Schlickeiser et al. (2012b) and Miniati & Elyiv (2013) attempted to describe the nonlinear relaxation of blazarinduced beams, reaching opposite conclusions regarding the ultimate fate of the beam energy. Schlickeiser et al. (2012a) assumed that the beam momentum distribution can be modeled as a delta function (i.e., they approximated the beam as monoenergetic and unidirectional), and they found that the relaxation phase occurs much faster than the IC cooling losses. From this, they argued that more than of the beam energy can be transferred to the IGM plasma. Miniati & Elyiv (2013) reached the opposite conclusion, when accounting for the effect of the finite transverse momentum spread of blazarinduced beams. They found that the beam energy is radiated by IC off the CMB well before the relaxation phase.
With our PIC simulations, we find that the relaxation phase occurs on a much longer timescale than the exponential oblique growth, at least by two orders of magnitude. However, it might be delayed even more for more extreme beam parameters, as suggested by both Schlickeiser et al. (2012b) and Miniati & Elyiv (2013). Even under the conservative assumption that the relaxation phase is faster than the IC cooling time, this does not imply that all of the beam energy is ultimately deposited into the IGM plasma. In short, the relaxation process being faster than the IC losses does not guarantee that it will also be efficient in heating the IGM electrons. Under the unrealistic assumption that blazardriven beams are born with a small longitudinal momentum spread , we find that only (rather than , as assumed by Broderick et al. (2012)) of the beam energy is transferred to the plasma. For the realistic spread of blazarinduced beams, the coupling will be much less efficient, with the heating efficiency as low as .
5. Summary and Discussion
The interaction of TeV photons from distant blazars with the extragalactic background light produces ultrarelativistic electronpositron pairs. The resulting pair beam is unstable to the excitation of beamplasma instabilities in the unmagnetized intergalactic medium (IGM). The ultimate fate of the beam energy is uncertain, and it is hard to capture with analytical tools. By means of 2D and 3D PIC simulations, we have investigated the linear and nonlinear evolution of ultrarelativistic dilute electronpositron beams. We have performed dedicated experiments with a broad range of beam Lorentz factors and beamtoplasma density contrasts , so that our results can be extrapolated to the extreme parameters of blazarinduced beams ( and ).
We find that the earliest stages of evolution of ultrarelativistic dilute beams are governed by the oblique instability. For cold beams, the oblique mode proceeds in the socalled reactive regime, where all the beam particles are interacting with the unstable waves, and the coupling between the beam and the plasma is most efficient. In this case, the instability grows at the reactive rate , where and is the plasma frequency of the IGM electrons. However, blazarinduced beams are not cold. The initial spread in transverse velocity is , so that the oblique mode proceeds in the kinetic (rather than reactive) regime. Here, the spectrum of unstable modes is broader, with only a few beam particles being in resonance with each mode. The growth rate is , which is smaller than the corresponding rate of cold beams, yet large enough such that the kinetic oblique instability has ample time to grow, before the beam loses energy to IC emission by scattering off the CMB.
On the other hand, the oblique instability growing faster than the IC losses does not guarantee that it will also be efficient in heating the IGM electrons. Due to selfheating of the beam in the transverse direction, we find that the exponential phase of the kinetic oblique instability terminates when only a minor fraction of the beam energy has been transferred to the IGM plasma. At the end of the kinetic oblique phase, the heating efficiency of IGM electrons is extremely poor, reaching .
Additional transfer of energy from the beam to the plasma occurs at later times (with a delay of two or more orders of magnitude, relative to the oblique growth time), and it is mediated by the quasilinear relaxation process. Here, the beam generates longitudinal electrostatic waves, which scatter the beam particles – thus broadening the beam momentum spectrum in the longitudinal direction – and heat the IGM electrons. The quasilinear relaxation can operate only if the beam momentum spectrum along the longitudinal direction is sufficiently narrow. More precisely, we find that if the beam distribution peaks at , the quasilinear relaxation requires the parallel momentum spread below the peak to be . In this case, a fraction of the beam energy is transferred to the plasma (rather than , as assumed by Broderick et al. 2012).
The constraint on the longitudinal dispersion at birth is hard to fulfill by blazarinduced beams. Since the pair creation cross section peaks slightly above the threshold energy, the pairs are born moderately warm, with a comoving temperature . This corresponds to a longitudinal momentum spread at birth (measured in the IGM frame) of , which is already prohibitive for the development of the quasilinear relaxation. Moreover, the momentum dispersion of blazarinduced beams might be even larger, since the spectrum of both the EBL and the blazar TeV emission are usually modeled as broad power laws (Miniati & Elyiv, 2013). For , the quasilinear relaxation will be suppressed, which results in a much lower efficiency of IGM heating – with a firm lower limit being the heating fraction at the end of the oblique phase .
A the end of the relaxation phase, the beam and plasma distributions are highly anisotropic, with the longitudinal momentum spread much larger than the transverse one. As a result, the system is prone to the Weibel instability (e.g., Weibel, 1959; Yoon & Davidson, 1987; Silva et al., 2002), which relaxes the beam and plasma anisotropy by generating transverse magnetic fields. The Weibel instability is predominantly magnetic, so it does not mediate any further exchange of energy from the beam to the plasma electrons. Yet, it might be a promising source for the generation of magnetic fields in the IGM, as discussed by Schlickeiser et al. (2012b). A multidimensional PIC investigation of the strength and scale of the magnetic fields resulting from blazarinduced beams will be presented elsewhere.
Our results have important implications for the ultimate fate of the energy of blazarinduced beams. Since at most of the beam energy is deposited into the IGM plasma, most of the energy () is still available for IC interactions with the CMB. Therefore, the Fermi nondetection of the IC scattered GeV emission around TeV blazars can be reliably used to probe the strength of the EBL and of the IGM magnetic fields. In particular, the lower bounds on the IGM field strength derived by various authors (e.g., Neronov & Vovk, 2010; Tavecchio et al., 2010) are still valid, despite the fast growth of beamplasma instabilities in the IGM. However, when computing the expected IC signature, one should take into account that beamplasma instabilities tend to broaden the beam distribution function. As a result, the reprocessed GeV emission will be spread over a wider frequency range (and consequently, with smaller flux), as compared to the case of monoenergetic beams.
The fraction of beam energy deposited into the IGM might have important cosmological implications, as argued by Chang et al. (2012) and Pfrommer et al. (2012). By assuming that all of the beam energy is transferred to the IGM electrons by beamplasma instabilities, they showed that blazar heating could dominate over photoheating in the lowredshift evolution of the IGM, by almost one order of magnitude. If the heating efficiency of blazarinduced beams is , rather than as they assumed, blazar heating could still contribute as much as photoheating to the thermal history of the IGM. However, we expect the heating efficiency of blazarinduced beams to be , for the following reasons:

As we have argued above, blazarinduced beams are born with a significant longitudinal spread in momentum, . In this case, the quasilinear relaxation process, which mediates efficient transfer of the beam energy to the IGM electrons up to , will be suppressed, and the heating fraction remains .

Density inhomogeneities associated with cosmic structure induce loss of resonance between the beam particles and the excited plasma oscillations, strongly inhibiting the growth of the unstable modes (Miniati & Elyiv, 2013).

A largescale IGM field might spread a monoenergetic unidirectional beam in the transverse and longitudinal directions, thus inhibiting the efficiency of the beamplasma interaction.
^{15} In the oblique phase, a sufficiently strong magnetic field might trigger the transition to the kinetic regime, if during the growth time of the reactive instability () it deflects the beam velocity sideways by more than the threshold