The Impact of Plasma Instabilities on the Spectra of TeV Blazars

TL;DR

This paper investigates how plasma instabilities affect the gamma-ray spectra of TeV blazars, potentially mimicking intergalactic magnetic field effects and influencing the interpretation of observational data.

Contribution

It provides a detailed analysis of various plasma instabilities and their impact on blazar spectra, offering insights into when these effects are significant or negligible.

Findings

Plasma instabilities can mimic IGMF effects in some blazar spectra.

For certain parameters, plasma instabilities have minimal impact on observed spectra.

The study helps refine methods to constrain intergalactic magnetic fields using gamma-ray observations.

Abstract

Relativistic jets from blazars are known to be sources of very-high-energy gamma rays (VHEGRs). During their propagation in the intergalactic space, VHEGRs interact with pervasive cosmological photon fields such as the extragalactic background light (EBL) and the cosmic microwave background (CMB), producing electron-positron pairs. These pairs can upscatter CMB/EBL photons to high energies via inverse Compton scattering, thereby continuing the cascade process. This is often used to set limits on intergalactic magnetic fields (IGMFs). However, the picture may change if plasma instabilities, arising due to the interaction of the pairs with the intergalactic medium (IGM), cool down the electrons/positrons faster than inverse Compton scattering. As a consequence, the kinetic energy lost by the pairs to the IGM could cause a hardening in the observed gamma-ray spectrum at energies below $\sim$100 GeV. Here we study several types and models of instabilities and assess their impact for interpreting observations of distant blazars. Our results suggest that plasma instabilities can describe the spectra of some blazars and mimic the effects of IGMFs, depending on parameters such as intrinsic spectrum of the object, the density and temperature of the IGM, and the luminosity of the beam. On the other hand, we find that for our fiducial set of parameters plasma instabilities do not have a major impact on the spectra of some of the blazars studied. Therefore, they may be used for constraining IGMFs.

Tables

Table 1. Table 1: List of the models used and the instability that dominates in each case. We also identify the treatment employed in the original references.

model	dominant instability	treatment	references
A	oblique	analytical	Broderick et al. (2012)
B	non-linear Landau damping	analytical	Miniati & Elyiv (2013)
C	modulation	analytical	Schlickeiser et al. (2012, 2013)
D	oblique	PIC	Sironi & Giannios (2014)
E	modulation	PIC	Vafin et al. (2018)

Table 2. Table 2: Tabulated values for 𝔗 B ​ ( D ) subscript 𝔗 B 𝐷 \mathfrak{T}_{\rm B}(D) from equation  17 .

${\log }_{10}\left(D/\mathrm{Mpc}\right)$	-0.05	0.14	0.35	0.59	0.79	0.96	1.17	1.40	1.57	1.77	1.99	2.20	2.41	2.60	2.80	3.00
${\log }_{10}{𝔗}_B\left(D\right)$	3.28	3.46	3.14	3.08	3.05	3.00	2.84	2.56	2.40	2.55	2.36	2.08	1.73	1.55	1.05	0.75

Table 3. Table 3: Objects used in this study, together with values for the parameters of their intrinsic spectrum that are compatible with observations, as taken from the corresponding references.

object	z	$\alpha$	$E_{\text{max}}$ [TeV]	observations	parameters
1ES 0229+200	0.140	1.2	5.0	Taylor et al. (2011)	Taylor et al. (2011)
1ES 0347-121	0.188	1.8	3.0	Bonnoli et al. (2015)	Costamante et al. (2018)
1ES 0414+009	0.287	1.9	2.0	Aliu et al. (2012)	Costamante et al. (2018)
1ES 1101-232	0.186	1.7	4.0	Abramowski et al. (2014)	Costamante et al. (2018)
1ES 1218+304	0.182	1.9	2.0	Taylor et al. (2011)	Costamante et al. (2018)
1ES 1312-423	0.105	1.9	1.3	HESS Collaboration et al. (2013)	HESS Collaboration et al. (2013)
RGB J0710+591	0.125	1.9	2.0	Taylor et al. (2011)	Costamante et al. (2018)
PKS 2155-304	0.116	1.9	10.0	Abramowski et al. (2014)	Abramowski et al. (2014)

Equations

$${\lambda }_{\text{pl}}\equiv \frac{2\pi c}{{\omega }_{\text{pl}}},$$

$${\omega }_{\mathrm{pl}}=\sqrt{\frac{e^2{n}_{\mathrm{IGM}}}{{ϵ}_0{m}_e}}$$

$$\mathcal{T}=\frac{2\gamma }{\sqrt{3n_{\mathrm{IGM}}}{\omega }_{\mathrm{pl}}}{\left(\frac{n_{\mathrm{IGM}}}{n_{\mathrm{beam}}}\right)}^{\frac{1}{3}},$$

$${\omega }_{\mathrm{pl}}=\sqrt{\frac{e^2{n}_{\mathrm{IGM}}}{{ϵ}_0{m}_e}},$$

$$1=\frac{n_{\mathrm{beam}}}{n_{\mathrm{IGM}}}\frac{m_{\mathrm{beam}}}{m_{\mathrm{IGM}}}\Psi \left(\Gamma ,\gamma ,\mathbf{k},\omega \right)-{\left(\omega -C\frac{kc\cos \theta n_{\mathrm{beam}}}{{\omega }_{\mathrm{pl}}{n}_{\mathrm{IGM}}}\right)}^{-1},$$

$${ϵ}_e=n_{\mathrm{beam}}{m}_e{c}^2(\gamma -1)$$

$${ϵ}_{\mathrm{crit}}=\frac{5}{3}\frac{n_{\mathrm{IGM}}{\left(k_{\mathrm{B}}{T}_{\mathrm{IGM}}\right)}^2}{m_e{c}^2}.$$

$$1=\frac{{\omega }_{\mathrm{pl}}^2}{{\omega }^2}+\frac{{\omega }^2{\sin }^2\theta \left(1-\frac{{\gamma }^2-1}{{\gamma }^2}{\cos }^2\theta \right)}{{\left(\omega -kc\cos \theta \sqrt{\frac{{\gamma }^2-1}{{\gamma }^2}}\right)}^2}.$$

$$\mathbf{B}=B\cos \left(ky\right){\mathbf{e}}_z.$$

$$\mathcal{T}=\sqrt{\frac{n_{\mathrm{beam}}{e}^2}{{ϵ}_0{m}_e\gamma }}.$$

$$\frac{\mathrm{d}{ϵ}_{\mathbf{k}}}{\mathrm{d}t}=2\frac{{ϵ}_{\mathbf{k}}}{{\mathcal{T}}_{\mathbf{k}}}-\frac{{ϵ}_{\mathbf{k}}{\omega }_{\mathrm{pl}}}{8(2\pi {)}^{5/2}{n}_e{m}_e{v}_e^2}\int \frac{(\mathbf{k}\cdot {\mathbf{k}}^{\prime })}{k^2{k}^{\prime 2}}\varphi (\mathbf{k},{\mathbf{k}}^{\prime }){ϵ}_{{\mathbf{k}}^{\prime }}\mathrm{d}{\mathbf{k}}^{\prime },$$

$$n_{\mathrm{IGM}}=n_{\mathrm{IGM},0}(1+z{)}^3,$$

$$n_{\mathrm{beam}}\simeq 3.7\times 10^{-22}\mathrm{c}{\mathrm{m}}^{-3}\left(\frac{\mathcal{L}}{10^{45}\text{erg}{\text{s}}^{-1}}\right)\left(\frac{E_e}{10^{12}\mathrm{eV}}\right){\left(\frac{1+z}{2}\right)}^{3\zeta -4},$$

$${\tau }_i={\mathcal{T}}_i,$$

