The origin of multi-TeV flares from the nearest blazar Markarian 421

Markarian 421 is a high-peaked BL Lac object and it has undergone many strong outbursts since its discovery as a TeV source in 1992. Markarian 421 has been studied intensively and was observed by various Cherenkov telescope arrays ever since. The outbursts of April 2004 observed by the Whipple telescope and of February 2010 by the HESS telescopes are explained well in this work by using the photohadronic model. To account for the attenuation of these high-energy gamma-rays by the extragalactic background light (EBL), we use template EBL models. The intrinsic spectrum of each epoch is different even though the high-energy protons have almost the same spectral index. We observe that this difference in intrinsic spectra is due to the change in the spectral index of the low-energy tail of the synchrotron self Compton (SSC) photons during different epochs of flaring. Our results show that the contemporaneous multiwavelength observations, particularly in the low-energy tail region of the SSC emission of the source, are important in explaining the flaring phenomenon.


Introduction
The extragalactic very high-energy (VHE, E γ > 100 GeV) gamma-rays undergo energy dependent attenuation en route to Earth by the intervening extragalactic background light (EBL) through electron-positron pair production [1]. This interaction process not only attenuates the absolute flux but also significantly changes the VHE emission spectrum. The a e-mail: sarira@nucleares.unam.mx b e-mail: albertoros4@ciencias.unam.mx c e-mail: shigehiro.nagataki@riken.jp diffuse EBL contains a record of the star formation history of the universe. A proper understanding of the EBL spectral energy distribution (SED) is very important for the correct interpretation of the deabsorbed VHE spectrum from the source. The direct measurement of the EBL is very difficult with high uncertainties mainly due to the contribution of zodiacal light [2,3], and galaxy counts result in a lower limit since the number of unresolved sources (faint galaxies) is unknown [4]. Keeping in mind the observational constraints, several approaches with different degrees of complexity have been developed to calculate the EBL density as a function of energy for different redshifts [5][6][7][8][9][10][11][12]. Mainly three types of EBL models exist: backward and forward evolution models and semi-analytical galaxy formation models with a combination of information as regards galaxy evolution and observed properties of galaxy spectra. In the backward evolution models [8], one starts from the observed properties of galaxies in the local universe and evolves them from cosmological initial conditions or extrapolating backward in time using parametric models of the evolution of galaxies. This extrapolation induces uncertainties in the properties of the EBL which increases at high redshifts. On the contrary, the forward evolution models [7,10] predict the temporal evolution of galaxies forward in time starting from the cosmological initial conditions. Finally, semi-analytical models have been developed which follow the formation of large scale structures driven by cold dark matter in the universe by using the cosmological parameters from observations. This method also accounts for the merging of the dark matter halos and the emergence of galaxies which form as baryonic matter falls into the potential wells of these halos.
Blazars are a subclass of active galactic nuclei (AGN) and they are the dominant extra galactic population in γrays [13]. These objects show rapid variability in the entire electromagnetic spectrum and have non-thermal spectra produced by the relativistic jet of plasma oriented close to the observers line of sight [14]. The jets are powered by matter accretion onto a supermassive black hole. The spectral energy distribution (SED) of these blazars has a double peak structure in the ν-ν F ν plane. The low-energy peak corresponds to the synchrotron radiation from a population of relativistic electrons in the jet and the high-energy peak is believed to be due to the synchrotron self Compton (SSC) scattering of the high-energy electrons with their self-produced synchrotron photons [15,16]. As leptons (e ± ) are responsible for the production of the SED, this is called the leptonic model and in general it is very successful in explaining the multiwavelength emission from blazars and FR I galaxies [17][18][19][20].
Blazars detected in VHE are predominantly high-energy peaked blazars (HBLs) whose synchrotron peak lies in the UV to X-ray energy range and the inverse Compton peak is in the GeV-TeV energy range [14,18]. Flaring seems to be a major activity of these blazars, which is unpredictable and switches between quiescent and active states involving different time scales and flux variabilities [21][22][23][24][25]. However, broadly speaking the flaring mechanism is unknown.
Mrk 421 is the first extragalactic source detected in the multi-TeV domain [45] and it is one of the fastest varying γ -ray sources. There are other BL Lac objects with lower redshifts than Mrk 421; however, these objects were never reported as TeV emitters and could have been misclassified [46]. Also, there are many BL Lacs with unknown redshifts.
Several large flares of Mrk 421 were observed in 2000-2001 [47][48][49] and 2003-2004 [36,50]. During April 2004, a large flare took place both in the X-rays and in the TeV energy band. The source was observed simultaneously at TeV energies with the Whipple 10 m telescope and at Xray energies with the Rossi X-ray Timing Explorer (RXTE) [36]. It was also observed simultaneously in radio and optical wavelengths. During the flaring period, the TeV flares had no coincident counterparts at longer wavelengths and it was observed that the X-ray flux reached its peak 1.5 days before the TeV flux did during this outburst. The orphan TeV flare in 1ES 1959 + 650 of 2002 [35] had also no simultaneous X-ray counterpart and the variation patterns in X-rays were similar to the one observed in Mrk 421. A strong outburst in multi-TeV energy in Mrk 421 was first detected by the VERITAS telescopes on 16th of February 2010 and follow up observations were done by the HESS telescopes during four subsequent nights [30].
In the framework of a 6-month long multi-instrument campaign, the MAGIC telescopes observed VHE flaring from Mrk 421 on 25th of April 2014 and the flux (above 300 GeV) was about 16 times brighter than usual. This triggered a joint ToO program by XMM-Newton, VERITAS, and MAGIC. These three instruments individually observed approximately 3 h each day on April 29, May 1, and May 3 of 2014 [51]. The simultaneous VERITAS-XMM-Newton observation was published recently and it was shown that the observed multiwavelength spectra are consistent with onezone synchrotron self-Compton model [51]. However, the details of the large flare observed on 25th April by MAGIC are not yet publicly available.
For Mrk 421, being the nearest HBL, the VHE gammarays can be attenuated the least and in principle presently operating Cherenkov telescopes will be able to observe higher-energy photons. At the same time, different EBL models can also be compared. So Mrk 421 can be used as an example to study the intrinsic emission mechanism in the multi-TeV energy range from the nearest blazars and test different EBL models.
From the above various VHE emission epochs, we shall analyze the following events: (1) the flare of April 2004 [36] and (2) the flare of February 2010 [52]. Both experimental data sets are public and have multiwavelength information, which is important to improve the photohadronic fitting of the flare. Also both flares were having a high flux, which is important for our study. The flare of 2004 was analyzed by Sahu et al. [53] where an EBL correction was not considered. In this work we include the EBL correction for the same 2004 data. For our analysis of these flaring events, we will be using the photohadronic model; a detailed account of this model is given in Refs. [53][54][55][56][57][58][59]. This model is successfully employed to explain the multi-TeV flaring from many high frequency blazars. A brief account of the model is given in the next section.

