1 Introduction

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 electron-positron pairs streaming through the intergalactic medium (IGM). The fate of the beam energy is uncertain. By means of two- and three-dimensional particle-in-cell simulations, we study the non-linear evolution of dilute ultra-relativistic pair beams propagating through the IGM. We explore a wide range of beam Lorentz factors and beam-to-plasma density ratios , so that our results can be extrapolated to the extreme parameters of blazar-induced 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 self-heating 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 quasi-longitudinal 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 blazar-induced 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 up-scattering of the Cosmic Microwave Background by the beam pairs.

Subject headings:
gamma rays: general – instabilities – intergalactic medium – plasmas – radiation mechanisms: non-thermal

1. 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 electron-positron 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 electron-positron 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 blazar-induced beams may be different. As they stream through the IGM plasma, the electron-positron 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 blazar-induced beams (i.e., dilute and ultra-relativistic), 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 beam-plasma instabilities, rather than IC cooling. In this case, the blazar TeV emission would not be reprocessed down to multi-GeV 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 non-linear 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 blazar-induced beams). The nature of the fastest growing instability can change as the beam-plasma system evolves, due to the suppression of temperature-sensitive modes as the beam heats up (e.g., Bret et al., 2008). Also, beam-plasma 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 first-principles particle-in-cell (PIC) simulations in two and three dimensions to study the non-linear stages and saturation of the instabilities generated as the blazar-induced pair beams propagate through the IGM. The non-linear effects of the beam-plasma interaction are hard to capture with analytical tools, and they require fully-kinetic simulations. We explore a wide range of beam Lorentz factors and beam-to-plasma density ratios , so that our results can be extrapolated to the extreme parameters of blazar-induced beams ( and , for the most powerful blazars). We find that, for ultra-relativistic dilute beams that start with a negligible thermal spread, electrostatic beam-plasma 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 blazar-induced 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 up-scattering 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 blazar-induced 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 ultra-relativistic cold beams) and we describe the complete evolution of the beam-plasma unstable system, from the early exponential phase up to the non-linear stages. In §4.2, we discuss the dependence of our findings on the beam Lorentz factor, the beam-to-plasma density contrast and the beam temperature. For the convenience of readers uninterested in the kinetic details of beam-plasma instabilities, in §4.3 we summarize our results in application to blazar-induced 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 Blazar-Driven Beams

In this section, we summarize the physical parameters of blazar-induced beams, including the density contrast to the IGM, the beam Lorentz factor and velocity spread. We present order-of-magnitude 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 electron-positron pairs (Neronov & Semikoz, 2009). Here, accounts for uncertainties in the intensity of the EBL, with models predicting for (e.g., Aharonian, 2001).1 Each particle moves along the direction of the incident TeV photon, and it carries about half of the photon energy, so the beam Lorentz factor is . Assuming that plasma instabilities in the IGM do not appreciably affect the beam propagation (an assumption that is correct a posteriori, as we demonstrate in this work), the pairs travel a distance of before cooling by IC scattering off the CMB (leading to many photons of energy per original TeV photon).

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 If the number density in the IGM is cm (but it may be a factor of several smaller in cosmological voids, which dominate the cosmic space at ), the density ratio between the streaming pairs and the background plasma is

(2)

The density contrast and the beam Lorentz factor are the two crucial parameters determining the plasma physics of the beam-IGM interaction. For blazar-induced beams, we expect that they should vary in the range and , respectively.

As discussed by Broderick et al. (2012), collective beam-plasma 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 skin-depth radius contains beam particles. For , we find that collective phenomena always play a role in the evolution of blazar-induced 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 non-linear evolution of blazar-induced beams.

3. Simulation Setup

We investigate blazar-driven plasma instabilities in the IGM by means of fully-kinetic PIC simulations. We employ the three-dimensional (3D) electromagnetic PIC code TRISTAN-MP (Spitkovsky, 2005), which is a parallel version of the publicly available code TRISTAN (Buneman, 1993), that was optimized for handling ultra-relativistic flows (see, e.g., Sironi & Spitkovsky 2011 and Sironi et al. 2013 for studies of ultra-relativistic collisionless shocks using TRISTAN-MP). We initialize a relativistic dilute pair beam that propagates along through an unmagnetized electron-proton 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 beam-plasma instabilities will grow from noise.