$$n_{\mathrm{crit},\mathrm{A}}=1.6\times 10^{-19}{\left(\frac{E_e}{10^{12}\mathrm{eV}}\right)}^{-2}\left(\frac{n_{\mathrm{IGM}}}{10^{-7}\mathrm{c}{\mathrm{m}}^{-3}}\right)\mathrm{c}{\mathrm{m}}^{-3}.$$

$${\tau }_{\mathrm{A}}=\{\begin{array}{l}\begin{array}{rl}& 7.1\times 10^7\mathrm{s}{\left(\frac{E_e}{10^{12}\mathrm{eV}}\right)}^{-1}{\left(\frac{n_{\mathrm{beam}}}{10^{-22}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{n_{\mathrm{IGM}}}{10^{-7}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{+\frac{1}{2}}\\ & \text{if}{n}_{\mathrm{beam}}\le n_{\mathrm{crit},\mathrm{A}}\text{(warm-plasma beam)},\end{array}\\ \begin{array}{rl}& 5.1\times 10^5\mathrm{s}{\left(\frac{E_e}{10^{12}\mathrm{eV}}\right)}^{\frac{1}{3}}{\left(\frac{n_{\mathrm{beam}}}{10^{-22}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-\frac{1}{3}}{\left(\frac{n_{\mathrm{IGM}}}{10^{-7}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-\frac{1}{6}}\\ & \text{if}{n}_{\mathrm{beam}}>n_{\mathrm{crit},\mathrm{A}}\text{(cold-plasma beam)},\end{array}\end{array}$$

$${\tau }_{\mathrm{B}}={\tau }_{\mathrm{IC}}{\mathfrak{T}}_{\mathrm{B}}(D)(1+z{)}^2=3.9\times 10^{13}\mathrm{s}(1+z{)}^{-2}{\left(\frac{E_e}{10^{12}\mathrm{eV}}\right)}^{-1}{\mathfrak{T}}_{\mathrm{B}}(D),$$

$${\tau }_{\mathrm{C}}=\{\begin{array}{l}\begin{array}{rl}& 5.2\times 10^{14}\mathrm{s}{\left(\frac{E_e}{10^{12}\mathrm{eV}}\right)}^{\frac{5}{3}}{\left(\frac{n_{\mathrm{beam}}}{10^{-22}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{\frac{1}{3}}{\left(\frac{n_{\mathrm{IGM}}}{10^{-7}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-\frac{5}{6}}\\ & \times {\left(\frac{T_{\mathrm{IGM}}}{10^4\mathrm{K}}\right)}^{-2}\text{if}{n}_{\mathrm{beam}}\le n_{\mathrm{crit},\mathrm{C}}\text{(weak blazar)},\end{array}\\ \begin{array}{rl}& 8.3\times 10^6\mathrm{s}{\left(\frac{E_e}{10^{12}\mathrm{eV}}\right)}^{\frac{1}{3}}{\left(\frac{n_{\mathrm{beam}}}{10^{-22}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-\frac{1}{3}}{\left(\frac{n_{\mathrm{IGM}}}{10^{-7}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-\frac{1}{6}}\\ & \times {\mathfrak{T}}_{\mathrm{C}}(n_{\mathrm{IGM}},T_{\mathrm{IGM}})\text{if}{n}_{\mathrm{beam}}>n_{\mathrm{crit},\mathrm{C}}\text{(strong blazar)},\end{array}\end{array}$$

$${\mathfrak{T}}_{\mathrm{C}}(n_{\mathrm{IGM}},T_{\mathrm{IGM}})=1+\frac{5}{4}\ln \left(\frac{T_{\mathrm{IGM}}}{10^4\mathrm{K}}\right)-\frac{1}{4}\ln \left(\frac{n_{\mathrm{IGM}}}{10^7\mathrm{c}{\mathrm{m}}^{-3}}\right)$$

$$n_{\mathrm{crit},\mathrm{C}}=2.5\times 10^{-25}{\left(\frac{E_e}{10^{12}\mathrm{eV}}\right)}^{-1}\left(\frac{n_{\mathrm{IGM}}}{10^{-7}\mathrm{c}{\mathrm{m}}^{-3}}\right){\left(\frac{T}{10^4\mathrm{K}}\right)}^2\mathrm{c}{\mathrm{m}}^{-3}.$$

$${\tau }_{\mathrm{D}}=\{\begin{array}{l}\begin{array}{rl}& 1.4\times 10^7\mathrm{s}{\left(\frac{E_e}{10^{12}\mathrm{eV}}\right)}^{-1}{\left(\frac{n_{\mathrm{beam}}}{10^{-22}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{n_{\mathrm{IGM}}}{10^{-7}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{+\frac{1}{2}}\\ & \text{if}{n}_{\mathrm{beam}}\le n_{\mathrm{crit},\mathrm{D}}\text{(warm-plasma beam)},\end{array}\\ 9.6\times 10^5\mathrm{s}{\left(\frac{E_e}{10^{12}\mathrm{eV}}\right)}^{\frac{1}{3}}{\left(\frac{n_{\mathrm{beam}}}{10^{-22}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-\frac{1}{3}}{\left(\frac{n_{\mathrm{IGM}}}{10^{-7}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-\frac{1}{6}}\\ \text{if}{n}_{\mathrm{beam}}>n_{\mathrm{crit},\mathrm{D}}\text{(cold-plasma beam)},\end{array}$$

$$n_{\mathrm{crit},\mathrm{D}}=8.0\times 10^{-20}{\left(\frac{E_e}{10^{12}\mathrm{eV}}\right)}^{-2}\left(\frac{n_{\mathrm{IGM}}}{10^{-7}\mathrm{c}{\mathrm{m}}^{-3}}\right)\mathrm{c}{\mathrm{m}}^{-3}.$$

$$\begin{array}{rl}{\tau }_{\mathrm{E}}& =1.9\times 10^{11}\mathrm{s}{\left(\frac{E_e}{10^{12}\mathrm{eV}}\right)}^{\frac{4}{3}}{\left(\frac{n_{\mathrm{beam}}}{10^{-22}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{\frac{1}{3}}{\left(\frac{n_{\mathrm{IGM}}}{10^{-7}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-\frac{1}{3}}\\ & \times {\left(\frac{T_{\mathrm{IGM}}}{10^4\mathrm{K}}\right)}^{-1}.\end{array}$$

$${\lambda }_{\text{PP}}\simeq 80\mathrm{Mpc}\frac{\kappa }{(1+z{)}^2}{\left(\frac{E_{\gamma }}{10^{13}\mathrm{eV}}\right)}^{-1},$$

$${\lambda }_{\text{IC}}\simeq 1.2\mathrm{kpc}(1+z{)}^{-3}.$$

$$P(E_e,x,z)\equiv -\frac{\mathrm{d}{E}_e}{\mathrm{d}x}(E_e,x,z)=\frac{E_e}{c\tau (E_e,x,z)},$$

$$\frac{dN}{dE}\propto E^{-\alpha }\exp \left(-\frac{E}{E_{\text{max}}}\right).$$

Full text

††pubyear: 2019††pagerange: The Impact of Plasma Instabilities on the Spectra of TeV Blazars–The Impact of Plasma Instabilities on the Spectra of TeV Blazars††volume: 489††pagerange: 3836-3849

1 Introduction

High-energy (i.e. GeV and TeV) and very-high-energy gamma rays (VHEGRs) have become important tools to probe astrophysical and cosmological phenomena. This includes studies of our own galaxy (Weekes et al., 1989; Strong et al., 2004; Aharonian et al., 2006b), BL Lac objects (Quinn et al., 1996), and gamma-ray bursts (Gehrels & Meszaros, 2012). Cosmological scales can also be probed with high-energy gamma rays, enabling constraints on the opacity of the Universe (Taylor et al., 2011), and heating of the intergalactic medium (IGM) (Chang et al., 2012), among others. The plethora of information that can be obtained with this messenger is possible mainly thanks to experimental efforts such as the High-Energy Stereoscopic System (H.E.S.S.) (Hofmann, 2000), the Very Energetic Radiation Imaging Telescope Array System (VERITAS) (Weekes et al., 2002), the Fermi Large Area Telescope (Fermi-LAT) (Atwood et al., 2009), and the Major Atmospheric Gamma Imaging Cherenkov (MAGIC) (Rico, 2017). Future air-Cherenkov telescopes such as the upcoming Cherenkov Telescope Array (CTA) (CTA Consortium, 2019) will provide further insights into the high-energy universe.

Of particular interest are gamma-ray observations of TeV blazars, i.e. active galactic nuclei (AGNs) whose jets are nearly parallel to the line of sight of emission. The standard model for the propagation of TeV gamma rays from the source to the observer is well-understood (Neronov & Semikoz, 2009). The emitted VHEGRs interact with background photons from pervasive radiation fields, predominantly the extragalactic background light (EBL), producing electron-positron pairs. The most relevant component of EBL in the context of this work is the infrared light, which consists mostly of diffuse photons stemming from star-formation processes, and starlight reprocessed by intergalactic dust (Hauser & Dwek, 2001). It is rather challenging to model the EBL (Fan et al., 2004), especially at high redshifts. For this reason, a number of different models exists (Kneiske et al., 2004; Stecker et al., 2006; Franceschini et al., 2008; Kneiske & Dole, 2010; Finke et al., 2010; Domínguez et al., 2011; Gilmore et al., 2012; Stecker et al., 2016). This, in turn, results in different estimates of the mean free path of photons. In addition, higher order processes such as double pair production (Heiter et al., 2018) and interactions with other pervasive photon fields like the radio background (Protheroe & Biermann, 1996) are possible. The pair-produced electrons (and positrons) can interact with low-energy photons from the cosmic microwave background (CMB), mostly via inverse Compton (IC) scattering. Here, the high-energy electrons upscatter the low-energy photons to energies in the GeV–TeV range. The mean free path of these electrons is typically $\lesssim 100\;\text{kpc}$, much less than that for pair production ($\sim 100\;\text{Mpc}$ at 10 TeV). Again, higher order processes like triplet pair production, even though highly suppressed, may be important for photons with energies $\gtrsim 10^{15}\;\text{eV}$ (Heiter et al., 2018).

Intergalactic magnetic fields (IGMFs) can deflect the charged particles produced in the cascade, i.e., electrons and positrons. Hence, the observed arrival direction of the secondary photons produced does not necessarily point back to the source. The result is that, instead of a point-like source (as would be the case in the absence of magnetic fields), a halo-like structure emerges (Aharonian et al., 1994; Neronov & Semikoz, 2007, 2009; Dolag et al., 2009; Ackermann et al., 2018; Chen et al., 2018). The size of the halo is determined by an interplay between the magnetic-field strength ($B$), its coherence length ($L_{B}$), and the electron energy. In general, $B$ relates to the deflection angle of electrons via the Lorentz force, whereas $L_{B}$ determines the propagation regime (diffusive when the Larmor radii of electrons is much smaller than the $L_{B}$, and ballistic otherwise). In addition, it has been shown that the shape of the halo is important as well as it may contain information about the magnetic helicity of IGMFs (Tashiro et al., 2014; Alves Batista et al., 2016a; Duplessis & Vachaspati, 2017), since helical fields produce spiral-like structures in the arrival directions.

The second effect of IGMFs on gamma-ray observables are time delays (Plaga, 1995; Fan et al., 2004; Murase et al., 2008; Neronov & Semikoz, 2009; Fitoussi et al., 2017). The effective propagation time of a secondary gamma ray which is observed as the result of the cascading processes described above is longer than the propagation time of a primary photon, i.e. a gamma ray which arrives at Earth without interacting with the EBL. Therefore, if two such photons are emitted at the same time in a flare, they may be observed at different times, i.e. the "reprocessed" photon arriving with a time delay with respect to the arrival time of the primary one.

Finally, another consequence of the propagation of VHEGRs in IGMFs is the suppression of the photon flux in the GeV range. This is a direct consequence of the energy dependence of the deflection angle. This angle increases dramatically with decreasing energy, such that GeV photons, indeed, are deflected away from the line of sight in such a way that they can not be observed at Earth, thus resulting in the aforementioned flux suppression. So far, this method has been used to derive only the lower limit on the IGMF strength (Dai et al., 2002; d’Avezac et al., 2007; Neronov & Vovk, 2010; Taylor et al., 2011).

While seemingly simple, the underlying mechanism of electromagnetic cascades is controversial. As first pointed out by Broderick et al. (2012), the pairs may interact with the IGM due to emerging plasma instabilities in the beam, cooling down before IC scattering can occur. This is a possible explanation for the absence of photons at GeV energies, observed in the spectra of some blazars. Yet, it should be noted that this interpretation, too, is disputed. While it is true that instabilities should occur when a relativistic beam traverses a background plasma, it is presently virtually impossible to directly simulate these scenarios, such that one has to rely on extrapolations which, however, are rather uncertain.

In this paper we study the spectra of some specific blazars for different plasma instability models. We first (section 2) provide a brief description of the types of instabilities considered, including a discussion of the relevant physical parameters and conditions, accompanied by a number of models found in the literature. Then, in section 3, we describe the implementation of the instabilities in our Monte Carlo code. In section 4 we study the effect of some relevant quantities on the expected gamma-ray fluxes for one specific scenario, while section 5 gives an account of the predicted fluxes of gamma rays for a selected list of blazars. In section 6 we discuss our results, and present our concluding remarks in section 7.

2 Plasma Instabilities in the IGM

Collective phenomena start to become relevant for high enough densities of electron-positron pairs. More precisely, the number of pairs within a spherical volume of radius equal to the wavelength of the most unstable mode – the plasma skin depth – which reads

[TABLE]

where

[TABLE]

gives the plasma frequency, and $n_{\text{IGM}}$ denotes the density of free-electrons in the IGM. In other words, plasma effects are relevant if the number of pairs in the beam within a volume of radius $\lambda_{\text{pl}}$ is much larger than unity, i.e., $\lambda_{\text{pl}}^{3}n_{\text{beam}}\gg 1$, wherein $n_{\text{beam}}$ is the number density of pairs in the beam.

For the densities involved in blazar beams interacting with the IGM, $n_{\text{beam}}/n_{\text{IGM}}\ll 1$, thus ensuring that plasma effects are relevant for the evolution of the beams.

In this section we present some instabilities that might occur due to interactions between the beam and the background plasma. In principle, all of them can be relevant for the interaction between ultrarelativistic particles and the IGM.

2.1 Types of instabilities

The two-stream instability, first described in Bohm & Gross (1949), appears if the wave vector of an electrostatic perturbation (a Langmuir wave) is parallel to the flow for a range of wave numbers $0<k\leq\omega_{\rm pl}/v_{\rm beam}$.

An intuitive picture of this instability can be obtained considering two flows moving towards each other. The lower-density flow – the beam – interacts with the higher-density background – the plasma – triggering the instability if the thermal velocity of the particles is smaller than the drift between the beams. The Langmuir waves arising in the background plasma resonate with those in the beam, leading to the growth of the instability.

In a warm111A beam is considered to be warm if the condition $T\gtrsim\frac{3}{2^{10/3}}\left(\frac{n_{\rm beam}}{n_{\rm IGM}}\right)^{\frac{2}{3}}\gamma^{\frac{1}{3}}\frac{m_{e}c^{2}}{k_{\rm B}}$ is fulfilled, according to Broderick et al. (2012). beam the velocity ($v_{\text{beam}}$) of some particles equals the phase velocity of the wave, hence being in resonance with it. Due to the form of the velocity distribution, more energy is transferred from the beam particles to the wave mode than vice-versa, thus resulting in an instability. For a cold-plasma beam the mechanism is slightly different, but still results in growing wave modes since the anisotropic electron and positron distribution functions interact with the co-moving background electrostatic wave, when the wave vector is parallel to the beam momentum, resulting in an instability growth time of (Broderick et al., 2012)

[TABLE]

where $n_{\rm IGM}$ and $n_{\rm beam}$ are the number densities of the IGM and the beam, respectively, $\gamma$ is the Lorentz factor of the particle and $\omega_{\rm pl}$ is the plasma frequency given by

[TABLE]

with $\epsilon_{0}$ being the electrical vacuum permittivity, and $e$ and $m_{e}$ are the electron charge and mass, respectively.

The oblique instability is, in principle, very similar to the two-stream instability described above. The main difference is the fact that in this case the wave vector is not parallel to the beam velocity; instead, the two are at an angle $\theta$ to each other. It can then be shown (Nakar et al., 2011) that in the astrophysical setting the dispersion relation is given by

[TABLE]

where $\mathbf{k}$ is the wave vector, $C$ is the fractional charge of the beam, $\Gamma$ is the Lorentz factor and $\Psi$ is a function of the given quantities. For $\theta=0$ this equation can be solved analytically and gives the growth time for the two-stream instability discussed above. To find the growth time for the oblique instability at an arbitrary angle $\theta$ one has, in general, to resort to numerical simulations. As pointed out by Broderick et al. (2012), these two growth rates can differ significantly due to the fact that the relevant quantity for the beam–wave interaction is the projected velocity of the beam onto the wave vector, which for large $\theta$ is rather small, such that it is easier to create deflections and, subsequently, instabilities.

The modulation instability in its simplest form can be explained as the result of ions in a turbulent medium scattering off oscillations caused by the beam. For this reason, they are transferred to smaller wave numbers, shifting the wave energy to higher (even superluminal) phase speeds (Galeev et al., 1977). This takes place if the energy density of the electrostatic fluctuations $\epsilon_{e}$, given by (Schlickeiser et al., 2002)

[TABLE]

becomes larger than the critical value (Galeev et al., 1975; Schlickeiser et al., 2013)

[TABLE]

If this condition is fulfilled, one can find the growth rate for the modulation instability by analyzing the dispersion relation

[TABLE]

It should be noted here that formally the modulation instability is nothing but an oscillating two-stream instability (Papadopoulos, 1975; Schamel et al., 1976; Schlickeiser et al., 2012).

The Weibel instability, named after the author of the seminal work on the topic, Weibel (1959), is the result of the interaction of several counter-streaming beams (Fried, 1959). In its simplest form, following Medvedev & Loeb (1999), we can consider two beams in opposite directions with a vanishing net current. Small anisotropies in velocity cause fluctuations of the magnetic field perpendicular to the direction of the beam. For example, without loss of generality, if we take the particles to move along the $x$-axis, then the magnetic fluctuation can be represented in the form

[TABLE]

The resulting Lorentz force deflects the electrons in the beams in such a way that two current sheets are produced. Such large-scale currents induce larger magnetic fields which, in turn, increase the force, thus enhancing the effect and consequently producing the instability. The instability growth time is then given by (Broderick et al., 2012)

[TABLE]

Note that while the Weibel instability does not contain oscillatory modes, they are present in the two-stream instability, which represents a fundamental difference between the two.

The Weibel instability grows until the magnetic field is amplified to a point where the Larmor radii of the particles become comparable to the plasma skin depth, when the growth saturates (Davidson et al., 1972). Moreover, small temperature gradients perpendicular to the beam direction can kill off this instability, sometimes giving it a short lifetime (Broderick et al., 2012).

Finally, another possible mechanism is non-linear Landau damping. As the name suggests, the effect is highly non-linear and may be described as induced scattering by thermal ions, such that the frequency and wave vector of a given plasma (Langmuir) wave are transformed to different values (Brejzman & Ryutov, 1974; Chang et al., 2014). The corresponding kinetic equation, which describes the wave transformation $(\mathbf{k},\omega)\rightarrow(\mathbf{k}^{\prime},\omega^{\prime})$, is given by

[TABLE]

where $\epsilon_{\mathbf{k}}$ is the spectral energy density and $\phi(\mathbf{k},\mathbf{k}^{\prime})$ is the overlap integral given in Kaplan & Tsytovich (1968). Here, we have $1/\mathcal{T}_{\mathbf{k}}=1/\mathcal{T}_{\rm inst}+1/\mathcal{T}_{\rm LD}$, where $\mathcal{T}_{\rm inst}$ denotes the instability growth time, and $\mathcal{T}_{\rm LD}$ denotes the growth time for linear Landau damping.

2.2 Physical Parameters

In this section we discuss the IGM and beam parameters which have an influence on the growth of plasma instabilities.

First, one should consider the number density of the IGM ($n_{\rm IGM}$), which evolves as

[TABLE]

wherein $z$ is the redshift at which the number density is calculated, and $n_{\text{IGM,0}}$ is the IGM density at present time. Our fiducial value for $n_{\rm IGM,0}$ is $n_{\text{IGM,0}}=10^{-7}\;\text{cm}^{-3}$ (Broderick et al., 2012). It should be noted that this is a rather simple model, not taking into account density fluctuations; these may be included if we multiply equation 12 by $(1+\delta)$, where $\delta$ is the overdensity.

The temperature of the IGM ($T_{\rm IGM}$) has many uncertainties associated with it. Depending on whether one considers cosmic voids or the intracluster medium, the value of $T_{\rm IGM}$ ranges from $\sim 10^{3}\,{\rm K}$ (Hui & Gnedin, 1997) in voids, to $\sim 10^{9}\,{\rm K}$ in clusters (Reimer & Böttcher, 2013; Perucho et al., 2011; Burns et al., 2010). In order to account for this range, we adopt $T_{\rm IGM}=[10^{3},10^{4},10^{5},10^{6}]$, with $T_{\rm IGM}=10^{4}\,{\rm K}$ being our fiducial value for the simulations presented below. This choice is fully justified if we consider that the volume filling factor of voids exceeds the one of clusters by orders of magnitude, so that the beam will interact predominantly with voids, hence the choice of temperatures. This range also encompasses the values considered by Broderick et al. (2012), Miniati & Elyiv (2013), Schlickeiser et al. (2012), Sironi & Giannios (2014), and Vafin et al. (2018).

When it comes to the source parameters, the first important one is the total isotropic-equivalent luminosity ($\mathcal{L}$). For TeV sources this value lies in the range $10^{41}\;\text{erg}\,\text{s}^{-1}$ and $10^{47}\;\text{erg}\,\text{s}^{-1}$ (Broderick et al., 2012). We adopt the fiducial value of $\mathcal{L}=10^{45}\;\text{erg}\,\text{s}^{-1}$.

Another important quantity for the calculation of plasma instabilities is the beam plasma number density ($n_{\rm beam}$). As this is a dynamic quantity, it is usually only possible to estimate it. In this work we follow Broderick et al. (2012), using an upper limit for $n_{\rm beam}$ which comes from analytical estimates of the cascade development assuming IC cooling and that the cooling is dominated by the kinetic oblique mode:

[TABLE]

wherein $\zeta=4.5$ for $z<1$ is a parameter that can be inferred from the analysis of the local star formation rate (Kneiske et al., 2004).

2.3 Models

Here we present and discuss the different models which we use to obtain blazar spectra in chronological order. In general, there are two different fundamental time scales to be considered. The first is the instability growth rate, which we label $\mathcal{T}_{i}$. The second is the energy loss time of the electrons/positrons due to plasma effects, $\tau_{i}$, where in both cases $i$ refers to the model considered (see below). For the calculation of the blazar spectra the relevant quantity is the electron energy loss time $\tau_{i}$. In order to obtain the most robust lower limit for the photon flux, we adopt the estimate

[TABLE]

even though in principle $\tau_{i}$ might be substantially larger than $\mathcal{T}_{i}$ (Grognard, 1975; Pavan et al., 2011).

Note that while any of the instabilities described above can operate, for each of the models studied here there is always a dominant one. We summarise this in table 1. Further details are discussed below.

2.3.1 Model A

One of the first works on this topic and the one on which we mostly base our model A, is Broderick et al. (2012), later on expanded and further analyzed by Chang et al. (2012); Chang et al. (2014, 2016); Broderick et al. (2014); Shalaby et al. (2017, 2018). The authors distinguish between two regimes referred to as the ‘warm’ and ‘cold’ plasma beam cases depending on whether the beam density, $n_{\rm beam}$, is below (warm) or above (cold) the critical value $n_{\rm crit,A}$:

[TABLE]

Broderick et al. (2012) come to the conclusion that in both cases the oblique instability is the dominant one, however resulting in different expressions which can be summarized as

[TABLE]

For typical values ($n_{\text{IGM}}\sim 10^{-7}\;\text{cm}^{-3}$ and $n_{\text{beam}}\sim 10^{-22}\;\text{cm}^{-3}$) we find ourselves in a regime where the second expression from equation 16 holds, assuming that the electrons have energies $E_{e}\sim\;\text{a few TeV}$. If the electrons have higher energies ($E_{e}\sim 10\;\text{TeV}$), then in overdense regions of the IGM $n_{\text{beam}}\gtrsim n_{\text{crit,A}}$, so that $\tau_{A}$ is given by the first expression of eq. 16.

2.3.2 Model B

Miniati & Elyiv (2013) consider the interplay between Langmuir waves and non-linear Landau damping to be the most relevant effect. They have taken into account the finite nature of the spread of the momentum component transverse to the beam, which were obtained from Monte Carlo simulations of the cascade development. They conclude that due to that the effect of plasma oscillations on energy losses of the electron-positron beams is mostly negligible, since inhomogeneities in the IGM can break the resonance between modes of the Langmuir and the beam waves. Nevertheless, for the sake of completeness, we consider it here and compare it with the other models.

The authors obtained the energy loss time, $\tau_{\rm B}$, through a combination of Monte Carlo simulations to determine the beam properties as it propagates, and analytic solutions of appropriate equations, to obtain the instability growth rate. We can write

[TABLE]

where $E_{e}$ is the electron/positron energy, $D$ is the co-moving distance to the source, and $\mathfrak{T}_{\rm B}(D)$ is tabulated as shown in table 2.

From table 2 one can see that non-linear Landau damping is more important closer to the blazar. On the other hand, the inhomogeneities in the IGM are more relevant at larger distances. The interplay between these two effects lead to an effective stabilisation of them beam, thus explaining why the instabilities are completely subdominant in this scenario.

2.3.3 Model C

This model is based on analytic calculations developed by Schlickeiser et al. (2012, 2013). In particular, in Schlickeiser et al. (2012) the authors find that the different effects contribute to the suppression of the electromagnetic cascade at different regimes. In the ‘strong blazar’ regime, the beam plasma density is lower than a certain critical value $n_{\rm crit,C}$, resulting in a dominance of the modulation instability. If $n_{\rm beam}>n_{\rm crit,C}$ – the ‘weak blazar’ case – the modulation instability does not set in. In this case, non-Linear Landau damping becomes relevant, such that energy is deposited in electrostatic and electromagnetic fluctuations in the background plasma.

The energy loss time is

[TABLE]

where $\mathfrak{T}_{\rm C}$ is

[TABLE]

and $n_{\rm crit,C}$ is given by

[TABLE]

Note that while the cooling rate of electrons depends quadratically on the temperature in the weak blazar case, in the strong blazar case this dependence is much weaker. This is a consequence of the onset of the modulation instability, which depends on both the beam density, $n_{\text{beam}}$, as well as on the IGMF temperature, $n_{\text{IGM}}$.

2.3.4 Model D

This model is based on Sironi & Giannios (2014). The authors distinguish two cases, the ‘cold-plasma beam’ and the ‘warm-plasma beam’, depending on whether the beam plasma density $n_{\rm beam}$ is above or below value $n_{\rm crit,D}$, respectively. For both cases they find that the oblique instability is the most relevant one, resulting in the energy loss time

[TABLE]

where $n_{\rm crit,D}$ is given by

[TABLE]

At first glance, this model is remarkably similar to model A, apart from some numerical factors (compare equations 16 and 21). However, these two models are intrinsically different. While Broderick et al. (2012) argued that the background plasma takes about $\sim 50\%$ of the beam energy as the beam relaxation time is much shorter than IC losses, the PIC simulations of Sironi & Giannios (2014) result in an energy-transfer factor of $\sim 10^{-1}$. This discrepancy stems, in part, from the treatment of the transverse momentum, which inherits from the momentum distribution of the pairs generated in the cascade.

Sironi & Giannios (2014) also note that the oblique exponential growth takes place on timescales much shorter than the time it takes for the relaxation phase to be reached. For this reason, there is an associated efficiency factor that indicates how much of the beam’s energy is used to heat the IGM, which depends on the spectrum of the component of the electrons’ momenta perpendicular to the beam. Here we take this factor to be 1, since we are interested in investigating the scenario that provides the maximum quenching of the cascade.

2.3.5 Model E

The final model we adopt is mostly taken from Vafin et al. (2018), who find that the modulation instability dominates the energy losses. Through a thorough analysis, they arrive at

[TABLE]

Vafin et al. (2018) claim that non-linear Landau damping is properly resolved in their PIC simulations; in this case, the peak frequency is $\sim 10^{-5}\omega_{\text{pl}}$, whereas the peak frequency for the linear growth is $\sim 10^{-3}\omega_{\text{pl}}$.

The electrostatic mode in stabilised by the modulation instability when the energy density stored in the waves is a fraction of that in the beam. This is related to an efficiency factor which, once again, we have assumed to be maximal for the purposes of this study.

3 Simulation Setup

For this analysis we employ the Monte Carlo code CRPropa 3 (Alves Batista et al., 2016b). It enables the propagation of high-energy particles in the intergalactic space. It includes the relevant interactions for the development of electromagnetic cascades, namely inverse Compton scattering, pair production, and adiabatic losses due to the expansion of the Universe. We consider the CMB and the EBL as target background fields and choose the EBL model by Gilmore et al. (2012) for the latter. This choice quantitatively affects the shape of the simulated spectra, but it is virtually irrelevant for the qualitative discussion that follows.

The most important background for pair production at the energies of interest ($E\lesssim 10^{14}\;\text{eV}$) is the EBL, whose interaction rates are shown in figure 1 for a few different models. Although this figure provides a detailed model for pair production, a rough analytical estimate of the mean free path for this process yields (Neronov & Semikoz, 2009)

[TABLE]

where $E_{\gamma}$ is the photon energy, $z$ is the redshift and $\kappa$ is an energy- and redshift-dependent coefficient to account for uncertainties in the EBL.

For electrons, the dominant background is the CMB. The energy-loss length, in this case, is nicely approximated by (Neronov & Semikoz, 2009)

[TABLE]

This approximate analytical expression agrees exceedingly well with what we obtain with our detailed Monte Carlo treatment. Note that $\lambda_{\text{IC}}\ll\lambda_{\text{PP}}$, implying that the charged component of the cascade does not live very long compared to photons.

We took advantage of the modular structure of CRPropa code to develop an extension222https://github.com/rafaelab/grplinst dedicated to calculate the effects of plasma instabilities in an approximate fashion. To this end, we convert the energy loss times given by equations 16, 17, 18, 21, and 23 to the energy loss per distance, $P(E,x,z)$, which in the ultrarelativistic limit can be written as

[TABLE]

where $x$ is the trajectory length propagated by a particle and $c$ is the speed of light.

We should stress that the values of $\tau$ for models B, C, and E are the actual energy loss times of the beams, as in our treatment; for models A and D, they correspond to the maximum instability growth time. As has been argued in the corresponding publications, the actual energy loss (and hence the value of $P$) can be significantly smaller. Therefore, all values of $\tau$ presented here should be interpreted as lower limits; this should be kept in mind for the analyses carried out in sections 4 and 5.

To calculate the energy loss function $P$ for the different models, we combine the respective values for $\tau$ from section 2.3 with the model of the beam density $n_{\rm beam}$ from equation 13 (with $\zeta=4.5$), as we are only interested in blazars with $z\ll 1.0$ and the value of $n_{\rm IGM}$ from equation 12.

In figure 2 we compare the electron cooling rates for the different models, for a typical combination of parameters. We also show the interaction rate for inverse Compton scattering. Note that the energy dependence across the different models differs considerably.

As a benchmark scenario, we consider the case of a source located at $z\approx 0.14$, which is approximately the redshift of the extreme blazar 1ES 0229+200, commonly used in studies of IGMFs (Neronov & Vovk, 2010; Vovk et al., 2012; Arlen et al., 2014; Yang & Dai, 2015; Yan et al., 2019). We model the spectrum injected by the blazar as a power law of the form

[TABLE]

Unless stated otherwise, we adopt the parameters from Vovk et al. (2012): $\alpha=1.2$ and $E_{\text{max}}=5\times 10^{12}\;\text{eV}$. We use the EBL model by Gilmore et al. (2012), although this choice should not significantly change the interpretation of the results.

4 Simulation results

Now we present the results of the simulations. Throughout this section, we present the results for our benchmark scenario, fixing $T_{\text{IGM},0}=10^{4}\;\text{K}$, $n_{\text{IGM},0}=10^{-7}\;\text{cm}^{-3}$, and $\mathcal{L}=10^{45}\;\text{erg/s}$ and varying one of the parameters at a time. As stated in section 2.2, we also fix $\alpha=1.2$ and $E_{\text{max}}=5\;\text{TeV}$.

4.1 Plasma instability models

We first compare all the models for our benchmark scenario. This is shown in figure 3. Note that the slope for all models but B are rather similar in the energy range of $\sim 10-300\;\text{GeV}$. This is expected because the spectral index of injection ($\alpha=1.2$ for the benchmark scenario) should be retrieved if virtually all secondary electrons are removed from the cascade. If plasma instabilities are not as efficient as we considered, then the picture would change and a fraction of the electrons would remain in the cascade. If identified, this slope would be a phenomenological signature of a very efficient process quenching the cascade; plasma instabilities would then be a candidate explanation. Departures from this behaviour may occur for different values of $n_{\text{IGM},0}$, $T_{\text{IGM},0}$, and $\mathcal{L}$, as will be discussed later. Nevertheless, for $30\lesssim E/\text{GeV}\lesssim 100$ all models A, C, D, and E led to a slope of $E^{-1.2}$ at energies $\lesssim 100\;\text{GeV}$ in most cases studied, the exception being model C for some particular combinations of parameters.

We confirm the results by Miniati & Elyiv (2013), who argue that plasma effects do not lead to noticeable changes in the observed gamma-ray flux. This is shown in figure 3 and can be confirmed by analysing figure 2, which shows that cooling rates for model B are orders of magnitude below the interaction rate for IC scattering at the energy range of interest ($E\lesssim 1\;\text{PeV}$). For this reason, we omit plots for model B when studying the effects of temperature, IGM density, beam luminosity, and spectral parameters.

4.2 Temperature of the IGM

We study the effects of the temperature on the spectrum for our benchmark scenario. We vary only the temperature, while keeping all other quantities ($n_{\text{IGM},0}=10^{-7}\;\text{cm}^{-3}$ and $\mathcal{L}=10^{45}\;\text{erg/s}$) fixed.

As shown in figure 4, the temperature dependence is more prominent for models C and E, while models A and D are virtually temperature-independent. This follows immediately from equations 16, 18, 21, and 23 as there are no explicit temperature dependencies in equations 16 and 21.

At $T_{\text{IGM,0}}\sim 10^{4}\;\text{K}$ the spectrum for model C presents a change in slope at energies $\lesssim 3\;\text{GeV}$. This may be an observational signature of this model, and if observed in multiple blazars, would allow us to draw conclusions about the temperature or correlated parameters.

4.3 Density of the IGM

Now we discuss how the density of the IGM affects the gamma-ray spectrum. We fix $T_{\text{IGM},0}=10^{4}\;\text{K}$ and $\mathcal{L}=10^{45}\;\text{erg/s}$, varying the density only. The predicted gamma-ray spectra for our benchmark object are shown in figure 5.

Models A, D, and E are unremarkable, in the sense that the cascade is completely suppressed by the plasma instabilities, as shown in figure 5. Model C presents clear spectral signatures of the IGM density. There is a bump at $E\sim 5\;\text{GeV}$ for $n_{\text{IGM},0}=10^{-7}\;\text{cm}^{-3}$. The onset of this feature seems to start at higher energies as $n_{\text{IGM},0}$ increases. With our simulations we cannot confirm if this bump is shifted to lower energies, as our minimum energy333The lower energy cutoff of 1 GeV was chosen to optimise the time required to run the simulations, while covering a range of energies that would allow for meaningful conclusions to be drawn. is set to 1 GeV. However, if we extrapolate the analysis of the scenario with $n_{\text{IGM},0}=10^{-6}\;\text{cm}^{-3}$, it is not unreasonable to expect this to be the case.

4.4 Beam luminosity

We proceed our study of the relevant parameters by investigating how the luminosity of the blazar beam affects the shape of the gamma-ray spectrum. Once again, we fix all other variables and vary solely the luminosity for $T_{\text{IGM},0}=10^{4}\;\text{K}$ and $n_{\text{IGM},0}=10^{-7}\;\text{cm}^{-3}$. The results are shown in figure 6.

As in the case of the temperature and the IGM density, spectral features arise for model C and, to a lesser degree, for models A and D. As the beam luminosity decreases, an increase in the flux at $E\sim 3\;\text{GeV}$ can be seen. The onset of this flux enhancement with respect to the full suppression of the cascade becomes visible at increasingly higher energies as the beam becomes dimmer.

The aforementioned spectral signatures are not surprising. In fact, equations 16, 18, 21, and 23 all depend on the beam density, as can be seen in equation 13.

The vast majority of blazars have isotropic-equivalent luminosities $\mathcal{L}\lesssim 10^{45.5}\;\text{erg/s}$. Therefore, the cooling due to plasma instabilities may not be as severe as for the case of very luminous objects.

4.5 Spectral parameters

We now investigate the effects of the spectral parameters ($E_{\text{max}}$ and $\alpha$) on the shape of the predicted gamma-ray spectra. To this end, we use our benchmark scenario fixing the following parameters: $T_{\text{IGM},0}=10^{4}\;\text{K}$, $n_{\text{IGM},0}=10^{-7}\;\text{cm}^{-3}$, and $\mathcal{L}=10^{45}\;\text{erg/s}$.

As previously mentioned, for model B there are no significant spectral changes. For the other models the spectral modifications are similar to each other, so that we choose only one model to illustrate this discussion: model A.

The energy-dependent quenching factor is defined as the ratio between $j_{\text{pl}}$, the spectrum obtained in the presence of plasma instabilities, and $j_{0}$, the corresponding spectra in their absence, where $j\equiv dN/dE$. This is plotted as a function of the energy in figure 7, for different spectral indices.

We find that the blazar spectral index plays a major role in the development of the cascades. Harder spectra ($\alpha\lesssim 1.5$) lead to quenching factors of $\sim 100$ at $E\sim 1\;\text{GeV}$ for $E_{\text{max}}\lesssim 10\;\text{TeV}$, whereas for $\alpha\gtrsim 2$ the quenching factors are $\lesssim 2$. The maximal energy also significantly affects the spectral shape. The higher the $E_{\text{max}}$, the larger the quenching factors for a fixed spectral index.

These results are expected. The contribution of events with higher energies is larger in the case of harder spectra (lower $\alpha$) than for softer spectra. Similarly, as $E_{\text{max}}$ increases, so does the quenching factor.

5 Application to selected blazars

In this section we discuss the consequences of plasma instabilities for VHEGR spectra of a few blazars. We select a few extreme BL Lacertae objects known to produce multi-TeV gamma rays. This list is not meant to be complete, being merely a sample of objects commonly used in studies aiming to constrain IGMFs with gamma rays. We list these object in table 3.

Note that the parameters $\alpha$ and $E_{\text{max}}$, presented in table 3, are used only to illustrate the effects of the plasma instabilities on the gamma-ray flux and for a qualitative discussion. It is beyond the scope of this work to perform a fit of the model to the data.

In figure 8 we compute the expected gamma-ray fluxes considering the effect of plasma instabilities. As discussed in section 4, the predictions for model B are compatible with a negligible action of the instabilities and are ignored in this section. We consider the range of parameters discussed in the section 4 and represent these uncertainties as bands.

From figure 8 we can draw our first conclusion, namely that the cooling due to the instabilities can cause a substantial hardening of the spectrum at $E\lesssim 100\;\text{GeV}$. The spectral signatures of this effect are similar to those of IGMFs.

We note that the spectral suppression due to the plasma instability depends on the intrinsic spectrum of the source. For instance, this effect is small for objects like 1ES 1218+304 and 1ES 1312-423, whose (fitted) spectral indices are $\alpha\approx 1.9$. On the other hand, for 1ES 0229+200, whose spectral index is $\alpha\approx 1.2$ (Taylor et al., 2011), the predicted flux varies by several orders of magnitude at $E\lesssim 100\;\text{GeV}$. For this reason, the spectra of most of the objects shown in figure 8 could change considerably, improving the agreement with the observations.

For 1ES 0229+200, Costamante et al. (2018) found that $\alpha=1.5$ and $E_{\text{max}}=12\;\text{TeV}$ adjusts the data satisfactorily. On the one hand, the increase in the spectral index would decrease the width of the bands and mitigate the effects of plasma instabilities; on the other hand, the increase in the maximal energy from 5 TeV to 12 TeV may increase this effect as $\alpha$ and $E_{\text{max}}$ are degenerate parameters.

We have considered blazars whose variabilities are small over the period of time of observation, with exception of PKS 2155-304. This object presents a state of enhanced emission during a flaring episode whose contribution is excluded from the analysis, as done by Abramowski et al. (2014).

The uncertainty due to the choice of the EBL model (Gilmore et al. (2012) in this work) does affect our results quantitatively. Qualitatively, the same conclusions as the ones presented here hold. In this case, the breadth of the band could be enlarged or reduced.

The luminosity of the beam emitted by the object can affect the observed gamma-ray flux, as shown in figure 6. This is particularly important for model C. Nevertheless, all objects listed here have very similar luminosities, of the order of $\mathcal{L}\sim 10^{43.5}-10^{45}\;\text{erg/s}$ (see e.g. Broderick et al. (2012)).

One could think that for the combinations of $\alpha$ and $E_{\text{max}}$ adopted, it would be possible to constrain some plasma instability models using the plots from figure 8. However, our calculations are meant to maximise the effects of the instabilities. A change in their growth time, for instance, could render the plasma cooling much more inefficient.

6 Discussion

The relevance of plasma instabilities for the development of electromagnetic cascades is a rather controversial issue. Many authors have constrained IGMFs neglecting the possible existence of the instabilities.

In their seminal work Broderick et al. (2012) claimed that the effects of plasma instabilities place stringent bounds on the strength of IGMFs, $B\lesssim 10^{-12}\;\text{G}$. This limit was derived by comparing the growth time of the oblique instability with the Larmor radius described by the electrons. This argument holds approximately, and other factors such as the magnetic power spectrum may play an important role in determining the propagation regime of the electrons. For instance, if the coherence length of the IGMF were much smaller than the cooling length of the charged component of the cascade, then plasma instabilities could not arise as the initial beam would be disrupted to the point where it ceases to be collimated, since the particles do not move in a specific direction, but rather propagate via (magnetic) diffusion.

An important issue raised by Durrer & Neronov (2013) is that Broderick et al. (2012), in their original work, have limited their analysis to the the linear approximation, neglecting the backreaction of the beam perturbations on the growth rate. As a consequence, oblique modes of Langmuir waves would be suppressed due to non-linear Landau damping, stabilising the beam, and thus minimising the role plasma instabilities would play on the development of the cascade. This argument explains why plasma effects on the gamma-ray spectrum are not visible for model B; they are in line with the findings of Miniati & Elyiv (2013). The same arguments are made by Schlickeiser et al. (2012, 2013), though the range of parameters for which it holds is different, explaining why in this scenario (model C) the cascade is quenched by the instabilities for most combination of parameters.

Yan et al. (2019) claimed that even if plasma instabilities are considered, in particular a fast-growing oblique instability, IGMFs could still be reliably constrained using electromagnetic cascades. In this case, the IGMF strength would change by a factor $\sim 10$. At first sight, these results do not seem to agree with ours. Nevertheless, one should bear in mind that our results are limiting cases aiming at maximising the quenching of the cascade by the plasma instabilities (equation 14). Therefore, if we assume that only a small fraction of the total energy is deposited in the IGM, i.e., if the instability is not as efficient as we considered, than our results do not conflict with those by Yan et al. (2019). More studies using three-dimensional simulations are required to draw further conclusions.

The modelling of the gamma-ray flux considering plasma instabilities is non-trivial in the presence of IGMFs. For instance, the Weibel instability (subdominant and not explicitly treated here) is known to generate and amplify magnetic fields, as observed in laboratory experiments (Huntington et al., 2015). This may disrupt the beam, increasing its transverse momentum and potentially suppressing the instability. In fact, the momentum distribution of the pairs seems to be crucial to estimate the evolution of the instabilities.

Care should be taken when stating whether or not plasma instabilities quench the development of electromagnetic cascades in the IGM. Multiple kinds of instabilities can arise, and under certain conditions their cooling rates may be subdominant with respect to inverse Compton scattering. This, for instance, is the case for model C if the luminosity of the blazar beam is $\lesssim 10^{44}\;\text{erg/s}$, or if the electron-positron pairs are created in a relatively cold region of the IGM ($T_{\text{IGM}}\sim 10^{3}\;\text{K}$), or if the density of the IGM is exceedingly high. Similarly, some dependence on the temperature can also be seen for model E, although in none of the cases studied the effects of the instabilities are negligible.

Caveats exist in our analysis. First, we consider that the cooling is completely effective, which may not be true if the instability grows slower than predicted. In this sense, our results should be taken as a limiting case wherein energy dissipation via plasma instabilities is fully efficient. Longer growth times would decrease the energy-loss rates assumed, modifying the cooling rates shown in figure 2 such that inverse Compton scattering could become the dominant cooling process. Our model is one-dimensional. The effects of IGMFs were purposely neglected in the this study in order to maximise the effects of the plasma instabilities. In the presence of IGMFs, the distribution of momentum of the electrons would suffer some change. This could, in principle, reduce the role played by the instabilities.

In our analysis only the EBL model by Gilmore et al. (2012) was considered. Although qualitatively the discussion would remain unaltered if we adopted another EBL model, quantitavely it would change. In particular, the contribution of events with $E\gtrsim E_{\text{max}}$ would be more pronounced for EBL models wherein the density of target infrared photons is lower. Ultimately, it all depends on an interplay between the injected spectrum ($E_{\text{max}}$ and $\alpha$), and the distance of the object, which follows immediately from figure 1. As for the GeV part of the observed spectrum, if the flux is completely dominated by cascade photons, we expect it to retain the slope corresponding to the spectral index of the source, regardless of the choice of EBL model.

Our results suggest that the quenching factors are strongly dependent on the spectral index and maximal energy of the blazar. These parameters are determined by the underlying acceleration mechanism in the source, as well as by its opacity. Extreme blazars are thought to have spectra $\alpha\gtrsim 1.5$ (Aharonian et al., 2006a; Cerruti et al., 2015; Bonnoli et al., 2015). Our results from figure 8 suggest that the effects of plasma instabilities are relatively small for our choices of $\alpha$ and $E_{\text{max}}$, with the exception of the object 1ES 0229+200, whose $\alpha$ is the smallest among the sources studied ($\alpha=1.2$). Therefore, since $1.5\lesssim\alpha\lesssim 2.0$ for most extreme blazars, the suppression factor may be lower than $\sim 10$; for $\alpha\gtrsim 1.8$ the quenching factor is relatively small and we conclude that, unless $E_{\text{max}}$ is high, then the spectral changes due to the instabilities are small.

Taking the parameters from table 3, we see that plasma instabilities do not explain the spectra of some blazars such as 1ES 0347-121, 1ES 0414+009, 1ES 1101-232, and 1ES 1218+304, as can be seen in figure 8. This leaves room for IC cooling of the pairs to dominate, and for IGMF signatures to be seen in blazar spectra. A detailed fit of the fluxes of these objects measured by the Fermi telescope and IACTs for different combinations of $n_{\text{IGM}}$ and $T_{\text{IGM}}$ could change the values of the best-fit of $\alpha$ and $E_{\text{max}}$; it is not within the scope of this work to perform such an analysis. Nevertheless, it is unlikely that the fitted maximal energy would change considerably if the model is altered to include plasma instabilities. As a consequence, as one can see in figure 7, the quenching factors would not be enough to explain the fluxes from the aforementioned objects for this set of parameters.

We have treated $n_{\text{IGM}}$ and $T_{\text{IGM}}$ as independent quantities. In reality, their evolutions are correlated as the temperature affects the reionisation of the gas and consequently the density of the IGM plasma (McDonald et al., 2001). This relation may be useful to decrease the number of degrees of freedom when attempting to probe the IGM with VHEGRs.

The energy dissipated by plasma instabilities may be absorbed by either the IGM, thus causing beam cooling, or by the beam itself, leading to its disruption (Broderick et al., 2012). This difference is not important when it comes to interpreting the spectra of TeV-emitting blazars, but it is of fundamental importance for discussing other observational consequences of this phenomenon such as its effect on structure formation (Pfrommer et al., 2012) and on IGM heating (Chang et al., 2012). In particular, the thermal history of the Universe depends on the fraction of the total energy that is used to heat the IGM; this efficiency factor is $\sim 10\%$ for model D, and $\sim 50$–$100\%$ assuming model A.

Our analysis focuses solely on blazars. In the context of the AGN paradigm, blazars have their jets pointing directly at Earth. Small deviations from the line of sight are not unexpected. Moreover, the jets have characteristic opening angles of $\sim 5^{\circ}$. This misalignment of the jet emission, too, may affect the growth of instabilities, since the density of the beam decreases as the distance to the axis of the jet increases. While this effect is likely small, a careful modelling of the cascade and three-dimensional simulations would be required to confirm this picture for the specific case of misaligned blazars.

For typical blazars, only after $\sim 300$ years (Broderick et al., 2012) can the instabilities grow enough to be able to cool the electron-positron pairs. Thus, flaring objects may not exhibit the same kind of behaviour as steady sources. As a consequence, they may still be reliable probes of IGMFs.

The spectral changes stemming from plasma instabilities or IGMFs are not unique. Essey et al. (2010, 2011a, 2011b) have suggested that ultra-high-energy cosmic rays (UHECRs) can induce electromagnetic cascades in the IGM, with similar observational signatures.

7 Conclusions and Outlook

We developed a numerical tool to estimate the cooling rate of electrons due to plasma instabilities caused by interactions between a beam of very energetic particles emitted by blazars. We compared different kinds of instabilities and various models found in the literature. We applied this tool to simulate the spectrum of several blazars.

Our main conclusion is that the blazar spectrum determines the impact of plasma instabilities on the development of electromagnetic cascades in the IGM. In other words, the relevance of the instabilities is related to an interplay between the blazar spectral index and maximal energy. Therefore, a detailed understanding of the VHE emission by blazars is required to determine whether plasma instability dominates over inverse Compton cooling.

We stress here, once more, that our results have to be understood as the limiting case under the assumption of maximal growth of the dominant instability for each model, as our goal was to derive robust lower limits for the observed photon flux. In the future, we intend to simulate more realistic cases and account for other processes which may decrease the influence of the instabilities on the development of electromagnetic cascades.

We have shown that for typical IGM parameters and beam luminosities some types of plasma instabilities may cool electron-positron pairs faster than inverse Compton scattering. As a consequence, these instabilities may lead to a hardening of the spectrum of blazars at energies $\lesssim 100\;\text{GeV}$. This effect resembles the suppression of the flux caused by IGMFs. Therefore, IGMF constraints may be compromised due to a possible dominance of plasma instability cooling over inverse Compton scattering during the development of electromagnetic cascades in the IGM. This result, however, depends on the hardness of the intrinsic spectrum of the blazar.

The existence of pervasive IGMFs cannot be excluded. Therefore, future studies combining both the effects of IGMFs and plasma instabilities are required to predict gamma-ray fluxes from blazars more accurately, especially for those sources for which the plasma instabilities (even overestimated) do not seem to be responsible for the hardening of the $\lesssim 100$ GeV part of the spectrum (i.e., 1ES 0347- 121, 1ES 0414+009, 1ES 1101-232, and 1ES 1218+304, as seen in figure 8). Moreover, efforts towards simulating the growth of these instabilities using magnetohydrodynamical or particle-in-cell simulations would be needed to unambiguously determine whether or not they can quench electromagnetic cascades.

An uncontroversial window of opportunity for constraining IGMFs with electromagnetic cascades remains. The search for echoes associated with transient events, i.e., magnetically-induced delays in the arrival secondary gamma rays from primary VHEGR emission, are still excellent probes of IGMFs because the effects of plasma instabilities may be small. This is the case for some particular energetic events in AGNs. Gamma-ray bursts seem very promising candidates for this purpose.

Acknowledgments

RAB is supported by grant #2017/12828-4, São Paulo Research Foundation (FAPESP). The work of AS was supported by the Russian Science Foundation under grant no. 19-71-10018, carried out at the Immanuel Kant Baltic Federal University. EMGDP is partially supported by FAPESP grant #2013/10559-5 and CNPq grant 306598/2009-4. The simulations were run at the cluster GAPAE of the group of Plasmas and High-Energy Astrophysics of IAG-USP (FAPESP grant #2013/10559-5).