THE COSMOLOGICAL IMPACT OF LUMINOUS TeV BLAZARS. I. IMPLICATIONS OF PLASMA INSTABILITIES FOR THE INTERGALACTIC MAGNETIC FIELD AND EXTRAGALACTIC GAMMA-RAY BACKGROUND

Sources of VHEGRs are necessarily attenuated by the production of pairs upon interacting with the EBL. Namely, when the energies of the gamma ray (E) and the EBL photon (E_EBL) exceed the rest mass energy of the e^± pair in the center-of-mass frame, i.e., 2EE_EBL(1 − cos θ) > 4m²_ec⁴, where θ is the relative angle of propagation in the lab frame, an e^± pair can be produced with Lorentz factor γ ≃ E/2m_ec² (Gould & Schréder 1967). The optical depth the universe presents to high-energy gamma rays depends solely on the energy of the gamma ray and the evolving spectrum of the EBL.⁹ Detailed estimates for the EBL spectrum and its evolution have produced an estimated mean free path of TeV photons of

$$\begin{equation} D_{\rm pp}(E,z) = 35\left(\frac{E}{1\,{\rm T}{\rm eV}}\right)^{-1} \left(\frac{1+z}{2}\right)^{-\zeta }\,{\rm M}{\rm pc}\,, \end{equation} \tag{ 1 }$$

where the redshift evolution is due to that of the EBL alone, is dependent predominately on the star formation history, and ζ = 4.5 for z < 1 and ζ = 0 for z ⩾ 1 (Kneiske et al. 2004; Neronov & Semikoz 2009).¹⁰ The resulting optical depth to VHEGRs emitted with energy E from an object located at z is

$$\begin{equation} \tau _E(E,z) \equiv \int _0^z \frac{c dz^{\prime }}{D_{\rm pp}\left[E(1+z^{\prime })/(1+z),z^{\prime }\right] H(z^{\prime })(1+z^{\prime })} \propto E\, \end{equation} \tag{ 2 }$$

(in which H(z) is the Hubble factor), shown in the top panel of Figure 1.¹¹ From this it is evident that above 100 GeV the universe is optically thick to sources at z > 1 (cf. Franceschini et al. 2008). A related and perhaps more appropriate measure of the optical depth is that across a Hubble length, τ_H ≡ c/(D_ppH), shown in the bottom panel of Figure 1, providing a sense of how opaque the universe is as a function of redshift. In all cases it is clear that at E ≳ 100 GeV photons will pair produce on the EBL for z ≲ 10. In the future this may not be the case, since for sufficiently small z, τ_H decreases with decreasing z due to the dramatic decrease in the EBL photon density associated with both the Hubble expansion and the slowing of star formation.

Zoom In Zoom Out Reset image size

The pairs produced by the TeV gamma rays are ultrarelativistic, with typical Lorentz factors of E/2m_ec² ≃ 10⁶(E/TeV). They also necessarily constitute a cold, highly anisotropic, dilute beam propagating through the IGM. This follows immediately from the intergalactic distances traversed by the gamma rays and the comparatively small EBL photon energies (as seen in the IGM frame). Here we estimate the properties of this plasma beam.

The momentum dispersion of the resulting pairs is set by the gamma-ray spectrum, geometry of the TeV source, distribution of the EBL photons, and heating due to pair production. Of these, only the last plays a significant role in setting the transverse momentum dispersion.¹² The center-of-mass frame, i.e., the "beam frame," momentum dispersion resulting from pair production is roughly m_ec. This results in an IGM-frame transverse momentum dispersion of p_⊥/p_∥ ≃ 10⁻⁶(E/TeV)⁻¹. With p_⊥, m_ec ≪ p_∥, the temperature¹³ associated with this transverse momentum dispersion is

$$\begin{equation} \frac{k T_b}{m_e c^2} \simeq \frac{p_\parallel }{2 m_e c^2} \left(\frac{p_\perp }{p_\parallel }\right)^2 \simeq 5\times 10^{-7}\left(\frac{E}{{\rm T}{\rm eV}}\right)^{-1}\,. \end{equation} \tag{ 3 }$$

Since the transverse momentum dispersion of the pairs is much smaller than that associated with the bulk motion of the beam (i.e., since p_⊥/p_∥ ≪ 1), we may safely assume that the beam is transversely kinematically cold.

The density of the pair beam at a given point within the IGM is set by the rate at which pairs are produced, duration of the TeV emission, advection of the pairs through the IGM, and the processes by which they lose their kinetic energy. That is, the evolution of the density of pairs per unit Lorentz factor, n_γ, is governed by the Boltzmann equation:

$$\begin{equation} \frac{\partial n_\gamma }{\partial t} + \frac{c}{r^2} \frac{\partial r^2 n_\gamma }{\partial r} + \dot{\gamma } \frac{\partial n_\gamma }{\partial \gamma } = \dot{n}_\gamma \,, \end{equation} \tag{ 4 }$$

where the left-hand side assumes that all the pairs are moving away from the TeV source relativistically (v^r = c and p^r = γm_ec), and the right-hand side corresponds to pair production. Generally, we may neglect advection, which alters n_γ over timescales of D_pp/c, much longer than any relevant timescale of interest here (i.e., (c/r²)∂r²n_γ/∂r may be neglected). Furthermore, for most of the potential sources we will consider (primarily TeV blazars), we will assume that the duration of the TeV emission is sufficiently long that n_γ reaches a steady state (i.e., ∂n_γ/∂t = 0). In this case, we have $\dot{\gamma } \partial n_\gamma /\partial \gamma \simeq \dot{n}_\gamma$.

Making further progress requires us to define the spectrum of the pairs, which itself depends on the spectrum of the gamma rays and the energy dependence of the cooling processes. Nevertheless, we may obtain an estimate of the beam density in the vicinity of a given Lorentz factor, n_b ≃ γn_γ, by setting ∂n_γ/∂γ ≃ −n_γ/γ = −n_b/γ².

The source term is given by twice (since each gamma ray produces two leptons) the rate at which high-energy gamma rays with energy E ≃ 2γm_ec² annihilate within the region of interest, i.e., $\gamma \dot{n}_\gamma = 2 (E dN/dE)/D_{\rm pp}= 2 F_E/D_{\rm pp}$, where N is the gamma-ray number flux, with units of photons cm⁻² s⁻¹. Thus, upon defining a cooling rate $\Gamma \equiv -\dot{\gamma }/\gamma$, we have

$$\begin{equation} n_{b}\simeq \frac{2 F_E}{D_{\rm pp}\Gamma }\,, \end{equation} \tag{ 5 }$$

i.e., the density of pairs of a given energy is determined by balancing cooling and pair creation.

Generally, Γ is a function of energy and beam density, as well as external factors (e.g., seed photon density, IGM density, etc.). Thus, this gives a nonlinear algebraic equation to solve for n_b, the particulars of which depend on the various mechanisms responsible for extracting the bulk energy of the beam. In practice, given expressions for Γ, associated with the processes discussed in the following section, we solve Equation (5) numerically to obtain n_b(E, F_E, z).

Which mechanism dominates the cooling of the beam depends on a variety of environmental factors and the properties of the pair beam itself. Nevertheless, inverse Compton cooling via the CMB provides a convenient lower limit on Γ and thus an upper limit on n_b. This is a function of z and γ alone, given by

$$\begin{equation} \Gamma _{\rm IC}= \frac{4\sigma _{\rm T} u_{\rm CMB}}{3 m_e c} \gamma \simeq 1.4\times 10^{-20}(1+z)^4\gamma \,\,{\rm s}^{-1}\,, \end{equation} \tag{ 6 }$$

where σ_T denotes the Thompson cross section. The strong redshift dependence arises from the rapid increase in the CMB energy density with z (u_CMB∝(1 + z)⁴). Furthermore, since it is ∝γ, inverse Compton cooling is substantially more efficient at higher energies. The associated cooling rate is shown as a function of E for a number of redshifts in Figure 2.

Zoom In Zoom Out Reset image size

When we set Γ ≃ Γ_IC, we obtain the following upper limit on the beam density:

$$\begin{eqnarray} n_{b}&\simeq & \frac{2 F_E}{D_{\rm pp}\Gamma _{\rm IC}} \simeq \frac{L_E}{2\pi D_{\rm pp}^3\Gamma _{\rm IC}}\simeq 3.7\times 10^{-22} \left(\frac{1+z}{2}\right)^{3\zeta -4} \nonumber\\ &&\times\left(\frac{E L_E}{10^{45}\ {\rm erg}\,{\rm s}^{-1}}\right) \left(\frac{E}{{\rm T}{\rm eV}}\right) \,{\rm cm}^{-3}\,, \end{eqnarray} \tag{ 7 }$$

where we have defined L_E to be the isotropic-equivalent luminosity (per unit energy) of a source located at a distance D_pp from the region in question. Setting EL_E to a typical value (10⁴⁵ erg s⁻¹) gives an idea of the typical pair beam densities. Note that despite the large blazar luminosities we consider, the associated beams are exceedingly dilute, a point that is of critical importance in the following section. Since Γ_IC is independent of n_b, this has no implication for inverse Compton cooling itself.

For collective phenomena to be relevant, it is necessary for many pairs to be present within each wavelength of the unstably growing modes. As we shall see, the relevant scale for the beam plasma instabilities we describe below is the plasma skin depth of the IGM,

$$\begin{equation} \lambda _P \equiv \frac{2\pi c}{\omega _P} = \sqrt{\frac{\pi m_e c^2}{e^2 n_{\rm IGM}}} = 2.3\times 10^{-9}(1+\delta)^{-1/2}(1+z)^{-3/2}\,{\rm pc}\,, \end{equation} \tag{ 8 }$$

