Multi-TeV Flaring from High-energy Blazars: An Evidence of the Photohadronic Process

Blazars are a subclass of active galactic nuclei (AGNs) that includes flat spectrum radio quasars (FSRQs) and BL Lacertae (BL Lac) objects (Romero et al. 2017). These objects are characterized by nonthermal spectra at all wavelengths, from radio to very high energy (VHE; >100 GeV) γ-rays and show flux variability on timescales ranging from months to a few minutes (Abdo et al. 2010a). The flux variability is produced in a highly relativistic jet pointing toward the observer. The spectral energy distributions (SEDs) of blazars are characterized by two nonthermal peaks (Abdo et al. 2010b). The first peak (low energy) is located between infrared to X-ray energies, produced from the synchrotron emission from the relativistic electrons in the jet. The general consensus is that the second peak (high energy) corresponds to the synchrotron self Compton (SSC) scattering of the high-energy electrons with their self-produced synchrotron photons. Depending on the location of the first peak, blazars are often subdivided into low-energy peaked blazars, intermediate-energy peaked blazars, and high-energy peaked blazars (HBLs; Abdo et al. 2010b). The leptonic model is very successful in explaining the multiwavelength emission from blazars (Tavecchio et al. 2011; Boettcher et al. 2013). The nearest HBL, Markarian 421 (Mrk 421), was the first to be detected in TeV energy by Whipple telescopes (Punch et al. 1992). In recent years, the highly sensitive Imaging Atmospheric Cerenkov Telescopes (IACTs) such as VERITAS (Holder et al. 2009), HESS (Hinton 2004), and MAGIC (Cortina 2005), have had great success in discovering many new extragalactic TeV sources, most of which are blazars. Consequently, blazars are an important class of objects to observe and study VHE gamma-ray astronomy.

Flaring in VHE γ-rays seems to be the major activity of many HBLs. This flaring is unpredictable and switches between quiescent and active states involving different timescales (Sentürk et al. 2013). It has been observed that while in some blazars a strong temporal correlation between X-ray and multi-TeV γ-ray exists, in others (except for VHE γ-rays) no low-energy counterpart is observed; such an anti-correlation is difficult to explain using a leptonic model (Krawczynski et al. 2004; Blazejowski et al. 2005). Different models have been developed to explain these flaring events (Giannios et al. 2010; Cerruti et al. 2015). Many simultaneous multiwavelength observations have been made to construct the SED of the flaring period to constrain different theoretical models (Sentürk et al. 2013; Ahnen et al. 2017).

The propagating VHE γ-rays undergo energy-dependent attenuation by the intervening extragalactic background light (EBL) through pair production (Ackermann et al. 2012), and the EBL significantly changes the shape of the VHE spectrum. Therefore, a proper understanding of the EBL SED is important for the calculation of the intrinsic spectrum. Well-known EBL models are used by the IACT collaborations in order to analyze the observed VHE γ-rays from objects of different redshifts (Franceschini et al. 2008; Domínguez et al. 2011).

By assuming that the multi-TeV emission in the HBLs are due to the photohadronic interaction in the jet (Sahu et al. 2017b; Sahu 2019), a simple relation between the observed VHE spectrum and the intrinsic spectrum is derived. We assume that during the VHE emission period the Fermi-accelerated protons with a power-law spectrum (Dermer & Schlickeiser 1993), dN/dE ∝ E^−α (the power index α ≥ 2), interact with the background seed photons in the jet to produce the Δ-resonance ($p\gamma \to {{\rm{\Delta }}}^{+}$), which subsequently decays to γ-rays via intermediate π⁰ and to neutrinos through π⁺. In a canonical jet scenario, the Δ production efficiency is very low due to the low photon density. Therefore, to explain the multi-TeV emission through this process, super-Eddington power in the proton is needed (Cao & Wang 2014). To circumvent this problem a double-jet structure scenario is proposed (Sahu 2019): a small compact cone is enclosed by a bigger one along the same axis, and the photohadronic interaction occurs in the inner jet region. The photon density ${n}_{\gamma ,f}^{{\prime} }$ in the inner compact region is much higher than that of the outer region ${n}_{\gamma }^{{\prime} }$ (where prime corresponds to jet comoving frame) and, due to the adiabatic expansion of the inner jet, its photon density decreases by crossing into the outer region. As the photon density is unknown in the inner jet region, we assume a scaling behavior of the photon densities in the inner and the outer jet regions, which essentially means that the spectra of the outer and the inner jets have the same slope. Using this scaling behavior, we can express the photon density in the inner region in terms of the photon density of the outer region, which is known from its observed SED.