Photohadronic scenario
In the photohadronic scenario, the Fermi accelerated protons, having a power-law spectrum, dN /dE p ∝ E −α p , interact with the background photons to produce the Δ-resonance, which subsequently decays to gamma-rays via an intermediate neutral pion and to neutrinos through a charged pion [54]. We assume that this interaction occurs within a compact and confined volume of radius R f inside the blob of radius R b (R f < R b ). Here the notation implies the jet's comoving frame. Geometrically this represents a double jet structure, one compact and smaller cone which is enclosed by a bigger one along the same axis; the geometry of this model is discussed in Fig. 1 of Ref. [53]. The inner compact region has a photon density much higher than the outer region. Due to the adiabatic expansion of the inner jet, the photon density will decrease when it crosses into the outer 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. Mathematically we can express this as i.e. the ratio of photon densities at two different background energies γ 1 and γ 2 in the flaring (n γ,f ) and in the non-flaring (n γ ) states remains almost the same. The photon density n γ in the outer region can be calculated in the usual way using the observed/fitted SED. Afterwards, by using the above relation in Eq. (1), we can express the unknown inner photon density in terms of the outer known density. Henceforth, for our calculation, we shall use n γ and its corresponding flux rather than the one from the inner jet region. The observed VHE γ -ray flux depends on the background seed photon density and the differential powerspectrum of the Fermi accelerated protons given as F γ ∝ n γ (E 2 p dN /dE p ). It is to be noted that the photohadronic process in a standard blazar jet environment (quiescent state) is inefficient due to the low seed photon density n γ . So to explain the multi-TeV emission from the flaring in the photohadronic scenario, jet kinetic power has to be increased to the super-Eddington limit [40,41]. However, the inner compact jet scenario evades this problem due to the higher photon density [55]. So, the assumption here is that the outer jet is always there and is responsible for the quiescent state of the blazar, while the inner jet is transient and is responsible for the flare. The photohadronic model provides a simple explanation for the multi-TeV observed events from the HBLs and it depends only on the seed photons from the low-energy tail of the SSC emission. On the other hand, the one-zone leptonic model does not fit well to the observed multi-TeV data and the multi-zone leptonic model needs many more free parameters (compared to the one-zone leptonic model) to fit it.
The observed γ -ray energy E γ from the π 0 decay and the background seed photon energy γ satisfy the kinematical condition following the resonance process pγ → Δ and the incident proton has the energy E p 10E γ . Here D is the Doppler factor and the bulk Lorentz factor of the blob is given as Γ .
For blazars we have Γ D. The different epochs of VHE flaring between 2004 and 2014 of Mrk 421 had different ranges of observed γ -rays, which correspond to different ranges of the seed photon energies. We observed that, irrespective of the E γ and epoch of flaring, the ranges of γ are always in the low-energy tail region of the SSC emission. This is the region in the valley where the SSC emission begins (see the red line in Fig. 2 where 0.26 ≤ E γ ≤ 100 MeV).
The SSC flux in this range of seed photon energy is exactly a power law, given by Φ SSC ∝ β γ with β > 0. Again, from the kinematical condition to produce Δ-resonance through pγ interaction, γ can be expressed in terms of E γ and can be written as From the leptonic model fit to the observed multiwavelength data (up to second peak) during a quiescent/flaring state we can get the SED for the SSC region from which Φ 0 and β can be obtained easily. One expresses the observed flux F γ in terms of the intrinsic flux F γ,in and the EBL correction as where the intrinsic flux is given in Ref. [58] as where A γ is a dimensionless normalization constant that can be fixed by fitting the observed VHE data. As discussed above, the power index β is fixed from the tail region of the SSC SED for a given leptonic model which fits the lowenergy data well. So the Fermi accelerated proton spectral index α is the only free parameter to fit the intrinsic spectrum.
To account for the EBL correction to the observed multi-TeV gamma-rays, here we use the EBL models of Dominguez et al. [5] (EBL-D) and Inoue et al. [6] (EBL-I) to interpret our results. Figure 1 presents the attenuation factor for these two models as functions of observed gamma-ray energy E γ for z = 0.031. Figure 1 shows that the difference between the two models becomes apparent at 600 GeV but continues to increase above 1 TeV; the difference is maximum around 2 TeV, and then it decreases until the models intersect at 6 TeV. At higher energies the models diverge and they converge again above 20 TeV with the same attenuation factor. The Fermi accelerated protons in the jet will emit synchrotron radiation but it will be suppressed by a factor of m −4 p , where m p is the proton mass. Also the emission from the ultra-high energy protons needs a stronger magnetic field. The same ultra-high energy protons can leak out from the jet region and can reach the Earth as ultra-high energy cosmic rays (UHECRs). The charged pions produced from the photohadronic process will decay to neutrinos which can in Fig. 1 The attenuation factor as a function of photon energy predicted by the two EBL models for z = 0.031 principle be detected on Earth [60][61][62]; however, the flux is too low.