where $\omega _P\equiv \sqrt{4\pi e^2n_{\rm IGM}/m_e}$ is the IGM plasma frequency and n_IGM = 2.2 × 10⁻⁷(1 + δ)(1 + z)³ cm⁻³ is the IGM free-electron number density. Generally, the growing mode must be uniform on scales considerably larger than λ_P, both longitudinally and transversely (otherwise, it is not well represented as a single Fourier component), and thus the volume in which many particles must be present is much larger than that defined by a sphere of diameter λ_P. Nevertheless, this gives us a conservative constraint, i.e., we require

$$\begin{equation} \frac{\pi }{6} \lambda _P^3 \,n_{b}\gg 1\,. \end{equation} \tag{ 9 }$$

With Equation (5) this gives a maximum plasma cooling rate:

$$\begin{equation} \Gamma _{\rm plasma} \ll \frac{\pi }{3} \frac{F_E}{D_{\rm pp}}\lambda _P^3\,. \end{equation} \tag{ 10 }$$

Plasma processes with cooling rates that exceed this limit necessarily saturate near this cooling rate. Such a super-critical process can potentially operate only until the beam density is driven below the value at which pairs can support collective phenomena, at which point they necessarily quench. However, after plasma cooling ceases, the plasma beam density rises again (since Γ_IC is always less than Γ_plasma in practice), and thus the super-critical plasma cooling may resume. Since the efficiency of the plasma cooling decreases smoothly to zero at the critical density, this sequence stabilizes near Γ_plasma. The associated excluded region is shown by the gray region in the upper left corner of Figure 2.

The upper limit on n_b obtained in Equation (7) implies an analogous limit on the relevant source isotropic-equivalent luminosities. Again, this arises from requiring a sufficient number of beam pairs to be present to support plasma processes. Since the beam density is linearly dependent upon the source luminosity, this gives a constraint on the latter:

$$\begin{eqnarray} E L_E &\gg & 12\,E \frac{D_{\rm pp}^3\Gamma _{\rm IC}}{\lambda _P^{3}} = 3.3\times 10^{38} (1+\delta)^{3/2}\nonumber\\ &&\times\left(\frac{1+z}{2}\right)^{8.5-3\zeta } \left(\frac{E}{{\rm T}{\rm eV}}\right)^{-1}\,, \end{eqnarray} \tag{ 11 }$$

effectively defining a luminosity cutoff, below which collective plasma phenomena may be ignored. This limiting luminosity is shown as a function redshift by the gray lines in Figure 3 for various gamma-ray energies. While the luminosity limit does depend on z, it is clear that all of the observed TeV blazars (listed in Table 1) are sufficiently luminous for their resulting pair beams to support collective phenomena at the redshifts of interest for active galactic nuclei (AGNs). The only source to fall marginally below the limiting luminosity at 1 TeV is the radio galaxy Cen A (the lowest point in Figure 3), the closest and dimmest object in Table 1.

Zoom In Zoom Out Reset image size

Within astrophysical contexts there are at least two well-known beam instabilities: two-stream and Weibel, the latter having been suggested as a mechanism for magnetizing strong shocks (Medvedev & Loeb 1999), and both implicated in the coupling at collisionless shocks (Spitkovsky 2008). These are, however, simply different limiting examples of the same underlying filamentary instability, for which the maximum growth rate exceeds either—the so-called oblique mode, which we discuss below (Bret et al. 2004, 2005a; Lemoine & Pelletier 2010). Note that since the pair beam is neutral, it contains its own return current, and thus beam instabilities that arise due to the electron return currents within the background plasma (e.g., the Bell and Buneman instabilities; see Bret 2009) are not relevant.

Here, we discuss the nature of these instabilities and their growth rates in the cold-plasma limit (i.e., mono-energetic beams). While the low-temperature approximation is unlikely to be applicable in practice for the beams of interest here, it provides a convenient limit in which to present the relevant processes within their broader context. For a similar reason, and because we have not found analogous derivations elsewhere in the astrophysical literature, we present the ultrarelativistic pair beam instability growth rates for the two-stream and Weibel instabilities in Appendixes A.1 and A.2, only summarizing the results here. We defer a discussion of the more directly relevant warm-plasma oblique instability to the following section. Generally, we find that the plasma instabilities are capable of dominating inverse Compton cooling as a means to dissipate the bulk kinetic energy of the pair beams from TeV blazars.

The pair two-stream instability arises due to the interaction of the anisotropic electron and positron distribution functions with the comoving background electrostatic wave (i.e., $\mbox{\boldmath $\rm k$}\parallel \mbox{\boldmath $\rm p$}$, where $\mbox{\boldmath $\rm k$}$ is the electromagnetic mode wavevector and $\mbox{\boldmath $\rm p$}$ is the beam momentum) with wavelength ∼λ_P for n_b/n_IGM ≪ 1, which is generally the case of interest here. The associated cooling rate¹⁴ in the cold-plasma limit is

$$\begin{equation} \Gamma _{\rm TS}= \sqrt{\frac{12\pi e^2 n_{\rm IGM}}{m_e}}\left(\frac{n_{b}}{2n_{\rm IGM}}\right)^{1/3}\frac{1}{\gamma }\,, \end{equation} \tag{ 12 }$$

which depends only weakly on the IGM density and the pair beam density, though it decreases rapidly as γ becomes large.

The Weibel instability is associated with coupling to a secularly growing, anharmonic transverse magnetic perturbation (i.e., $\mbox{\boldmath $\rm k$}\perp \mbox{\boldmath $\rm p$}$). The most rapidly growing wavelength is again that associated with the background plasma skin depth, ∼λ_P, with associated cooling rate in the cold-plasma limit

$$\begin{equation} \Gamma _{\rm W}= \sqrt{\frac{16\pi e^2 n_{b}}{m_e \gamma }}\,, \end{equation} \tag{ 13 }$$

which depends only on the pair beam density and the pair Lorentz factor. At large γ this is suppressed more weakly than the two-stream instability, though it is a moderately stronger function of n_b.

The computation of the growth rates of the two-stream and Weibel instabilities is greatly simplified by the particular geometries of the coupled electromagnetic waves and beam momenta. However, a generalized oblique treatment has shown the presence of continuum of unstable modes (Bret et al. 2004, 2005a, 2010b; Bret 2009; Lemoine & Pelletier 2010), characterized by the orientation of their wavevector relative to the bulk beam velocity. Of these, neither the two-stream nor Weibel is generally the most unstable. Rather, generically the most robust and fastest-growing mode occurs at oblique wavevectors and is thus referred to as the oblique instability by Bret et al. (2010b). The cold-plasma cooling rate of this maximally growing mode is

$$\begin{equation} \Gamma _{\rm M}\simeq \sqrt{\frac{12\pi e^2n_{\rm IGM}}{m_e}}\left(\frac{n_{b}}{2\gamma n_{\rm IGM}}\right)^{1/3} = \Gamma _{\rm TS}\gamma ^{2/3} \,. \end{equation} \tag{ 14 }$$

This can be much larger than the two-stream and Weibel growth rates when n_b/n_IGM ≪ 1 and γ ≫ 1.

The cold-plasma oblique instability cooling rates dominate inverse Compton cooling by orders of magnitude over the region of interest. However, at the very dilute beam densities of relevance here, the cold-plasma approximation requires exceedingly small beam temperatures. Above an IGM-frame temperature of roughly

$$\begin{eqnarray} \frac{kT_{\rm crit}}{m_e c^2} &\simeq & \frac{3}{2^{10/3}} \left(\frac{n_{b}}{n_{\rm IGM}}\right)^{2/3} \gamma ^{1/3}\nonumber\\ &\simeq & 1.0\times 10^{-9} \left(1+\delta \right)^{-2/3} \left(\frac{1+z}{2}\right)^{2\zeta -14/3}\nonumber\\ &&\times \left(\frac{E L_E}{10^{45}\ {\rm erg}\,{\rm s}^{-1}}\right)^{2/3} \left(\frac{E}{{\rm T}{\rm eV}}\right)\,, \end{eqnarray} \tag{ 15 }$$

pairs can traverse many wavelengths of the unstable modes over the cold-instability growth timescale, a situation commonly referred to as the kinetic regime. As a consequence, significant phase mixing can occur, substantially reducing the effective growth rate (Bret et al. 2010a).

The way in which finite beam temperatures limit the growth rate depends sensitively on the nature of the velocity dispersion and the modes of interest (see, e.g., Bret et al. 2005b). For example, the two-stream instability, associated with wavevectors parallel to the beam, enters the strongly suppressed "quasi-linear" regime when the parallel momentum dispersion is large, though it is insensitive to even large transverse dispersions. Conversely, the Weibel instability, associated with wavevectors orthogonal to the beam, is sensitive to even small transverse velocity dispersions but unaffected by large parallel velocity dispersions. This is simply because a given mode can tolerate large velocity dispersions within but not across the phase fronts of the unstable electromagnetic modes. For the situation of interest here, dilute beams and cool IGM (i.e., kT_IGM/m_ec² ≪ 1), the oblique modes are nearly transverse and thus sensitive primarily to large transverse velocity dispersions.

Nevertheless, even for the small temperatures we have inferred for the pair beams, we find ourselves in the kinetic regime. In this case the oblique instability cooling rate has been numerically measured to be

$$\begin{equation} \Gamma _{\rm M,k}\simeq 0.4 \frac{m_e c^2}{k T_b} \frac{n_{b}}{n_{\rm IGM}} \omega _P \simeq 0.4 \gamma \sqrt{\frac{4\pi e^2 n_{b}^2}{m_e n_{\rm IGM}}}\,, \end{equation} \tag{ 16 }$$

where we have set the beam temperature in the beam frame to m_ec²/k, and thus kT_b = m_ec²/γ (Bret et al. 2010a). Both the cold and hot growth rates have been verified explicitly using particle-in-cell (PIC) simulations, though for somewhat less dilute beams than those we discuss here (Bret et al. 2010b).