The kinematical condition to produce the Δ-resonance is given by Sahu (2019)

$$\begin{eqnarray}&&{E}_{\gamma }{\epsilon }_{\gamma }=0.032\,{\rm{\Gamma }}{ \mathcal D }{\left(1+z\right)}^{-2}\,{\mathrm{GeV}}^{2},\end{eqnarray} \tag{ 1 }$$

where E_γ, _γ, Γ, ${ \mathcal D }$, and z are the observed VHE γ-ray, seed photon energy in the observer's frame, bulk Lorentz factor, Doppler factor, and redshift, respectively. For an HBL, ${\rm{\Gamma }}\simeq { \mathcal D }$ is satisfied. The observed VHE γ-ray flux depends on the Fermi-accelerated proton flux F_p and the background seed photon density ${F}_{\gamma ,\mathrm{obs}}(={E}_{\gamma }^{2}\,{{dN}}_{\gamma }/{{dE}}_{\gamma })\propto {F}_{p}\,{n}_{\gamma ,f}^{{\prime} }$. Also ${F}_{p}\propto {E}_{\gamma }^{-\alpha +2}$, and using the scaling behavior we can express ${n}_{\gamma ,f}^{{\prime} }\propto {\rm{\Phi }}({\epsilon }_{\gamma }){\epsilon }_{\gamma }^{-1}$, where Φ is the observed/fitted flux corresponding to seed photon energy _γ. Previously, the photohadronic model has been successfully used to explain many flaring HBLs and has found that, for all the cases studied so far, Φ lies in the tail region of the SSC SED (Sahu et al. 2017a). But this region of the SED is not observed/measured due to technical difficulties. Mostly leptonic models are used to calculate the flux in this region, and different leptonic models predict different fluxes. In the same HBL, the flux in this region varies during different flaring states and also different epochs. However, irrespective of the model used, the predicted flux in the tail region of the SSC SED is a power law given by ${\rm{\Phi }}\propto {\epsilon }_{\gamma }^{\beta }$ and, using the above kinematical condition, we can re-express it as ${\rm{\Phi }}\propto {E}_{\gamma }^{-\beta }$. Putting everything together and taking into account the EBL correction, the observed VHE γ-ray spectrum can be expressed as the product of the intrinsic flux F_γ,int and the attenuation factor due to e⁺ e⁻ pair production as

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\gamma ,\mathrm{obs}}({E}_{\gamma }) & = & {F}_{\gamma ,{int}}({E}_{\gamma }){e}^{-{\tau }_{\gamma \gamma }({E}_{\gamma },z)}\\ & = & \ {F}_{0}{\left(\displaystyle \frac{{E}_{\gamma }}{{TeV}}\right)}^{-\delta +3}{e}^{-{\tau }_{\gamma \gamma }({E}_{\gamma },z)},\end{array}\end{eqnarray} \tag{ 2 }$$

where, F₀ is the normalization constant, and δ = α + β. The optical depth τ_γγ is a function of E_γ and z. F₀ and δ are the only parameters to be adjusted to fit the observed spectrum. However, strictly speaking the normalization constant is not a free parameter that can be fixed from the observed data. It is not necessary to know a priori the value of β, but it can be constrained by fitting the observed data with the parameter δ. Moreover, the spectral index of the intrinsic differential spectrum can be defined as δ_int = −δ + 1.

The stability of the inner jet on large scales can be estimated from the ratio σ of the magnetic stress (Poynting flux) and the kinetic stress and for BL Lac objects σ ≲ 1. By considering the generic values of the parameters magnetic field B ∼ 1 G, proton density n_p ∼ 10⁻¹–10⁻² cm⁻³, and bulk Lorentz factor Γ ∼ 10, we obtain σ ∼ 0.4, which corresponds to a stable inner jet (Cavaliere et al. 2017). The photon density within the inner jet region can be constrained by comparing the jet expansion timescale ${t}_{d}^{{\prime} }$ with the pγ interaction timescale ${t}_{p\gamma }^{{\prime} }$ and assuming that the high-energy proton luminosity is smaller than the Eddington luminosity (Sahu et al. 2016).