Two flaring episodes
As discussed in the introduction, several epochs of major multi-TeV emission/flaring are observed from Mrk 421 and we shall analyze two of these multi-TeV flares observed by different Cherenkov telescope arrays and use the photohadronic model to interpret these flaring events. As input to the photohadronic model we use the leptonic SSC SED fitted to the multiwavelength observations.

The flare of April 2004
The multi-TeV flare of April 2004 was the first flare observed in multiwavelength and it was difficult to explain by the onezone leptonic model [36]. The Whipple telescope observed the flare in the energy range 0.25 TeV(6.0 × 10 25 Hz) ≤ E γ ≤ 16.85 TeV(4.1 × 10 27 Hz). We use the one-zone leptonic model of Ref. [36] (lep-1) as input for the photohadronic model to explain the observed TeV emission. Here the bulk Lorentz factor associated with the lep-1 model is Γ = D = 14. The Fermi accelerated proton energy lies in the range 2.5 TeV ≤ E p ≤ 168 TeV and the corresponding background photon energy is in the range 23.6 MeV(5.7×10 21 Hz) ≥ γ ≥ 0.35 MeV (8.4×10 19 Hz). This range of γ is in the low-energy tail region of the SSC SED and the corresponding flux follows an exact power law  Fig. 2.
For the EBL correction to the observed data, we use EBL-D and EBL-I. The best fit to the observed multi-TeV flare data is obtained for the spectral index α = 2.7, and the normalization constant A γ for EBL-D is 2.3 and for EBL-I is 2.5, respectively, which are shown in Fig. 3. For E γ < 4 TeV both EBL models fit the data very well, and above this energy there is a slight difference due to the change in the attenuation factor. We observe that the EBL-D fit is better than the EBL-I fit. Again above 20 TeV the EBL models behave in the same way. For comparison we have also shown in the same figure (red curve) the photohadronic model fit without EBL correction but with an exponential cut-off [53], and the multi-zone leptonic model fit (magenta curve) [36]. It can be seen from Fig. 3 that the multi-zone fit is not so good as other fits for E γ ≤ 15 TeV. However, for higher energy it has the same behavior as EBL-D and EBL-I. With the exponential cut-off scenario, a good fit is obtained for the spectral index α = 2.7 and the cut-off energy E c = 6.2 TeV. Again comparing this with the EBL corrected models, below 4 TeV, all these fits are exactly the same. However, above 4 TeV we observe some discrepancy among these fits and above 10 TeV the fits of EBL-D and EBL-I fall faster than the exponential cut-off scenario. It is clear from the comparison of EBL-I with the exponential cut-off scenario that, for E γ ≤ 10 TeV, e −E γ /E c and e −τ γ γ are almost the same, and above 10 TeV the attenuation factor falls faster than the exponential cut-off for which the EBL-D and EBL-I curves fall fast. photohadronic model using two different EBL models are shown. It is also compared with the power law with exponential cut-off without EBL correction fit [53] and with the multi-zone leptonic fit [36]. The multizone leptonic model accounts for the attenuation of the very high-energy gamma-rays by the diffuse infrared background. The intrinsic fluxes for the two EBL models are also shown The intrinsic fluxes in EBL-D and EBL-I are also shown. They are almost the same and show a power-law behavior with F γ,in ∝ E −0.18 γ . Even though all these models fit quite well to the observed data below the 20 TeV energy range, the deviation is appreciable above 20 TeV between the EBL corrected plots and the exponential cut-off. So observation of VHE flux above ∼ 30 TeV will be a good test to constrain the EBL effect on the VHE gamma-rays from Mrk 421 and to see whether the energy cut-off scenario is necessary or not.
In Ref. [53], by comparing the expansion time scale, the interaction time scale of the pγ interaction and the high-energy proton luminosity are expected to be smaller than the Eddington luminosity in the inner jet region R f 3 × 10 15 cm (the size of the outer jet is R b 0.7 × 10 16 cm), and the range of the optical depth for the Δ-resonance production is estimated as 0.02 < τ pγ < 0.13. This corresponds to a photon density in the inner jet region of 1.3 × 10 10 < n γ,f < 8.9 × 10 10 cm −3 . The TeV photons produced from the neutral pion decay will mostly encounter the SSC photons in the energy range 0.35 ≤ γ ≤ 23.6 MeV. The pair production cross section for γ ≥ 0.35 MeV is very small, σ γ γ ≤ 10 −30 cm −2 , which corresponds to a mean free path of λ γ γ ≥ 10 19 cm for the multi-TeV gamma-rays, larger than the outer jet size. So the TeV photons will not be attenuated much due to the e + e − pair production. The parameters used in the photohadronic to fit the 2004 data are summarized in Table 1.  From the kinematical condition given in Eq. (2), the above range of the observed gamma-ray corresponds to the proton energy in the range 16.7 ≤ E p ≤ 210 TeV. Using lep-1, where Γ = D = 14, the seed photon energy lies in the range 0.28 MeV (6.8 × 10 19 Hz) ≤ γ ≤ 3.53 MeV (8.5 × 10 20 Hz), which is again in the tail region of the SSC SED as shown in Fig. 2. A very good fit to the multi-TeV spectrum is obtained by using the EBL-D and EBL-I. The best fit parameters are, respectively, α = 3.1 and A γ = 58.0 for EBL-D and α = 3.2 and A γ = 28.0 for EBL-I, which correspond to a very soft spectrum and the intrinsic spectrum is also soft (between −0.68 and −0.58). Even though these parameters fit well to the observed multi-TeV spectrum, in the low-energy limit the spectrum shoots up to become very high, as shown in Fig. 4, which is certainly not observed by the HESS telescopes. So we can ignore these soft power-law fits to the observed flare data and look for α < 3. This soft power-law problem arises because β = 0.48 is small. So we can use the leptonic model, which has β > 0.48 and we will Fig. 4 The SED of lep-1 is shown along with the power-law fit to the SSC tail region with β = 0.48. The best fit to the flare data of 2010 using EBL-D and EBL-I are shown. A sharp increase in the flux is observed at lower energy. We have also shown the power law with exponential cut-off fit for comparison [30] be able to get α < 3. The time averaged differential energy spectrum of this observation is also fitted with a power law with exponential cut-off having four parameters [30].
In the context of the one-zone leptonic model lep-2 with Γ = D = 25, the observed flare energy range 1.67 TeV (4.0 × 10 26 Hz) ≤ E γ ≤ 20.95 TeV (5.0 × 10 27 Hz) corresponds to the seed photon energy in the range 0.90 MeV (2.17 × 10 20 Hz) ≤ γ ≤ 11.26 MeV (2.72 × 10 21 Hz), which is again in the tail region of the SSC SED as can be seen from Fig. 5. This is fitted with a power law with β = 1.1 and Φ 0 = 4.37 × 10 −9 TeV cm −2 s −1 (the pink line). Here again, we have used EBL-D and EBL-I to fit the 2010 flare data in the photohadronic model, which is also shown in Fig. 5. The best fit parameters are α = 2.6 and A γ = 2.5 for EBL-D and α = 2.6 and A γ = 2.4 for EBL-I, respectively. We note that the observed data is fitted well and, at the same time, the flux decreases towards the low-energy regime as expected with a peak flux of F γ,peak ∼ 2.6 × 10 −10 TeV cm −2 s −1 at E γ ∼ 18 GeV. The EBL-D and the EBL-I corrections to the photohadronic model give practically the same result. The parameters used in the photohadronic model to fit the 2010 data are summarized in Table 1.
Normally, the leptonic models fit the multiwavelength data well. But as shown in Figs. 2 and 5, the leptonic models lep-1 and lep-2 do not fit to the radio (low-energy) data. It is possible that the low-energy behaviors of these models have to be modified to accommodate these radio data. The SSC Fig. 5 The SED of lep-2 [52] is shown along with the power-law fit to the SSC tail region with β = 1.1. The best fits to the flare of 2010 using EBL-D and EBL-I are also shown. For comparison, we have also shown the SED of lep-1, the power-law fit to the SSC tail region with β = 0.48 and the best fit to the flare data of 2004 by the Whipple telescope [36]. The low-energy observed data are taken from Refs. [36,52] photons are produced from the inverse Compton scattering of high-energy electrons with the synchrotron photons which obey the relation SSC 1.3γ 2 e syn , where SSC is the SSC photon energy, syn is the synchrotron photon energy and γ e is the electron Lorentz factor. From the energies of the synchrotron peak and the SSC peak in lep-1, we obtain γ e ∼ 2740. This implies that the 0.35 MeV SSC photon can be produced from the boosting of the ∼ 3.5 × 10 −11 GeV (∼ 8.5 × 10 12 Hz) synchrotron photon, which is in the infrared band. This clearly shows that low-energy photons in the radio band will not affect the prediction of multi-TeV gamma-rays.
Even though lep-1 fits well to the multi-TeV data, in the low-energy regime, the flux increases drastically which is clearly shown in Fig. 4. This behavior is absent with the lep-2 fit. It is to be noted that lep-1 corresponds to the observation during the year 2003-2004 and lep-2 is the recent one of January 2013. So we believe that during each observation period, the photon density distribution in the jet changes and this change in the seed photon changes the spectral behavior of the observed multi-TeV gamma-rays. This clearly shows that almost simultaneous observation in multiwavelength is essential to fit the observed data.
In Fig. 6 we have plotted only the 2010 flare data along with the EBL-D and EBL-I fits and their corresponding intrinsic fluxes. We can observe a minor difference between EBL-D and EBL-I predictions for 0.6 TeV ≤ E γ ≤ 20 TeV. The intrinsic flux is a power law with F γ,in ∝ E −0.  shows that the different spectral shapes of the observed events are solely due to the diversity in the shape of the seed photon density distribution (particularly in the SSC tail region) even though we have the same acceleration mechanism (α 2.6) of protons in the blazar jet.

Summary
We used the photohadronic scenario complemented by two EBL models (EBL-D and EBL-I) to interpret the multi-TeV flares observed in April 2004 and February 2010. These EBL models are equally good to explain the observed data, which can also be seen from the comparison of the attenuation factors in Fig. 1. The photohadronic scenario with both EBL models can fit these multi-TeV spectra very well and it is observed that the intrinsic spectrum of each epoch is different even though the Fermi accelerated high-energy protons have almost the same spectral index α 2.6. This difference in intrinsic spectra is attributed to the change in the spectral index of the seed photon in the SSC tail region of the jet in different epochs of flaring. The same photohadronic models were earlier employed to explain the multi-TeV emission from other HBLs [57,59]. We suggest that the same flaring mechanism may be acting in all HBLs. The differences between the flaring behavior of HBL blazars and in one blazar from flare to flare is determined by the changes in flux and spectrum of the seed photons in the SSC region. It is important to note that simultaneous or quasi-simultaneous observa-tion of the tail region of the SSC SED with TeV observation would be most useful to constrain the photohadronic model.