When the kinetic oblique mode dominates the beam cooling, the beam density is given by

$$\begin{eqnarray} n_{b} &\simeq & 2.1\times 10^{-24} \left(1+\delta \right)^{1/4} \left(\frac{1+z}{2}\right)^{(6\zeta +3)/4}\nonumber\\ &&\times \left(\frac{E L_E}{10^{45}\,{\rm erg}\,{\rm s}^{-1}}\right)^{1/2} \left(\frac{E}{{\rm T}{\rm eV}}\right)^{1/2} {\rm c}{\rm m}^{-3}\,. \end{eqnarray} \tag{ 17 }$$

The associated cooling rate is then

$$\begin{eqnarray} \Gamma _{\rm M,k} &\simeq & 3.6\times 10^{-11} \left(1+\delta \right)^{-1/4} \left(\frac{1+z}{2}\right)^{(6\zeta -3)/4}\nonumber\\ &&\times \left(\frac{E L_E}{10^{45}\,{\rm erg}\,{\rm s}^{-1}}\right)^{1/2} \left(\frac{E}{{\rm T}{\rm eV}}\right)^{3/2} {\rm s}^{-1}\,. \end{eqnarray} \tag{ 18 }$$

This is a stronger function of gamma-ray energy than inverse Compton cooling, implying that it will eventually dominate at sufficiently high energies, assuming a flat TeV spectrum. In addition, it is a very weak function of δ, being only marginally faster in lower-density regions, and thus the cooling of the pairs is largely independent of the properties of the background IGM.

The rates obtained by numerically solving Equation (5) for n_b, with Γ = Γ_IC + Γ_{M, k}, are shown for a number of redshifts in Figure 2. For the luminosity shown (EL_E = 10⁴⁵ erg s⁻¹, typical of the bright TeV blazars) at z ≲ 4, plasma cooling dominates inverse Compton above a TeV. In the present epoch, Γ_{M, k} is roughly two orders of magnitude larger than Γ_IC for bright TeV blazars.

The luminosity dependence of Γ_{M, k} implies a luminosity limit below which inverse Compton does dominate the linear evolution of the pair beam,

$$\begin{equation} E L_E \gtrsim 7.4\times 10^{40}\, (1+\delta)^{1/2} \left(\frac{1+z}{2}\right)^{9.5-3\zeta } \left(\frac{E}{{\rm T}{\rm eV}}\right)^{-1} \end{equation} \tag{ 19 }$$

(note that at this luminosity Γ = Γ_IC + Γ_{M, k} = 2Γ_{M, k}, and thus n_b is half the value shown in Equation (17)). This limit is shown as a function of redshift for a number of different energies by the black lines in Figure 3. Note that these lie above the corresponding limits associated with the applicability of the plasma prescription, suggesting that in practice the beam does support collective phenomena. The critical luminosity ranges from ∼10⁴⁰ erg s⁻¹ to 10⁴⁵ erg s⁻¹, depending on redshift and energy of interest. It also depends on the overdensity as roughly ∝(1 + δ)^1/2, and thus at low densities the critical luminosity is moderately smaller. The only two sources in Table 1 that fall below the plasma cooling luminosity limit are the radio galaxies M87 and Cen A, both of which are detectable only as a result of their close proximity. At 1 TeV, all but a handful of the remaining sources lie more than an order of magnitude above this limit, and in the case of the two sources that dominate the TeV flux at Earth, more than two orders of magnitude above this limit.

Given the technical nature of our discussion above, it is useful to have a qualitative understanding of these instabilities. We caution, however, that intuitive pictures of plasma processes frequently fail to capture all the relevant physics. Hence, generalization of these intuitive pictures beyond their limited range of applicability is potentially misleading. This is explicitly illustrated by our examples here: all of the instabilities discussed in this work belong to the same family, i.e., instabilities that arise from interpenetrating plasmas, but the underlying qualitative pictures for each differ substantially.

We begin with the Weibel (or filamentation) instability (Weibel 1959), for which a mechanical viewpoint (an initial magnetic perturbation deflects particles into opposing current streams that reinforce the perturbed field) can be found in Medvedev & Loeb (1999), to which we refer the interested reader. Here, we present a picture that, although having the virtue of being simpler, is not entirely correct: because like currents attract, small-scale current perturbations arising out of the fluctuations within the interpenetrating plasmas will coalesce preferentially to produce increasingly larger-scale currents. These induce stronger magnetic fields and thus larger attractive Lorentz forces between neighboring currents; a positive feedback loop develops, leading to instability. This process continues until the associated magnetic field strengths become sufficiently large to disrupt the currents (in the case of equal-density beams) or until the transverse velocity of the constituent particles is large enough to efficiently migrate between current structures on the linear growth timescale, i.e., enter the kinetic regime.

We should note that the aggregation of currents in the Weibel instability does not rely on oscillatory waves.¹⁵ Hence, there is no oscillatory component to this instability—it is a purely growing mode. In contrast, the two-stream instability is an overstable mode, where the oscillatory components are Langmuir (or plasma) waves with a wavevector parallel to the beam velocity. These are longitudinal waves, associated with local charge oscillations, and are completely described by a propagating perturbation in the electrical potential. As particles in the beam traverse Langmuir waves in the background plasma, they experience successive periods of acceleration and deceleration, with electrons (positrons) collecting in minima (maxima) of the electric potential, where the particle speeds are at their smallest.¹⁶ This charge-separated bunching of the beam plasma enhances the background electric perturbation, potentially growing the background Langmuir wave.

The bunching within the beam is simply an excitation of Langmuir waves within the beam plasma itself. Thus, the growth of the charge perturbations in the background and beam plasmas corresponds to the resonant coupling between Langmuir waves in the background and beam. When the beam density is much less than the background density, as is the case here, background and beam Langmuir waves only overlap in frequency, and therefore satisfy the conditions for resonance, when the wavevector of the latter is parallel to the beam velocity (see the discussion above Equation (A9)). This is always satisfied for the family of comoving beam Langmuir waves (i.e., waves which in the beam frame moves counter to the background plasma). However, of particular importance for the two-stream instability are the counterpropagating beam Langmuir waves (i.e., waves which in the beam frame moves parallel to the background plasma). If the beam velocity exceeds the phase velocity of these waves (as seen in the beam frame), the counterpropagating wave will be dragged in the direction of the beam (as seen in the background frame) and therefore has a wavevector that satisfies the resonant condition.

Nevertheless, the counterpropagating Langmuir wave still carries momentum in the direction opposite to the beam. Thus, as the counterpropagating wave grows, the momentum, and therefore energy, of the beam-wave system necessarily decreases. This implies that the counterpropagating Langmuir wave is also a negative-energy mode as seen in the background frame. As a consequence, the resonant coupling can transfer energy to the positive-energy background wave from the negative-energy beam wave, while growing the amplitudes of both, thereby leading to instability. Note that if the phase velocity of the counterpropagating Langmuir wave is larger than the beam velocity, it no longer is dragged in the direction of the beam and no longer satisfies the necessary resonant condition. For a distribution of particles, this constraint on the velocities within the beam corresponds to the familiar Penrose criterion (Sturrock 1994; Boyd & Sanderson 2003) and is satisfied in the pair beams resulting from VHEGRs.

The oblique instability encompasses the two-stream and Weibel instabilities, though the most unstable mode is most similar to the former in that the intuitive picture focuses solely on the electrostatic forces, ignoring electromagnetic forces.¹⁷ In practice, this is a relatively good approximation and can be used to calculate the growth of these modes in idealized situations (e.g., Nakar et al. 2011; Bret et al. 2004). The qualitative picture proceeds similarly to that for the two-stream instability described above, with the minor modification that the perturbing background Langmuir waves now move at an angle θ relative to the relativistic beam. As a consequence, the resonant Langmuir waves have a phase velocity such that v_k ≈ ccos θ. From the simple intuitive picture above, if θ is selected such that the projected beam velocity is slightly faster than the nearly resonant Langmuir wave, an instability develops.

Finally, to understand why the growth rates between the two-stream instability above and the oblique instability differ, it is useful to make the approximation that the electric fields generated are in the direction of the k-vector. In the two-stream case, the electric field must slow down (or speed up) particles. This gets progressively harder for more relativistic particles. In the oblique case, the electric field deflects the particles, changing their projected velocity. While this is also more difficult for more relativistic particles, this is not nearly as hard as changing the particles' parallel (along the beam) velocity. Hence, the oblique instability more easily drives charge density enhancements (and therefore instabilities) at large θ, i.e., easier deflection, than the two-stream instability.

We have thus far only treated the linear development of the relativistic two-stream, Weibel, and oblique instabilities. However, the impact pair beams have on the IGM, gamma-ray cascade emission, and measures of the IGMF will ultimately depend on their nonlinear development. To address this, however, we are currently forced to appeal to analytical and numerical calculations of systems in somewhat different (and less extreme) parameter regimes.

Motivated by the applicability of the Weibel instability in the context of gamma-ray bursts (GRBs), the nonlinear saturation of the relativistic Weibel instability for equal-density plasma beams is well understood analytically and numerically. Initially, the Weibel instability rapidly grows until the energy density of the generated magnetic field becomes of order the kinetic energy of the two beams (Silva et al. 2003; Frederiksen et al. 2004; Chang et al. 2008). Analytically, Davidson et al. (1972) argued that the Weibel instability would saturate when the generated magnetic fields become so large that the Larmor radius of the beam particles is of order the skin depth, i.e., when the energy of generated magnetic fields is equal to the kinetic energy of two equal-density relativistic beams (see also Medvedev & Loeb 1999). The particles rapidly isotropize with a Maxwellian distribution (Spitkovsky 2008), i.e., heat, and the magnetic energy then rapidly decays within an order of a few tens of skin depths (Chang et al. 2008). Hence, for two equal-density, relativistic, interpenetrating beams, the Weibel instability converts anisotropic kinetic energy into heat. However, as we have already noted, for the pair beams of interest here the Weibel instability is completely suppressed for tiny transverse beam temperatures. Hence, while this instability may initially operate, it may quickly become suppressed if it results in significant transverse heating of the beam.