Using Equation (2), we fitted the observed VHE spectra of 42 emission epochs of 23 HBLs of different redshifts very well with the free parameter δ in the range 2.5 ≤ δ ≤ 3.0. Depending on the value of δ, we roughly classify these flaring states in three different categories: (i) low state, when δ = 3.0, (ii) high state, when 2.6 < δ < 3.0, and (iii) very high state, when 2.5 ≤ δ ≤ 2.6. We know a priori that α ≥ 2, so during the simultaneous observation period in multiwavelength, we must have 0.0 ≤ β ≤ 1.0. The three different emission states are discussed through four examples with HBLs of different redshifts, and the EBL model of Franceschini et al. (2008) is used for our analysis.

3.1. 1ES 0229+200

The 1ES 0229+200 is an HBL at a redshift of z = 0.1396 that was discovered in the Einstein IPC Slew Survey in 1992 (Schachter et al. 1993). It was observed by VERITAS telescopes during a long-term observation over three seasons between 2009 October and 2013 January, for a total of 54.3 hr (Aliu et al. 2014a), and an excess of 489 γ-ray events were detected in the energy range 0.29 TeV ≤ E_γ ≤ 7.6 TeV. Using the proton-synchrotron model and the lepto-hadronic model, dominated by emission from the secondary particles from pγ interactions, the observed multiwavelength SEDs of several HBLs are fitted by Cerruti et al. (Cerruti et al. 2015). However, these numerical models use about 19 parameters to fit the entire SED. We show their fit to 1ES 0229+200 in Figure 1. Alternatively, using the photohadronic model an excellent fit is obtained for δ = 2.6 and F₀ = 3.5 × 10⁻¹² erg cm⁻² s⁻¹. According to the previously discussed classification scheme, this corresponds to very high emission state and the intrinsic flux ${F}_{\mathrm{int}}\propto {E}_{\gamma }^{0.4}$. Similarly, the extracted differential spectrum ${({{dN}}_{\gamma }/{{dE}}_{\gamma })}_{\mathrm{int}}\propto {E}_{\gamma }^{-1.6}$, which is not hard. This HBL has the central black hole of mass M_BH ∼ 1.4 × 10⁹ M_⊙ and outer blob size ${R}_{b}^{{\prime} }\sim {10}^{16}-{10}^{17}$ cm (Zacharopoulou et al. 2011). Assuming that the high-energy proton flux corresponding to E_γ = 7.6 TeV is smaller than the Eddington flux and comparing ${t}_{d}^{{\prime} }$ (inner blob size ${R}_{f}^{{\prime} }\sim 4\times {10}^{15}$ cm) with ${t}_{p\gamma }^{{\prime} }$, we obtain the photon density in the range $4\times {10}^{8}\,{\mathrm{cm}}^{-3}\lt {n}_{\gamma ,f}^{{\prime} }\,\lt 2.5\,\times {10}^{11}\,{\mathrm{cm}}^{-3}$.

Zoom In Zoom Out Reset image size

From this HBL between 2005 and 2006 the HESS telescopes also observed VHE γ-rays (Aharonian et al. 2007a), whose time-averaged spectrum is in the energy range 0.5 TeV ≤ E_γ ≤ 11.5 TeV and is very similar to the one discussed above. This spectrum is fitted with the hadronic model of Essey et al. (2010, 2011b). Using the photohadronic model a very good fit is obtained for δ = 2.5 (see Figure 5.3 of Figure Set 5 for details). By reducing 10% to the hadronic model of Essey et al. the spectrum of VERITAS can be fitted well, which is shown in Figure 1 for comparison.

3.2. 1ES 0347-121

The 1ES 0347-121 is a HBL at a redshift of z = 0.188. The HESS telescopes observed this blazar between 2006 August and 2006 December for a total of 25.4 hr (Aharonian et al. 2007a) when an excess of 327 VHE gamma-ray events were detected in the energy range 0.25 TeV ≤ E_γ ≤ 3 TeV and no flux variability was detected in the data set.