To follow the beam-plasma 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 in-plane components of the velocity, current and electric field, and only the out-of-plane 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 beam-plasma 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 non-relativistic, 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 charge-neutralizing background.

The beam consists of electron-positron pairs propagating with Lorentz factor along the direction. To our knowledge, our PIC simulations are the first to address the evolution of an electron-positron beam. All of the previous studies have focused on the case of an electron beam propagating through an electron-proton 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 beam-to-plasma density ratio and the beam Lorentz factor expected for blazar-induced 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 ultra-relativistic 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 .3 For numerical convenience, the density ratio between the beam and the plasma is established by initializing the same number of beam and plasma computational particles, with the beam particles having a weight . We have tested that, by choosing a different weight (yet, keeping the same physical density contrast), our results do not change. In addition to studying the dependence on and , we also compare the evolution of cold beams (with comoving temperature at initialization ) with the case of warm beams, up to the limit of mildly relativistic thermal spreads most relevant for blazar-induced beams.

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.4 In 2D runs, the simulation plane is typically a square with cells () on each side, but we have checked that our results do not substantially change when employing a larger box, that is long (in the direction of beam propagation) and wide. In 3D we employ a box with cells () in the transverse direction and cells () along the longitudinal direction.5 To capture the linear and non-linear stages of the beam-plasma evolution, we follow the system up to unprecedentedly long times, in 2D up to , or equivalently timesteps, and in 3D up to , or timesteps.

The number of beam particles is kept constant during the evolution of the beam-plasma system, since the photon-photon interactions that would introduce fresh electron-positron 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 beam-plasma 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, ultra-relativistic particles will emit numerical erenkov radiation. This might artificially slow down the beam, even in the absence of physical beam-plasma instabilities. We have assessed that the results reported below arise from a physical instability (as opposed to the numerical erenkov mode), by comparing our beam-plasma 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 non-linear evolution of ultra-relativistic dilute pair beams by means of 2D and 3D PIC simulations. In §4.1, we describe the different stages of evolution of the beam-plasma system, for a representative choice of beam parameters in the regime of ultra-relativistic 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 beam-to-plasma density contrast, the beam Lorentz factor and the beam temperature at birth. The reader that is not interested in the kinetic details of the beam-plasma interaction might proceed to §4.3, where we extrapolate the findings of our PIC simulations to the extreme parameters of blazar-induced beams.

Figure 1.— Temporal evolution of the beam-plasma interaction, from the 2D simulation of a cold beam with and . We follow the evolution of the system through the exponential phase of the oblique mode (panels (a)-(d)) until the relaxation stage (panels (e)-(h)). Panels (a) and (e): fraction of beam kinetic energy transferred to the plasma electrons (orange), to the longitudinal and transverse electric fields (red and blue, respectively) and to the transverse magnetic fields (green). In panel (a), the dotted orange line shows the growth rate of the oblique mode expected from linear dispersion analysis. Panels (b) and (f): temporal evolution of the momentum dispersion of the beam (solid) and plasma (dashed) electrons, along the beam (red) or transverse to the beam (blue). The momenta are in units of . Panels (c),(d) and (g): 2D plot of the longitudinal electric field , in units of . The electric field is shown at three different stages of evolution, as marked by the red arrows at the bottom of panels (a) and (e). Panel (h): 2D plot of the transverse magnetic field , in units of , at the time marked by the green arrow at the bottom of panel (e).
Figure 2.— 3D structure of the longitudinal electric field , from the 3D simulation of a cold beam with and . The electric field is in units of . The three snapshots are taken at the same times of panels (c), (d) and (g) in Figure 1, and they show that the 3D physics of the relevant electrostatic beam-plasma instabilities can be correctly captured by our 2D runs. In the 3D simulations, the temporal evolution of the fraction of beam kinetic energy transferred to the background electrons and to the electromagnetic fields (not shown here), as well as the time evolution of the beam and plasma momentum dispersions (still not shown here), closely follow the 2D results presented in Figure 1(a),(b),(e) and (f). The equivalence of 2D and 3D results for cold beams is indeed expected on analytical grounds (e.g., Bret et al., 2010b).