Unlike the Weibel instability, the two-stream and kinetic oblique instabilities continue to operate, though more slowly than implied by their cold-plasma limits, in the presence of substantial beam temperatures. Due to the geometry of the coupled modes, the oblique instability is primarily sensitive to transverse beam-velocity dispersions, though it shares the resistance to beam heating with its two-stream cousin (Bret et al. 2005b, 2010b; Bret 2009; Lemoine & Pelletier 2010). Unlike the two-stream instability, the oblique instability in the parameter range of interest here is also largely insensitive to longitudinal heating. Thus, generally, it appears that the oblique instability is substantially more robust than its more commonly discussed brethren.

What is less clear a priori is if these instabilities primarily heat the beam or primarily heat the background plasma. Here, we appeal to the numerical simulations of Bret et al. (2010b), where a mildly relativistic beam (γ = 3) penetrating into a hot, dense background plasma (beam-to-background density ratio of 0.1) was studied. In these the oblique instability resulted in a significant fraction (∼20%) of the beam energy heating the background plasma before the heating of the beam suppressed the oblique instability in favor of the two-stream instability. The relative effectiveness with which the beam heats the background plasma is due to the efficiency with which the longitudinal electrostatic modes are dissipated in the background plasma; unlike electromagnetic modes (e.g., those generated by the Weibel instability), electrostatic modes are rapidly dissipated via Landau damping.

We note that Bret et al. (2010b) found that the heating of the beam by the oblique instability eventually led to its suppression, allowing the two-stream instability to grow, continuing the dissipation of the beam kinetic energy. As a result, in their simulations a total of 30% of the energy was deposited into the background plasma via a combination of the oblique and two-stream instabilities, i.e., an additional 10% of the beam energy was thermalized via the two-stream instability during and after the suppression of the oblique instability.

In our case, we expect that much more beam energy (more than 20%, and possibly up to ∼100%) will be deposited into the background IGM because we are much deeper into the regime in which the oblique instability dominates, i.e., γ ≫ 1 and n_b/n_IGM ≪ 1. In particular, for the case of interest here, γ ∼ 10⁶ and

$$\begin{eqnarray} \frac{n_{b}}{n_{\rm IGM}} &\simeq & 1.3\times 10^{-18} \left(1+\delta \right)^{-3/4} \left(\frac{1+z}{2}\right)^{(6\zeta -9)/4} \nonumber\\ &&\times \left(\frac{E L_E}{10^{45}\,{\rm erg}\,{\rm s}^{-1}}\right)^{1/2} \left(\frac{E}{{\rm T}{\rm eV}}\right)^{1/2} {\rm c}{\rm m}^{-3}. \end{eqnarray} \tag{ 20 }$$

The extremely large Lorentz factor and tiny density ratio make it computationally prohibitive to assess the beam evolution with numerical PIC directly. Nevertheless, because we find ourselves in a regime in which the oblique instability is much more strongly dominant than that simulated in Bret et al. (2010b), we expect the linear growth of the kinetic oblique instability to continue for much longer before the beam changes character and moves out of the oblique-dominated regime—potentially once ∼100% of the beam energy has been dissipated. However, this remains to be studied in future work.

The effect of nonlinear processes might also affect the evolution of the linear instability. For instance, Lesch & Schlickeiser (1987) argued that for the relativistic electrostatic two-stream instability, nonlinear coupling to daughter modes arrests the growth of the linearly unstable mode at a very low mode energy. Hence, they claim that the electrostatic two-stream instability can only bleed energy from the beam at a slow rate.

Utilizing the order-of-magnitude estimates in Lesch & Schlickeiser (1987) or the expression for nonlinear Landau damping in Melrose (1980), we have found that the damping of the pair beam via the relativistic two-stream instability could be highly suppressed. The oblique instability may be similarly suppressed; however, this effect is much more marginal due to the instability's much larger growth rate. Nevertheless, determining its behavior for the parameters relevant here is an important unanswered question that is left for future work. For the purposes of this paper, we will rely on the intuition developed from numerical studies of the oblique instability and argue that the beam ends up heating the IGM primarily.

In what follows, we will presume that the nonlinear evolution of the "oblique" or related plasma instabilities leads to the heating of the background IGM, i.e., beam cooling. However, we note that the beam itself may be the primary recipient of this kinetic energy, i.e., beam disruption. Most of our results depend critically on beam cooling and not beam disruption. This includes the effect of blazar heating on the IGM temperature–density relation studied in Paper II, the effects on structure formation studied in Paper III, the excellent reproduction of the statistical properties of the high-redshift Lyα forest found in Puchwein et al. (2012), and the implications on the EGRB and the redshift evolution of TeV blazars studied in Section 5. However, our conclusions on the inapplicability of IGMF constraints determined from the non-observation of GeV emission from blazars still remain in the presence of beam disruption. This is because the self-scattering of pairs in the beam would suppress this GeV emission in a similar manner to an IGMF.

The beam instabilities we have discussed have been analyzed primarily within the context of unmagnetized plasmas. However, for a variety of theoretical reasons a weak IGMF is not unexpected. For example, a field strength of 10⁻¹⁵ G is sufficient to explain the observed nG galactic fields via compression and winding alone. Here, we consider the implications that an IGMF has for the instability growth rates we have described above.

A strong IGMF causes the ultrarelativistic pairs to gyrate and therefore to isotropize, suppressing the growth of instabilities that feed on the beam anisotropy (e.g., those we have described above). However, for this to efficiently quench the growth of the plasma beam instabilities, this isotropization must occur on a timescale comparable to the instability growth time, i.e., the Larmor frequency must be comparable to the cooling rate, Γ. This condition gives a lower limit on IGMF strengths sufficient to appreciably suppress the growth of plasma beam instabilities of

$$\begin{equation} B \gtrsim 10^{-12} \left(\frac{\gamma }{10^6}\right) \left(\frac{\Gamma }{10^{-4}\,{\rm yr}}\right)\,{\rm G}\,, \end{equation} \tag{ 21 }$$

considerably larger than those both typical of primordial formation mechanisms (Widrow 2002) and implied by galactic magnetic field estimates, assuming galactic fields are produced by contraction and winding alone.

Because the VHEGRs emitted by the TeV blazars travel cosmological distances prior to producing pairs, an IGMF capable of suppressing the plasma beam instabilities must necessarily have a volume-filling fraction close to unity. In particular, it must permeate the low-density regions, where most of the cooling occurs. However, the Alfvén velocity within those areas is extraordinarily small, roughly

$$\begin{equation} v_{\rm A}\simeq 4\times 10^{-3} \left(\frac{B}{10^{-12}\,{\rm G}}\right) \left(1+\delta \right)^{-1/2}\,{\rm k}{\rm m}\,{\rm s}^{-1}. \end{equation} \tag{ 22 }$$

The IGM sound speed,

$$\begin{equation} c_s \simeq 10 \left(\frac{T}{10^4\,{\rm K}}\right)^{1/2}\,{\rm k}{\rm m}\,{\rm s}^{-1}\,, \end{equation} \tag{ 23 }$$

is considerably larger, implying that convection is much more efficient. Nevertheless, even after substantial heating via the thermalization of the TeV-blazar emission (Paper II) is taken into account, magnetic fields will have propagated ≲ 100 kpc over a Hubble time via convection and much less via diffusion, implying that a pervasive, sufficiently strong magnetic field cannot be produced via ejection from galactic dynamos. While galactic winds can produce much faster outflows, v_W ∼ 10³ km s⁻¹, unless they inject a mass comparable to that contained in the low-density regions over a Hubble time, they are rapidly slowed via dissipation at shocks in the IGM, again limiting the spread of galactic fields. Moreover, we note that since the magnetic field must be volume filling to suppress the plasma beam cooling substantially, it is insufficient to produce pockets of strong fields, and thus any galactic origin powered by winds requires a nearly complete reprocessing of the low-density regions. However, the ultimate thermalization of such fast, dense winds would raise the IGM temperature to ∼10⁸ K, in conflict with the Lyα forest data and exceeding the entire bolometric output of quasars by at least a factor of two. For these reasons, we conclude that a volume-filling strong IGMF would demand a primordial origin.

Blazars dominate the extragalactic gamma-ray sky and thus represent the best-studied VHEGR source class. Since we are interested exclusively in those objects responsible for the bulk of the VHEGR emission, we necessarily concentrate on the subset of blazars that are luminous VHEGR emitters. Of the 28 TeV sources listed in Table 1, 23 are peaked at very high energies (the HBL and hard IBL sources; for a full definition see below), and constitute what we will call collectively TeV blazars. This necessarily is limited to z ≃ 0.1 as a consequence of the annihilation of VHEGRs on the EBL. Nevertheless, we will find that this is remarkably similar to the local quasar luminosity function (ϕ_Q) and thus attempt to construct a blazar luminosity function by analogy:

$$\begin{equation} \phi _B(z,L)\equiv \frac{d^2 \Phi _B}{dz\, d\log _{10}L}\,, \end{equation} \tag{ 25 }$$

in which Φ_B is the comoving number density of blazars at redshifts less than z with isotropic-equivalent luminosities below L. This is different from previous efforts to empirically constrain ϕ_B from the EGRB (see, e.g., Narumoto & Totani 2006; Inoue & Totani 2009) in at least two ways. First, we begin with an empirically determined local ϕ_B and attempt to extend this by placing the TeV blazars within the broader context of accreting supermassive black holes, rather than beginning with the EGRB and working backward to infer an acceptable ϕ_B. Second, we are primarily concerned with ϕ_B of TeV blazars and specifically do not consider the contributions from other kinds of objects. While this does not represent a significant oversight in terms of the high-energy contributions to the EGRB, which almost certainly arises from the VHEGR-emitting blazars, it does mean that our conclusions regarding the luminosity function of the TeV blazars do not necessarily apply to all Fermi AGNs (e.g., the flat-spectrum radio quasars (FSRQs); see below).