In a hadronic model scenario, ultra-high energy protons escaping from the jet produce secondary VHE gamma-rays by interacting with the cosmic microwave background (CMB) and/or EBL (Essey et al. 2011b). Assuming this scenario the spectra of 1ES 0347-121, 1ES 0229+200, and 1ES 1101-232 are explained well (Essey et al. 2011b). However, this scenario requires protons in the energy range 10⁸–10¹⁰ GeV, which are not easily produced in the jet environment, as well as a weak extragalactic magnetic field in the range 10⁻¹⁷ G < B < 10⁻¹⁴ G to produce the observed gamma-ray spectrum along the line of sight (Essey et al. 2011a). In an alternative scenario, Cerruti et al. (2015) have applied the proton-synchrotron and lepto-hadronic models to fit the spectrum of 1ES 0347-121. Using the photohadronic model we found an excellent fit to the spectrum with δ = 2.7 and F₀ = 6.0 × 10⁻¹² erg cm⁻² s⁻¹, which is a high-state emission. As a γ-ray carries ∼10% of the proton energy (Sahu 2019), E_γ = 3 TeV corresponds to 30 TeV cosmic ray proton energy, which can easily be produced and accelerated in the blazar jet. In Figure 2 we compare our result with Essey et al. (2011b) and Cerruti et al. (2015) and found that below 1 TeV all have similar behaviors. However, above 1 TeV our result differs substantially from the others, particularly that of Essey et al., which uses the EBL model of Stecker et al. (high EBL; Stecker et al. 2006). Comparison of the EBL models of Franceschini et al. (2008) and Stecker et al. (2006) shows a significant difference in the attenuation factor above 1 TeV.

Zoom In Zoom Out Reset image size

3.3. 1ES 0806+524

The 1ES 0806+524 is at a redshift of z = 0.138 and in 2008, the VERITAS telescopes discovered this in VHE γ-rays (Acciari et al. 2009). A multiwavelength observation was performed by MAGIC telescopes from 2011 January to 2011 March for 13 nights for about 24 hr (Aleksić et al. 2015b) and observed a flaring event on February 24. Within 3 hr of observation excess events above 250 GeV were recorded in the energy range 0.17 TeV ≤ E_γ ≤ 0.93 TeV, when the flux increased by a factor of about 3 from the mean flux level and no intra-night variability was observed. The flaring data and the remaining MAGIC observations are analyzed separately using the photohadronic scenario, which are shown in Figure 3. Using one-zone SSC model, the broadband SEDs during the flaring (high) and the remaining period (low) are explained using about 14 free parameters. The electron Lorentz factor for the high state is double the one for the low state, and the remaining parameters are the same (Aleksić et al. 2015b). With the photohadronic scenario, the flaring state can be fitted very well with δ = 2.9 and F₀ = 1.2 × 10⁻¹¹ erg cm⁻² s⁻¹ corresponding to a high emission state and the average of the remaining flux can be fitted with δ = 3.0 and F₀ = 4.0 × 10⁻¹² erg cm⁻² s⁻¹, which is a low state. The intrinsic fluxes in high and low states, respectively, are proportional to ${E}_{\gamma }^{0.1}$ and ${E}_{\gamma }^{0}$. Comparison of both the models is shown in the Figure 3. The SSC model does not fit well to the low-state spectrum. Although both models explain well the flaring data, a significant difference is observed in the predictions above 1 TeV.

Zoom In Zoom Out Reset image size

3.4. HESS J1943+213

HESS J1943+213 is a VHE gamma-ray point source discovered by HESS (Abramowski et al. 2011) that is identified as an extreme HBL (EHBL). In VHE, it was observed by VERITAS telescopes from 2014 May 27 to 2014 July 2 and from 2015 April 20 to 2015 November 9, a total exposure time of 37.2 hr, and no flux variability was observed (Archer et al. 2018). The time-averaged spectrum of both the observations is presented in Figure 4. Currently, the redshift of HESS J1943+213 is not known and indirect limits (0.03 < z < 0.45) were set by Peter et al. (2014). Improved gamma-ray spectra of Fermi-LAT and VERITAS were used to derive a conservative upper limit of z < 0.23 (Archer et al. 2018). Using the photohadronic model and different redshifts, we derived more stringent lower and upper limits on the redshift (0.14 ≤ z ≤ 0.19), which are shown in Figure 4. However, the best fit is obtained for z = 0.16 and δ = 2.9, corresponding to a high-state emission from the source. Two additional examples are discussed in the supplementary materials.

Zoom In Zoom Out Reset image size