4.1. The Linear and Non-Linear Evolution of Ultra-Relativistic 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 non-linear relaxation. We find that the exponential phase of the oblique mode (which is the fastest growing instability for dilute ultra-relativistic beams) terminates due to self-heating 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 quasi-longitudinal 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 quasi-longitudinal 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 beam-plasma 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 ultra-relativistic dilute beams. The reactive phase of the instability governs the evolution of the system for , where is marked as a dash-dotted 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).6 The oblique mode is also captured in 3D simulations, as shown in the 3D structure of the longitudinal electric field of Figure 2(a).

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).7 From the pattern in Figure 1(c), we infer , so the growth rate of the fastest growing oblique mode will be

(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 quasi-electrostatic, 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 quasi-electrostatic, the magnetic component will be sub-dominant relative to the electric fields. In agreement with the analytical considerations of Lemoine & Pelletier (2010), we find that . Given that for ultra-relativistic 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 quasi-isotropically, 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 ultra-relativistic 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 dash-dotted 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 quasi-electrostatic (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 self-consistently 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 self-heating 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 self-excited 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 ).

Figure 3.— Temporal evolution of the longitudinal momentum spectrum , for a beam-plasma system with and . The momentum is in units of . The beam, which can be identified with the isolated peak at , stops evolving after , in agreement with Figure 1(f). The beam spectrum at late times does not relax to the so-called “plateau” distribution , which is indicated as a black dotted line. We find that the beam approaches the plateau distribution only if we artificially inhibit the evolution of the beam transverse momentum, i.e., we force the beam relaxation to proceed in a quasi-1D configuration (red dashed line). As a result of the beam relaxation, the plasma develops a high-energy tail in the forward direction. For comparison, the plasma distribution in the backward direction (i.e., opposite to the beam) is shown as a dotted red line.

The Longitudinal Relaxation Phase

After the saturation of the oblique instability, further evolution of the beam-plasma system proceeds via quasi-longitudinal 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 quasi-longitudinal, 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 quasi-longitudinal oscillations shown in Figure 1(g) are a characteristic signature of the quasi-linear relaxation of the beam (see, e.g., Grognard, 1975; Lesch & Schlickeiser, 1987; Schlickeiser et al., 2002; Pavan et al., 2011). In the quasi-linear 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 quasi-linear 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 quasi-linear 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 quasi-linear relaxation (see §4.2 for further details).

The quasi-linear 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 so-called 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 beam-to-plasma density contrast, a generic by-product of the relaxation of cold ultra-relativistic dilute beams is the conversion of of their energy into plasma heating.8

The quasi-linear relaxation significantly affects the shape of the beam and plasma longitudinal momentum spectrum, as shown in Figure 3. As a result of the quasi-longitudinal relaxation, the plasma distribution at late times () develops a pronounced high-energy tail (at ), that bridges the main thermal peak of the plasma electrons (at ) with the beam particles (that populate the isolated high-energy bump at in Figure 3). This high-energy 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 quasi-linear relaxation, the beam spectrum evolves from a quasi-monoenergetic 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 so-called “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 beam-plasma system appears stable.9 In turn, the fact that the beam spectrum is harder than explains why the amount of beam energy transferred to the plasma is only (it should be for a plateau distribution extending up to , see Thode & Sudan 1975).

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 self-consistent 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 multi-dimensional evolution of ultra-relativistic beams, but it only occurs when the beam transverse dispersion stays sufficiently small, so the system is quasi-1D.

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 beam-plasma system is achieved when the longitudinal velocity spread of the ultra-relativistic 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 beam-plasma system.

As shown in Figure 1(a), the growth of the quasi-electrostatic 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 ultra-relativistic 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).10 As a result of the Weibel instability, the magnetic field energy increases (see the green line in Figure 1(e), at ), and the beam and plasma anisotropy is reduced by increasing the transverse momentum spread (see the solid and dashed blue lines at in Figure 1(f)).

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 large-scale 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 ultra-relativistic dilute pair beams from TeV blazars will be presented elsewhere.