5.1.1. Placing TeV Blazars in Context

5.1.2. An Empirical Local ϕ_B(z,L)

While we have attempted to place the observed TeV blazars, which we have identified with the HBLs and VHEGR-emitting IBLs, into the broader context of AGNs using the unified model, due to the substantial distinctions in accretion rate, emission properties, object morphology, and geometry, it is not obvious that any of the properties of TeV blazars should be similar to those of AGNs more generally. Nevertheless, evidence for a simple connection between the two populations can be found in the similarity between their luminosity functions (a fact we will exploit later in estimating the redshift evolution of the TeV blazars). Here, we define the luminosity for the purposes of defining ϕ_B to be the isotropic-equivalent value associated with emission between 100 GeV and 10 TeV. While this may be considered to be a VHEGR luminosity, because most TeV blazars are peaked within this band, this corresponds to the majority of the emission from these sources.

The objects listed in Table 1 were chosen because they have well-defined SEDs, based on a combination of VERITAS, H.E.S.S., and MAGIC observations. These 28 sources have VHEGR spectra that are well fit by the form

$$\begin{equation} \frac{dN}{dE} = f_0 \left(\frac{E}{E_0}\right)^{-\alpha }, \end{equation} \tag{ 26 }$$

where f₀ is the normalization in units of cm⁻² s⁻¹ TeV⁻¹. The gamma-ray energy flux is trivially related to dN/dE by F_E = EdN/dE∝E^{1 − α}, from which we obtain a VHEGR flux

$$\begin{equation} F = E_0 f_0 \int _{100\,{\rm G}{\rm eV}}^{10\,{\rm T}{\rm eV}} dE\,\left(\frac{E}{E_0}\right)^{1-\alpha }\,, \end{equation} \tag{ 27 }$$

and for sources with a measured redshift a corresponding isotropic-equivalent luminosity, L = 4πD²_LF, where D_L is the luminosity distance.²² The resulting f₀, E₀, α, F, and L are collected in Table 1. In addition, we list the redshift, inferred comoving distance, and absorption-corrected intrinsic spectral index, defined at E₀, obtained via

$$\begin{equation} \hat{\alpha } = {-}\left.\frac{d\ln E^{-\alpha } e^{\tau _E[E(1+z),z]}}{d\ln E}\right|_{E_0} \simeq \alpha - \tau _{E}\left[E_0(1+z),z\right]\,. \end{equation} \tag{ 28 }$$

For high-redshift sources $\hat{\alpha }$ can be substantially less than 2, implying that an intrinsic spectral upper cutoff must exist.

To produce ϕ_B, we must account for a variety of selection effects inherent in the sample listed in Table 1. The objects in Table 1 were originally selected for study for a variety of source-specific reasons, e.g., existing well-known sources, extremely high X-ray-to-radio flux ratio in the Sedentary High-energy-peaked BL Lac catalog, hard spectrum sources in the Fermi point-source catalog, and flagged as promising by the Fermi-LAT collaboration. In addition, the source selection suffered from the usual problems associated with surveys (e.g., scheduling conflicts with other targets, moon, bad weather, etc.). As a consequence, this sample is somewhat inhomogeneous. Nevertheless, in lieu of a less-biased sample, we will treat it as homogeneous and correct for the selection effects where possible, focusing on those due to the sky coverage and duty cycle of TeV observations and those due to sensitivity limits of current imaging atmospheric Cerenkov telescopes.

To estimate the sky completeness and duty cycle of this set of objects, we rely on the all-sky GeV gamma-ray observations of HBL and IBL sources by Fermi (Abdo et al. 2010d). Outside of the Galactic plane, Fermi observes 118 high-synchrotron-peaked (HSP) blazars and a total of 46 intermediate-synchrotron-peaked (ISP) blazars. Roughly half of the latter are likely to emit VHEGRs as indicated by their flat spectral index between 0.1 and 100 GeV (α ≲ 2; see the spectral index distribution (SID) of Figure 14 in Abdo et al. 2009). Of these potential 141 TeV blazars, only 22 have also been coincidentally identified as TeV sources, whereas there are a total of 33 known TeV blazars (29 HBL, 4 IBL). If these 141 sources are all TeV emitters but have not been detected due to incomplete sky coverage of current TeV instruments, then the selection factor is η_sel = 141/33 = 4.3. In addition, the duty cycle of coincident GeV and TeV emission is η_duty = 33/22 = 1.5.²³ Finally, by excluding the Galactic plane for galactic latitudes b < 10°, this is an underestimate by roughly η_sky = 1.17.

Zoom In Zoom Out Reset image size

We make a rough attempt to correct for the sensitivity limit of the TeV instruments, inferring a flux limit by fitting the upper envelope of the TeV-blazar flux–luminosity distance distribution. Somewhat surprisingly, despite the heterogeneous nature of the TeV observations, they are remarkably well described by a single flux limit: 4.19 × 10⁻¹² erg cm⁻² s⁻¹. This process and the associated limit are shown explicitly in the inset of Figure 5. From this we obtain a maximum redshift, and therefore comoving volume, associated with each luminosity. We then construct ϕ_B by counting the number of TeV blazars per unit log₁₀L in each logarithmic luminosity bin and dividing by the comoving volume. The resulting ϕ_B, weighted by luminosity (and therefore showing the luminosity density in comoving units), is shown in Figure 5. It peaks at ∼4 × 10⁴⁴ erg s⁻¹, implying that, as expected, these objects are systematically dimmer than most other AGNs, and exhibits the broken power-law shape typical of AGN luminosity functions. More importantly, despite a handful of sources with redshifts ∼0.5, the objects in Table 1 are all nearby, and therefore our ϕ_B corresponds to that in the local universe, i.e., z ∼ 0.1. Based on this, the inferred present-day TeV-blazar luminosity density is roughly 1 × 10³⁸ erg s⁻¹ Mpc⁻³.

Also shown in Figure 5 is the Lϕ_Q obtained by Hopkins et al. (2007). After rescaling Lϕ_Q to lower luminosities (0.55) and lower-luminosity densities (2.1 × 10⁻³), it provides a remarkably good fit to our Lϕ_B (reduced χ² of 0.13 with 3 degrees of freedom), i.e., we find

$$\begin{equation} \phi _B(0.1,L) \simeq 3.8\times 10^{-3} \phi _Q(0.1,1.8 L)\,, \end{equation} \tag{ 29 }$$

where an explicit expression for ϕ_Q is given in Appendix B. While there is considerable uncertainty in Lϕ_B, especially at low luminosities, it clearly does not fit Lϕ_Q at higher z. This suggests two immediate conclusions.

• 1.  
  The bolometric output of TeV blazars and quasars is regulated by similar mechanisms, presumably accretion, despite the large difference in luminosity and the details of the emission processes between the two populations.
• 2.  
  TeV blazars and quasars are contemporaneous elements in a single AGN distribution; specifically, TeV-blazar activity does not lag that of quasars.

5.1.3. Extending ϕ_B(z, L) to High-z

Based on the strong similarities between ϕ_B and ϕ_Q and the associated implications, we make the conservative theoretical assumption that the redshift evolution of the TeV blazars follows that of quasars. That is, we suppose that Equation (29) holds at all z. This implies that the integrated comoving TeV-blazar isotropic-equivalent luminosity density is given by

$$\begin{equation} \Lambda _B(z) \equiv \int _{-\infty }^\infty L \phi _B(z,L)\,d\log _{10}L = \eta _B \Lambda _Q(z)\,, \end{equation} \tag{ 30 }$$

where the constant of proportionality, η_B ≃ 2.1 × 10⁻³, is then set by the comparison between ϕ_B and ϕ_Q at z = 0.1. As a consequence, in comoving units, the TeV-blazar luminosity density would be roughly an order of magnitude larger at z ≃ 2, ∼1 × 10³⁹ erg s⁻¹ Mpc⁻³.

The redshift distribution of the Fermi BL Lac sample (i.e., the First LAT Catalog, 1LAC; Abdo et al. 2010d) is peaked at z < 0.2 and falls rapidly thereafter. Inasmuch as the Fermi HSP and ISP counts directly probe the low-z ϕ_B, it may appear that they are inconsistent with the rapidly evolving ϕ_B described in the previous section. Indeed, an analysis of the first three months of Fermi observations suggested that the BL Lac population does not grow substantially with redshift (Abdo et al. 2009). However, as seen in Figure 6, this is not necessarily the case.

Zoom In Zoom Out Reset image size

Assuming, as we have, that the number of TeV blazars is equal to the number of hard Fermi gamma-ray blazars, corresponding to ISPs with α ≲ 2 (composing roughly half of that class) and HSPs, the expected number observed inside of a given redshift is

$$\begin{eqnarray} {\mathcal N}_B(z) &=& \int _0^z dz^{\prime } \frac{dD}{dz^{\prime }} 4\pi D_A^2\nonumber \\ &&\times \int _{\log _{10}L_{\rm min}(z^{\prime })}^{\log _{10}L_{\rm max}} (1+z^{\prime })^3 \phi _B(z^{\prime },L) \,d\log _{10}L{,}\quad\quad \end{eqnarray} \tag{ 31 }$$

where D ≡ ∫cdt' = ∫cdz'/[H(z')(1 + z')] and D_A = D_L/(1 + z)² are the proper and angular diameter distances, respectively, L_max ≃ 2 × 10⁴⁶ erg s⁻¹ is the intrinsic upper cutoff of the TeV-blazar isotropic-equivalent luminosity function, and