The HBLs are known to undergo episodes of VHE flaring in gamma-rays involving different timescales, and the flaring mechanism is not well understood. Also, the VHE gamma-rays are attenuated by EBL background. Here we have derived a simple relation between the observed VHE flux and the intrinsic flux from the flaring HBLs. This was accomplished by assuming that during flaring, Fermi-accelerated high-energy protons interact with the seed photons in the inner compact region of the jet to produce Δ-resonance, which subsequently decays to gamma-rays and neutrinos from intermediate π⁰ and π⁺, respectively. These gamma-rays can be observed. To account for the EBL effect we consider the well-known EBL model of Franchesccini et al. and analyzed 23 HBLs of different redshifts with a total of 42 different emission epochs. For detailed analysis we only used five emission epochs of four HBLs, and the rest of the flaring states are summarized in Table 1. Some of these are briefly discussed in the supplementary materials to strengthen further the validity of the photohadronic origin of multi-TeV flaring events.

From the analysis we observed that the free parameter δ is constrained to be in the range 2.5 ≤ δ ≤ 3.0. The intrinsic flux for the low state is a constant, but for high and very high state it is a power law proportional to ${E}_{\gamma }^{\eta }$, where 0 < η ≤ 0.5. We could not find any flaring state that has δ < 2.5. Some flaring spectra can be fitted well with δ > 3. However, it is important to note that for these cases −δ + 3 is positive (a very soft spectrum) and in the low-energy limit the spectrum shoots up very high, which is certainly not observed. Consequently, the soft power-law fits are ignored (Dwek & Krennrich 2005; Sahu et al. 2018) and we always adhere to δ ≤ 3.0. From the analysis we observed that about 48% are low-state emissions, 38% are high-state emissions, and 14% are very-high-state emissions. This implies that low and high emission states constitute the major part of the flaring in HBLs.

Although the photohadronic scenario works well for E_γ ≳ 100 GeV, there are contributions from the leptonic processes to the observed spectrum in this energy regime, so in the low-energy regime our model may not fit the data very well. In some cases, we have observed that the averaging of long-term VHE observations are difficult to explain by photohadronic model for the following reasons: gamma-rays from the leptonic processes contribute to the spectrum in the low-energy regime, and the averaging of many unobserved short flares with the low emission periods contaminates the data.

Several models explain well the observed broadband SEDs but require many assumptions and many free parameters, some of which are difficult to realize in the jet environment (Essey et al. 2011b; Boettcher et al. 2013; Aleksić et al. 2015b; Cerruti et al. 2015). On the other hand, the photohadronic scenario is based on very simple assumptions that are very likely to be realized in the jet during the VHE emission period. Another important aspect of our model is that the assumption of the power-law behavior of the background seed photon is sufficient to fit the observed spectrum, and it is not necessary to have simultaneous multiwavelength observations. Moreover, the exact simultaneous multiwavelength observation during a flaring event in a HBL is usually limited to a few. In our case, the IACTs observations are sufficient. From the fitting to the observed spectrum and using α ≥ 2, the seed photon spectral index β can be constrained. For example, an excellent fit to the flaring of PG 1553+113 is obtained for δ = 2.5, which shrinks the β value in the interval 0–0.5 (see Table 1). Nevertheless, the fact that we can explain very well the VHE spectra of 42 epochs of 23 HBLs with a single parameter provides strong evidence that VHE gamma-rays are produced mostly through the photohadronic process with the intermediate Δ-resonance. In addition, it is important to mention that for HBLs of unknown redshifts, whose multi-TeV spectra are known, stringent bounds on the redshifts can be placed using the photohadronic model.

The work of S.S. is partially supported by DGAPA-UNAM (Mexico) project No. IN103019. S.N. is partially supported by "JSPS Grants-in-Aid for Scientific Research $\langle$ KAKENHI $\rangle$ (A) 19H00693," "Pioneering Program of RIKEN for Evolution of Matter in the universe (r-EMU)," and "Interdisciplinary Theoretical and Mathematical Sciences Program of RIKEN."

Due to space constraints in the main article, we only analyzed five flaring states of four HBLs in the context of the photohadronic model and compared it with other available models. However, to further support the validity of our model and its predictions, we provide in Figure Set 5 eleven additional flaring states of HBLs of different redshifts. Particularly, our best fits to the flaring events of 1ES 0229+200 and 1ES 1101+232 are compared with other existing leptonic and hadronic models, where we observed that our results are as good as or better than these models. The redshifts of the HBLs PKS 1440-389 and RGB J2243+203 are unknown and, using different observations, limits were set to the redshifts. We have shown that the predicted photohadronic model limits are more stringent than the existing ones. The references to all the additional HBLs given in Table 1 are shown in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (11 panels)

Tarball of figure set images