Figure 4.— Dependence of the beam-plasma evolution on the beam Lorentz factor (from in black up to in red, in each panel; see the legend in panels (d)-(f)) and on the beam-to-plasma density contrast ( for the leftmost column, for the middle column, and for the rightmost column). Panels (a)-(c): fraction of the beam kinetic energy transferred to the plasma electrons (in the inset, a zoom-in on the earliest phases of evolution). Panels (d)-(f): beam momentum dispersion in the longitudinal direction, normalized to the initial beam momentum (i.e., ). Panels (g)-(i): beam (solid) and total (dashed) momentum spectrum in the longitudinal direction, at the time indicated with the dotted black lines in the upper rows. In panels (g)-(i), the slope expected for the plateau distribution is shown as a dotted black line.

4.2. Dependence on the Beam Parameters

In this section, we discuss the dependence of the beam-plasma evolution on the beam parameters. In §4.2.1, we consider the case of cold beams, and we show that the quasi-longitudinal relaxation leads to a beam momentum spread along the direction of motion , regardless of the beam Lorentz factor or the beam-to-plasma 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 beam-to-plasma energy transfer. We find that if the initial dispersion in longitudinal momentum satisfies (as typically expected for blazar-induced 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 beam-plasma system depends on the beam Lorentz factor (that we vary from up to ) and on the beam-to-plasma 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 ultra-relativistic dilute electron-positron beams, as appropriate for blazar-induced beams.

For each choice of and , we follow the beam-plasma system from the oblique phase until the quasi-linear 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 self-heating 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 beam-plasma 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 long-term evolution of the system is controlled by the quasi-longitudinal relaxation, which operates on a timescale , where is the characteristic e-folding time of the oblique mode. In Figure 4, the quasi-linear 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 quasi-linear 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 quasi-linear 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 quasi-linear 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.

Figure 5.— Temporal evolution of a beam-plasma system with and , for different beam comoving temperatures at initialization (as shown by the legend in panel (b), where the beam temperature is in units of . Panel (a): fraction of the beam kinetic energy deposited into the background electrons, with the inset showing the evolution at early times. Panel (b): temporal evolution of the beam longitudinal momentum spread, in units of . Panel (c): beam (solid) and total (dashed) momentum spectra in the longitudinal direction, at the time indicated in panels (a) and (b) with the vertical black dotted line. In panel (c), the dotted oblique line shows the slope expected for a 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 .11 The growth rate in the kinetic regime is reported in Equation (9). The plasma thermal energy grows exponentially, until the self-generated electric fields increase the transverse dispersion in beam velocity beyond the initial value . At this point, the exponential growth at the rate in Equation (9) will necessarily terminate. From the Lorentz force applied to the beam particles, we find that this will happen when

(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 ).12

Regardless of the character of the oblique mode (reactive or kinetic), the long-term evolution of the beam-plasma system is controlled by the quasi-linear 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 non-relativistic. 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 quasi-linear 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 quasi-linear 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 quasi-linear phase does not mediate any further increase in the heating efficiency , beyond the early oblique phase (see the black line in Figure 5(a)).13 In short, if the initial momentum spread is , the plasma heating efficiency stays fixed at the value attained at the end of the reactive oblique phase – or , if the oblique instability proceeds in the kinetic regime.

The longitudinal momentum spectrum in Figure 5(c) clarifies why the quasi-linear 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 quasi-linear 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 quasi-linear 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 quasi-relaxation 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 quasi-linear 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 so-called 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 quasi-linear 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 power-law distribution has a negligible spread in the transverse direction (which might not be the case for blazar-induced 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 mono-energetic beams. As we have discussed above, the effectiveness of the quasi-linear 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 low-energy cutoff at , then the quasi-linear relaxation does operate, and the beam deposits of its energy into the background electrons. However, if the low-energy end of the beam distribution is broader, the quasi-linear 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 quasi-linear relaxation is suppressed if , with the case corresponding to the plateau distribution.

4.3. Implications for Blazar-Induced Beams

We now analyze the implications of our findings for the evolution of blazar-induced beams in the IGM. As we have discussed in §3, the Lorentz factors and density contrasts of blazar-driven 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 ultra-relativistic dilute beams depends on the beam Lorentz factor and the beam-to-plasma density ratio. Our results can then be confidently extrapolated to the extreme parameters of blazar-induced beams.

We have demonstrated that the oblique instability, which governs the earliest stages of evolution of ultra-relativistic 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 photon-photon 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 blazar-induced beams (). The oblique instability will grow at the kinetic rate in Equation (9), which for , as appropriate for blazar-induced beams, reduces to . For the parameters of blazar-induced 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.14

Due to self-heating 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 blazar-induced 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 quasi-linear relaxation phase. We find that the quasi-linear 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 quasi-linear 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 quasi-linear relaxation (see Schlickeiser et al., 2012b).

In §4.2.2, we have demonstrated that the quasi-linear 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 power-law tail, this would be the low-energy cutoff), the quasi-linear relaxation can operate only if the parallel momentum spread below the peak satisfies . This constraint is hard to fulfill by blazar-induced 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 quasi-linear relaxation. Moreover, the momentum dispersion of blazar-induced 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 quasi-linear 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 quasi-linear relaxation phase, the beam spectrum will be harder than the so-called 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 quasi-linear 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 blazar-induced beams in the IGM. As we have emphasized in §3, our PIC simulations are the first to address the evolution of dilute ultra-relativistic electron-positron beams, as appropriate for blazar-induced 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 blazar-induced 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 blazar-induced beams is controlled by the oblique (rather than longitudinal) mode, transverse thermal effects are indeed important for the beam-plasma interaction at early times.

The evolution of blazar-induced beams at late times is governed by non-linear plasma processes, which are extremely hard to capture with analytical tools. Schlickeiser et al. (2012b) and Miniati & Elyiv (2013) attempted to describe the non-linear relaxation of blazar-induced 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 mono-energetic and uni-directional), 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 blazar-induced 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 blazar-driven 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 blazar-induced 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 ultra-relativistic electron-positron pairs. The resulting pair beam is unstable to the excitation of beam-plasma 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 non-linear evolution of ultra-relativistic dilute electron-positron beams. We have performed dedicated experiments with a broad range of beam Lorentz factors and beam-to-plasma density contrasts , so that our results can be extrapolated to the extreme parameters of blazar-induced beams ( and ).

We find that the earliest stages of evolution of ultra-relativistic dilute beams are governed by the oblique instability. For cold beams, the oblique mode proceeds in the so-called 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, blazar-induced 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 self-heating 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 quasi-linear 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 quasi-linear 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 quasi-linear 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 blazar-induced 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 quasi-linear relaxation. Moreover, the momentum dispersion of blazar-induced 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 quasi-linear 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 multi-dimensional PIC investigation of the strength and scale of the magnetic fields resulting from blazar-induced beams will be presented elsewhere.

Our results have important implications for the ultimate fate of the energy of blazar-induced 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 non-detection 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 beam-plasma instabilities in the IGM. However, when computing the expected IC signature, one should take into account that beam-plasma 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 mono-energetic 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 beam-plasma instabilities, they showed that blazar heating could dominate over photo-heating in the low-redshift evolution of the IGM, by almost one order of magnitude. If the heating efficiency of blazar-induced beams is , rather than as they assumed, blazar heating could still contribute as much as photo-heating to the thermal history of the IGM. However, we expect the heating efficiency of blazar-induced beams to be , for the following reasons:

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

  • Even if the relaxation process were to be operating, the beam relaxation timescale might be much longer than the IC loss time, as argued by Miniati & Elyiv (2013) (but see Schlickeiser et al. (2012b), for opposite conclusions).

  • 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 large-scale IGM field might spread a mono-energetic uni-directional beam in the transverse and longitudinal directions, thus inhibiting the efficiency of the beam-plasma 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