$$\begin{eqnarray} L_{\rm min}(z) &=& F_{\rm min} 4\pi D_L^2(z) (1+z)^{\alpha -2}\nonumber \\ &\simeq & 4\times 10^{46}\,\eta _{\rm min}\left(\frac{1+z}{2}\right)^{\alpha -2} \left[\frac{D_L(z)}{D_L(1)}\right]^{2} \,{\rm erg}\,{\rm s}^{-1}\quad\quad \end{eqnarray} \tag{ 32 }$$

is the lower limit set by the Fermi flux limit (see Figure 23 of Abdo et al. 2010d, and surrounding discussion). The factor η_min is a correction relating the 100 GeV–10 TeV isotropic-equivalent luminosities we employ to define ϕ_B to the 100 MeV–100 GeV luminosities used by Fermi to define the flux limit. The dependence on α arises from the limited spectral coverage of both definitions, though here we fix α = 3 for all objects based on the sources that dominate the TeV flux at Earth. Note that since the ϕ_Q from Hopkins et al. (2007) diverges at small L, the lower-luminosity cutoff is critical to getting both the total number of TeV blazars and the shape of their redshift distribution correct. The resulting ${\mathcal N}_B$ is shown in Figure 6 for three different choices of η_min: 0.78, 1.6, and 3.1, of which η_min = 1.6 is most similar to the 1LAC hard gamma-ray blazars.

Generally, our ϕ_B does an excellent job of reproducing the overall number of 1LAC hard gamma-ray blazars and the dominance of nearby objects in their redshift distribution. This is despite the strong redshift evolution of ϕ_B implied by its relationship with ϕ_Q. The reason for this is the flat distribution at low L, the steep drop-off at high L (a result of which is that the shape is only marginally sensitive to the cutoff at L_max), and the rapidly growing L_min due to the fixed flux limit.

However, we note that the comparison between ϕ_B and the 1LAC hard gamma-ray blazar statistics assumes that those with measured redshifts, composing roughly half of the 1LAC hard gamma-ray blazar sample, are representative of the hard gamma-ray blazar population as a whole. This appears not to be the case; as we discuss in detail in Appendix C, it is clear that the 1LAC HSPs with and without measured redshifts are not drawn from the same underlying α-distribution. Based on a similar analysis for 1LAC BL Lac objects generally, Abdo et al. (2010d) argued that the objects without redshifts are more consistent with that z > 0.5 population. However, this conclusion does not extend to HSPs, for which there are only three sources in the clean 1LAC HSP sample with z > 0.5. Rather, in Appendix C, we show that the 1LAC HSPs without redshifts are distributed in both spectral index and flux much more similarly with nearby HSPs (z < 0.25) than more distant HSPs (z ⩾ 0.25). If this is because these sources are intrinsically underluminous, nearby objects, it could marginally improve the already remarkable comparison between the number of objects implied by our ϕ_B and those observed. However, because the number of predicted 1LAC hard gamma-ray blazars is strongly dependent on the flux limit, even a marginal increase in either ϕ_B at low luminosities or the effective flux limit results in a $d{\mathcal N}_B/dz$ that is more strongly weighted at low z, potentially allowing even more dramatically evolving luminosity functions. For this reason, we conclude that our ϕ_B is broadly consistent with the 1LAC hard gamma-ray blazar distribution and that Fermi does not, at present, exclude a rapidly evolving ϕ_B in the recent past.

While the statistics of the 1LAC hard gamma-ray blazar sample probes the evolution of ϕ_B at z ≲ 0.5, the Fermi EGRB is sensitive to TeV blazars at high z. The contribution to the EGRB from TeV blazars given our ϕ_B is shown in Figure 7. Because ICCs are suppressed, computing the EGRB flux requires only summing the individual intrinsic GeV spectra of the TeV blazars. For simplicity, to estimate the TeV-blazar contribution to the EGRB, we assume that all TeV blazars have identical intrinsic spectra, given by a broken power law,

$$\begin{equation} F_E = L\hat{F}_E \propto \frac{1}{\left(E/E_b\right)^{\alpha _L-1} + \left(E/E_b\right)^{\alpha -1}}\,, \end{equation} \tag{ 33 }$$

where the normalization is set such that the $\int _{100\,{\rm G}{\rm eV}}^{10\,{\rm T}{\rm eV}} \hat{F}_E dE = 1$ (i.e., $\int _{100\,{\rm G}{\rm eV}}^{10\,{\rm T}{\rm eV}} L\hat{F}_E dE = L$, corresponding to the luminosity used to define ϕ_B and thus determine η_B), E_b is the energy of the spectral break, and α_L < α are the low- and high-energy spectral indexes, respectively.²⁴ With this, after performing the integral over L, the EGRB flux is then

$$\begin{eqnarray} E^2 \frac{dN}{dE}(E,z) &=& \frac{\eta _B}{4\pi } \int _z^\infty \! dz^{\prime } \,\frac{dD}{dz^{\prime }} \,4\pi D_A^2 \frac{\tilde{\Lambda }_Q(z^{\prime }) \hat{F}_{E^{\prime }}}{4\pi D_L^2}{\rm e}^{-\tau _E[E^{\prime },z^{\prime }]}\nonumber\\ &=& \frac{\eta _B}{4\pi } \int _z^\infty \! dz^{\prime } \frac{c\tilde{\Lambda }_Q(z^{\prime })}{H(z^{\prime })(1+z^{\prime })^5} \,\hat{F}_{E^{\prime }}\,{\rm e}^{-\tau _E[E^{\prime },z^{\prime }]}\,,\nonumber\\ \end{eqnarray} \tag{ 34 }$$

where E' = E(1 + z') and $\tilde{\Lambda }_Q = \Lambda _Q (1+z)^3$ is the physical quasar luminosity density. We choose fiducial values for the spectral parameters of α_L = 1.33, E_b = 1 TeV, and α = 3, typical of the TeV blazars (see Table 1 and Ghisellini 2011). Finally, since Fermi resolves individual gamma-ray blazars with isotropic-equivalent luminosities ∼10⁴⁵ erg s⁻¹, roughly the location of the peak in ϕ_B, for z ≲ 0.25, we construct the EGRB from sources at larger redshift.

Zoom In Zoom Out Reset image size

If E_b is sufficiently high (≳ 2 × 10² GeV), the EGRB is remarkably insensitive to the particular values of E_b and α (see the bottom two panels of Figure 7). This is a direct result of the annihilation of the VHEGRs on the EBL, effectively removing the relevant portion of the blazar SED from the EGRB. More important is the low-energy spectral index. However, for values of α_L that are consistent with the Comptonization models for the VHEGR emission, the anticipated EGRB is consistent with the Fermi result. Thus, we find that despite our dramatically evolving ϕ_B, we are able to satisfy the limits imposed by the Fermi EGRB. Moreover, for reasonable spectral-parameter values it is possible to accurately reproduce both the magnitude and shape of the high-energy EGRB (i.e., at energies >10 GeV).²⁵ At first this may appear in conflict with other studies that have performed more sophisticated analyses of the gamma-ray emission from blazars (e.g., Narumoto & Totani 2006; Kneiske & Mannheim 2008; Inoue & Totani 2009; Venters 2010) that have typically found that such a strongly evolving ϕ_B overproduces the Fermi EGRB. However, this is not the case for three reasons.

First, the AGN populations typically considered in computations of the EGRB include the FSRQs and LBLs, which dominate at low energies. The TeV blazars of interest here are only the high-energy tail of the blazar distribution, and it is these alone that we fix to the quasar luminosity function (indeed, it is for only these objects that the arguments surrounding Figure 5 have been made). Moreover, the TeV blazars themselves are generally dimmer, both bolometrically and within the Fermi LAT energy band (200 MeV–100 GeV), than the LBLs and FSRQs, which have typical isotropic-equivalent luminosities of 3 × 10⁴⁷ erg s⁻¹ (Narumoto & Totani 2006).

Second, and perhaps most importantly, the suppression of the ICCs means that the VHEGR emission is lost to heat in the IGM, not reprocessed to below 100 GeV. Since the TeV-blazar flux is presumed to be dominated by emission above 100 GeV, the non-radiative dissipation of this emission substantially reduces the impact of the TeV blazars on the EGRB. For all of the spectral parameters we considered the unabsorbed fluxes exceed the 30 GeV EGRB by a significant margin, implying that the suppression of the ICCs is crucial to bringing the TeV-blazar contributions in line with the Fermi EGRB. Thus, the lack of the ICCs generally appears necessary to reconcile the blazar and quasar luminosity functions.

Third, the VHEGR spectra from the most luminous TeV blazars necessarily peak near or above TeV energies. For a number of the sources in Table 1, the intrinsic VHEGR spectrum (adjusted for annihilation on the EBL) is inverted, indicating that this break is above a TeV. The brightest TeV blazar in Table 1, Mkn 421, has an inverted spectrum below 100 GeV (Abdo et al. 2010d), implying E_b ∼ 500 GeV. The second brightest, 1ES 1959+650, is consistent with being flat below 100 GeV and appears to peak near 1 TeV when in the high state (Daniel et al. 2005). Similarly, 1ES 2344+514 peaks near ∼300 GeV (Albert et al. 2007c) and Mkn 501 is inverted below 100 GeV, implying a peak above that value (Abdo et al. 2010d). Thus, the sources likely to dominate the VHEGR background, and by extension the EGRB, all turn over near 100 GeV–1 TeV. This is expected if the VHEGR emission arises from inverse Compton scattering the synchrotron bump (see, e.g., Ghisellini 2011, and references therein). The consequence of the spectral break is the suppression of the contribution from the TeV blazars to the EGRB below E_b. Were the TeV-blazar contribution to the EGRB dominated by sources with E_b ≲ 2 × 10² GeV, it would exceed the Fermi EGRB at E ≲ 10 GeV.

A more complete analysis of the EGRB should include a variety of spectra, including a distribution of break energies, smoothly connecting the HBL, IBL, and LBL populations. However, this is beyond the scope of this paper, which is primarily concerned with the fate of the VHEGR emission, absent in the LBLs and FSRQs. Nevertheless, a recent comprehensive model, which includes the contributions of the FSRQs and BL Lac objects observed by Fermi, as well as starburst galaxies, has had considerable success fitting the Fermi EGRB with an even more extremely evolving ϕ_B, fixing it to the luminosity function of radio galaxies (Cavadini et al. 2011; Stecker & Venters 2011). Of particular relevance here is that Cavadini et al. (2011) and Stecker & Venters (2011) explicitly ignored the potential contributions of ICCs. Unlike their analysis, however, we find that the TeV blazars are capable of reaching the highest Fermi bands despite the annihilation on the EBL.

The two-stream instability arises due to the excitation of negative-energy electrostatic waves in the beam and target plasmas. These waves carry away both the energy and momentum of the beam. Specifically, we compute the growth rate of the electrostatic wave moving in the direction of the beam in the absence of a background magnetic field.²⁶ In this case, the Vlasov equations for the electrons and positrons are

$$\begin{eqnarray} \frac{D f^-}{Dt} - e\mbox{\boldmath $\rm E$}\cdot \frac{\partial f^-}{\partial \mbox{\boldmath $\rm p$}} &=& 0\nonumber\\ \frac{D f^+}{Dt} + e\mbox{\boldmath $\rm E$}\cdot \frac{\partial f^+}{\partial \mbox{\boldmath $\rm p$}} &=& 0\,, \end{eqnarray} \tag{ A1 }$$

where $D/Dt\equiv \partial /\partial t + \mbox{\boldmath $\rm v$}\cdot \mbox{\boldmath $\rm \nabla $}$ is the Stokes derivative and $\mbox{\boldmath $\rm E$}$ is the net electric field. Linearizing these and Fourier transforming in t and $\mbox{\boldmath $\rm x$}$ gives

$$\begin{eqnarray} {-}i(\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}) f^-_1 - e \mbox{\boldmath $\rm E$}_1\cdot \frac{\partial f^-_0}{\partial \mbox{\boldmath $\rm p$}} &=& 0\nonumber\\ {-}i(\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}) f^+_1 + e \mbox{\boldmath $\rm E$}_1\cdot \frac{\partial f^+_0}{\partial \mbox{\boldmath $\rm p$}} &=& 0{,} \end{eqnarray} \tag{ A2 }$$

where f^∓₁ are the perturbations to the electron and positron distribution functions, $\mbox{\boldmath $\rm E$}_1$ is the electrostatic wave field, and we have assumed that the background distributions, f^∓₀, are isotropic. As a result, we have

$$\begin{equation} f^\mp _1 = \frac{\pm ie}{\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}} \mbox{\boldmath $\rm E$}_1\cdot \frac{\partial f^\pm _0}{\partial \mbox{\boldmath $\rm p$}}\,. \end{equation} \tag{ A3 }$$

At this point we may compute the dielectric tensor associated with the plasma response; however, it will suffice to consider Gauss's law. Thus, we now compute the perturbed charge density:

$$\begin{eqnarray} \rho _1 &=& \int \left(e f^+_1 - e f^-_1\right) d^3\!p\nonumber\\ &=& -ie^2\mbox{\boldmath $\rm E$}_1\cdot \int \frac{1}{\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}} \left(\frac{\partial f^-_0}{\partial \mbox{\boldmath $\rm p$}} + \frac{\partial f^+_0}{\partial \mbox{\boldmath $\rm p$}} \right)d^3\!p\nonumber\\ &=& ie^2\mbox{\boldmath $\rm E$}_1\cdot \int \left(f^-_0 + f^+_0 \right) \left(\frac{\partial }{\partial \mbox{\boldmath $\rm p$}} \frac{1}{\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}} \right)d^3\!p \nonumber\\ &=& \frac{ie^2}{m_e}\mbox{\boldmath $\rm E$}_1\cdot \int \left(f^-_0 + f^+_0 \right) \frac{(\mbox{\boldmath $\rm k$}- \mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}\mbox{\boldmath $\rm v$})}{\gamma (\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$})^2}\,d^3\!p\,. \end{eqnarray} \tag{ A4 }$$

It is now necessary to specify the f^∓₀. We idealize the target ionic plasma as cold and the pair beam plasma as mono-energetic, yielding

$$\begin{equation} f^-_0 = n_t \delta ^3(\mbox{\boldmath $\rm p$}) + \frac{n_{b}}{2}\delta ^3(\mbox{\boldmath $\rm p$}-\mbox{\boldmath $\rm p$}_b) \quad {\rm and}\quad f^+_0 = \frac{n_{b}}{2} \delta ^3(\mbox{\boldmath $\rm p$}-\mbox{\boldmath $\rm p$}_b)\,, \end{equation} \tag{ A5 }$$

where n_t and n_b are the lepton number densities in the target and beam, respectively. After performing the trivial integrals, we then obtain

$$\begin{eqnarray} \rho _1 &=& \frac{i e^2}{m_e} \mbox{\boldmath $\rm E$}_1 \cdot \left(n_t \frac{\mbox{\boldmath $\rm k$}}{\omega ^2} + n_{b}\frac{\mbox{\boldmath $\rm k$}-\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}_b\mbox{\boldmath $\rm v$}_b}{\gamma _b(\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}_b)^2} \right)\nonumber\\ &=& \frac{i e^2}{m_e} \mbox{\boldmath $\rm E$}_1\cdot \mbox{\boldmath $\rm k$}\left(\frac{n_t}{\omega ^2} + \frac{n_{b}}{\gamma _b^3(\omega -k v_b)^2} \right)\,, \end{eqnarray} \tag{ A6 }$$

where in the final equality we used the fact that $\mbox{\boldmath $\rm k$}\parallel \mbox{\boldmath $\rm v$}_b$. From Gauss's law we have $i \mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm E$}_1 = 4\pi \rho _1$, and therefore

$$\begin{equation} 1 - \frac{\omega _{P,t}^2}{\omega ^2} - \frac{\omega _{P,b}^2}{\gamma _b^3(\omega -kv_b)^2} = 0\,, \end{equation} \tag{ A7 }$$

where ω²_{P, t} ≡ 4πe²n_t/m_e and ω²_{P, b} ≡ 4πe²n_b/m_e are the plasma frequencies associated with the target and beam plasmas.

This explicitly provides the dispersion relation, quadratic in ω (one electrostatic wave traveling in each direction for each plasma). When n_b = 0, ω = ω_{P, t}. When n_b ≪ n_t, as is the case of interest here, we may solve the dispersion relation perturbatively. We do this by setting ω = ω_{P, t}(1 + η), with η ≪ 1, which gives

$$\begin{equation} 2 \eta - \frac{\omega _{P,b}^2}{\gamma _b^3(\omega _{P,t}-kv_b+\eta \omega _{P,t})^2} =0\,, \end{equation} \tag{ A8 }$$

which is no longer independent of k. When |ω_{P, t} − kv_b| ≫ |ηω_{P, t}|, η is real and thus there is no instability. On the other hand, where k ≃ ω_{P, t}/v_b, we have

$$\begin{equation} 2\eta - \frac{\omega _{P,b}^2}{\gamma _b^3\omega _{P,t}^2\eta ^2} = 0\,, \end{equation} \tag{ A9 }$$

and therefore η³ = n_b/2n_tγ³_b. This has three solutions:

$$\begin{equation} \eta = \frac{1}{\gamma _b} \left(\frac{n_{b}}{2n_t}\right)^{1/3} \left\lbrace 1\,,\quad -\frac{1}{2}-\frac{\sqrt{3}}{2}i\,,\quad -\frac{1}{2}+\frac{\sqrt{3}}{2}i \right\rbrace , \end{equation} \tag{ A10 }$$

the first of which is oscillatory, the second is decaying, and the third is growing with timescale

$$\begin{equation} \Im \left(\omega \right) = \Im (\omega _{P,t}\eta ) = \frac{\sqrt{3}}{2}\left(\frac{n_{b}}{2n_t}\right)^{1/3}\frac{\omega _{P,t}}{\gamma _b}. \end{equation} \tag{ A11 }$$

This differs from that associated with non-relativistic, ionic beams only by the factor of 1/γ_b, arising due to time dilation within the beam.

Note that since the energy within the electrostatic wave is proportional to E²₁, the rate at which energy is removed from the beam is Γ_TS ≡ 2ℑ(ω).

The relativistic version of the Weibel instability has been discussed in some detail in the literature for the case of anisotropic but symmetric beams (Weibel 1959; Fried 1959; Yoon & Davidson 1987; Medvedev & Loeb 1999). Here, we consider a case of a dilute pair beam incident upon a much denser target plasma. This is essentially identical to the two-stream instability discussed previously, though coupling instead to an electromagnetic mode with $\mbox{\boldmath $\rm k$}\perp \mbox{\boldmath $\rm v$}_b$.

Again we begin with the linearized Vlasov equations for the electrons and positrons:

$$\begin{eqnarray} \frac{Df^-_1}{Dt} - e \left(\mbox{\boldmath $\rm E$}_1+\mbox{\boldmath $\rm v$}\times \mbox{\boldmath $\rm B$}_1\right)\cdot \frac{\partial f^-_0}{\partial \mbox{\boldmath $\rm p$}} &=& 0\nonumber\\ \frac{Df^+_1}{Dt} + e \left(\mbox{\boldmath $\rm E$}_1+\mbox{\boldmath $\rm v$}\times \mbox{\boldmath $\rm B$}_1\right)\cdot \frac{\partial f^+_0}{\partial \mbox{\boldmath $\rm p$}} &=& 0\,, \end{eqnarray} \tag{ A12 }$$

where we have $\mbox{\boldmath $\rm B$}_1 = \mbox{\boldmath $\rm k$}\times \mbox{\boldmath $\rm E$}_1/\omega$ from Faraday's law. Fourier transforming in t and $\mbox{\boldmath $\rm x$}$ and solving for f^∓₁ gives

$$\begin{eqnarray} f_1^\mp &=& \pm \frac{ie\left[\mbox{\boldmath $\rm E$}_1+\mbox{\boldmath $\rm v$}\times \left(\mbox{\boldmath $\rm k$}\times \mbox{\boldmath $\rm E$}_1\right)/\omega \right]}{\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}} \cdot \frac{\partial f^\mp _0}{\partial \mbox{\boldmath $\rm p$}}\nonumber\\ &=& \pm \frac{ie}{\omega } \left(\mbox{\boldmath $\rm E$}_1 + \frac{\mbox{\boldmath $\rm v$}\cdot \mbox{\boldmath $\rm E$}_1}{\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}}\mbox{\boldmath $\rm k$}\right) \cdot \frac{\partial f^\mp _0}{\partial \mbox{\boldmath $\rm p$}}\,. \end{eqnarray} \tag{ A13 }$$

The associated perturbation to the current is

$$\begin{eqnarray} \mbox{\boldmath $\rm j$}_1 &=& \int \left(e f^+_1 - e f^-_1 \right)\mbox{\boldmath $\rm v$}d^3\!p \nonumber\\ &=& -\frac{ie^2}{\omega }\int \frac{\partial (f^+_0+f^-_0)}{\partial \mbox{\boldmath $\rm p$}} \cdot \left(\mbox{\boldmath $\rm E$}_1 + \frac{\mbox{\boldmath $\rm v$}\cdot \mbox{\boldmath $\rm E$}_1}{\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}}\mbox{\boldmath $\rm k$}\right) \mbox{\boldmath $\rm v$}d^3\!p \nonumber\\ &=& \frac{ie^2}{\omega }\int \frac{f_0}{m_e \gamma } \Bigg[ \hat{\mbox{\boldmath $\rm 1$}} + \frac{\mbox{\boldmath $\rm k$}\mbox{\boldmath $\rm v$}+\mbox{\boldmath $\rm v$}\mbox{\boldmath $\rm k$}}{\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$}}\nonumber\\ &&- \frac{\omega ^2-k^2}{(\omega -\mbox{\boldmath $\rm k$}\cdot \mbox{\boldmath $\rm v$})^2} \mbox{\boldmath $\rm v$}\mbox{\boldmath $\rm v$}\Bigg] \cdot \mbox{\boldmath $\rm E$}_1 d^3\!p\,, \end{eqnarray} \tag{ A14 }$$

in which we have defined f₀ ≡ f⁺₀ + f⁻₀.

Choosing the f^∓₀ given by Equation (A5) and assuming $\mbox{\boldmath $\rm k$}\parallel \mbox{\boldmath $\rm v$}_b\parallel \mbox{\boldmath $\rm E$}_1$, we recover the standard two-stream instability. For our purposes here, however, the computations may be substantially simplified by boosting into a frame in which the target plasma is not at rest. In this frame we have f₀ given by

$$\begin{equation} f_0 = n_t^{\prime } \delta ^3(\mbox{\boldmath $\rm p$}^{\prime }-\mbox{\boldmath $\rm p$}_t^{\prime }) + n_{b}^{\prime } \delta ^3(\mbox{\boldmath $\rm p$}^{\prime }-\mbox{\boldmath $\rm p$}_b^{\prime })\,, \end{equation} \tag{ A15 }$$

with $\mbox{\boldmath $\rm p$}_t^{\prime },\,\mbox{\boldmath $\rm p$}_b^{\prime }\parallel \mbox{\boldmath $\rm v$}_b^{\prime }$, where we will take care at the end to properly relate all of these quantities to their target-frame analogs. In particular, we choose this frame such that $\tilde{n}_t v_t^{\prime }=\tilde{n}_b v_b^{\prime }$, where $\tilde{n}_t$ and $\tilde{n}_b$ are the proper densities of the target and beam plasmas, respectively. If v_b is the beam plasma velocity in the target (or lab) frame, this gives

$$\begin{equation} v_t^{\prime } = \xi v_b^{\prime } = \xi \frac{v_b-v_t^{\prime }}{1-v_b v_t^{\prime }}\,, \end{equation} \tag{ A16 }$$

where $\xi \equiv \tilde{n}_b/\tilde{n}_t$ and is in our case much less than unity. Note that this is only the center of momentum frame if ξ = 1 (for which v'_t = v_b', the case most commonly discussed). If we choose $\mbox{\boldmath $\rm k$}^{\prime }\perp \mbox{\boldmath $\rm v$}_b^{\prime }$, this choice of frame causes the terms linear in $\mbox{\boldmath $\rm v$}_b^{\prime }$ in Equation (A14) to vanish identically, yielding

$$\begin{eqnarray} \mbox{\boldmath $\rm j$}_1^{\prime } &=& \frac{ie^2}{m_e\omega ^{\prime }} \Bigg[\Bigg(\frac{n_t^{\prime }}{\gamma _t^{\prime }}+\frac{n_{b}^{\prime }}{\gamma _b^{\prime }} \Bigg)\, \hat{\mbox{\boldmath $\rm 1$}} - \frac{\omega ^{\prime 2}-k^{\prime 2}}{\omega ^{\prime 2}}\nonumber\\ &&\times\Bigg(\frac{n_t^{\prime }}{\gamma _t^{\prime }}\mbox{\boldmath $\rm v$}_t^{\prime }\mbox{\boldmath $\rm v$}_t^{\prime } + \frac{n_{b}^{\prime }}{\gamma _b^{\prime }}\mbox{\boldmath $\rm v$}_b^{\prime }\mbox{\boldmath $\rm v$}_b^{\prime }\Bigg) \Bigg]\cdot \mbox{\boldmath $\rm E$}_1^{\prime }\,. \end{eqnarray} \tag{ A17 }$$

Using the inhomogeneous wave equation, obtained from Faraday's and Ampere's law,

$$\begin{equation} (k^{\prime 2} - \omega ^{\prime 2}) \mbox{\boldmath $\rm E$}_1^{\prime } = 4 \pi i \omega ^{\prime } \mbox{\boldmath $\rm j$}_1^{\prime }, \end{equation} \tag{ A18 }$$

and choosing $\mbox{\boldmath $\rm E$}_1^{\prime }\parallel \mbox{\boldmath $\rm v$}_b^{\prime }$ produces the dispersion relation

$$\begin{equation} \omega ^{\prime 2}-k^{\prime 2} - \frac{\omega _{P,t}^{\prime 2}}{\gamma _t^{\prime 3}} - \frac{\omega _{P,b}^{\prime 2}}{\gamma _b^{\prime 3}} - \frac{k^{\prime 2}}{\omega ^{\prime 2}} \left(\frac{\omega _{P,t}^{\prime 2} v_t^{\prime 2}}{\gamma _t^{\prime }} + \frac{\omega _{P,b}^{\prime 2} v_b^{\prime 2}}{\gamma _b^{\prime }}\right) = 0\,. \end{equation} \tag{ A19 }$$

At this point we make use of the fact that in our case ξ ≪ 1, which allows a perturbative solution of Equation (A16): $v_t^{\prime }=\xi v_b + \mathcal {O}(\xi ^2)$ and therefore $v_b^{\prime }=v_b(1-\xi)+\mathcal {O}(\xi ^2)$. Furthermore, $\omega ^{\prime } = \omega + \mathcal {O}(\xi ^2)$ and k' = k for the geometry under consideration. As consequence, ω'_{P, t}v_t'²/γ'_t ≃ ω_{P, t}ξ²v²_b = ξω_{P, b}v²_b/γ_b ≪ ω_{P, b}v²_b/γ_b and ω'²_{P, t}/γ_t'³ ≃ ω²_{P, t} = ξ⁻¹ω_{P, b}² ≫ ω'²_{P, b}/γ_b'³, which implies that

$$\begin{equation} \omega ^2 - k^2 - \omega _{P,t}^2 - \frac{k^2}{\omega ^2} \frac{\omega _{P,b}^2 v_b^2}{\gamma _b}(1-2\xi) = 0\,. \end{equation} \tag{ A20 }$$

This has a purely imaginary solution:

$$\begin{equation} \omega = i \sqrt{\frac{k^2+\omega _{P,t}^2}{2}} \sqrt{ \left[ 1 + \frac{4 k^2 \omega _{P,b}^2 v_b^2(1-2\xi)}{\big(k^2+\omega _{P,t}^2\big)^2\gamma _b} \right]^{1/2} - 1 }{.} \end{equation} \tag{ A21 }$$

For k ≪ ω_{P, t} this rises linearly with k, saturating for $k>\omega _{P,t}\gg \omega _{P,b}v_b^2/\sqrt{\gamma _b}$, giving the growth rate

$$\begin{equation} \Im \left(\omega \right) \simeq \frac{\omega _{P,b} v_b}{\sqrt{\gamma _b}}(1-\xi) + \mathcal {O}(\xi ^2)\,. \end{equation} \tag{ A22 }$$

Given that for our application ξ ≪ 1, we will drop terms first order in ξ as well. Note that the plasma frequency that enters into the growth rate is that of the beam, which is much lower than that associated with the target plasma. Nevertheless, the scale of the rapidly growing perturbations is limited to that associated with the target plasma, which is generally much smaller than that associated with the perturbations in the beam alone. Again, the rate at which energy is sapped from the beam is then Γ_W = 2ℑ(ω).