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 stud- ied intensively and was observed by various Cherenkov tele- scope 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.


I. 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 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 about 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 evolve 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 AGN and the dominant extra galactic population in γ-rays [13].
These objects show rapid variability in the entire electromagnetic spectrum and have nonthermal 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 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 believed to be due to the synchrotron self Compton (SSC) scattering of the high energy electrons with their selfproduced synchrotron photons [15,16]. As leptons (e ± ) are responsible for the production of the SED, this is called the leptonic model and in general 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.
During April 2004, 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 X-ray 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 pattern 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 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 six month long multi-instrument campaign, the MAGIC telescopes Mrk 421 being the nearest HBL, the VHE gamma-rays 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 the 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 high flux which is important for our study. The flare of 2004 was analyzed by Sahu et al. [53] where 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.

II. PHOTOHADRONIC SCENARIO
In the photohadronic scenario, the Fermi accelerated protons having a power-law spec- interact with the background photons to produce the ∆-resonance which subsequently decays to gamma-rays via intermediate neutral pion and to neutrinos through 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 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 the scaling behavior of the photon densities in the inner and the outer jet regions which essentially means, the spectra of both 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 power-spectrum 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 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 evade 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 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 Γ. 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. By expressing 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 and 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 low energy 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  1 presents the attenuation factor for these two models as functions of observed gamma-ray energy E γ for z = 0.031. The Fig. 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 than decreases until the models intersect at 6 TeV. At higher energies the models diverge and 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 to the Earth as ultra high energy cosmic rays (UHECRs). The charged pions produced from the photohadronic process will decay to neutrinos which can in principle be detected on Earth [60][61][62], however, the flux is too low. (3) with Φ 0 = 6.0 × 10 −10 erg cm −2 s −1 and β = 0.48 (red curve).

III. TWO FLARING EPISODES
As discussed in the introduction, several epochs of major multi-TeV emission/flaring are observed from Mrk 421 and below we analyze two of these multi-TeV flares observed by different Cherenkov telescope arrays and use photohadronic model to interpret these flaring events. As input to the photohadronic model we use the leptonic SSC SED fitted to the multiwavelength observations.  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 γ < 4T eV both these 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 EBL-D fit is better than the EBL-I. Again above 20 TeV both the EBL models behave the same. 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 compared to 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 both e −Eγ /Ec 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.
The intrinsic flux in EBL-D and EBL-I are also shown. Both are almost the same and having power-law behavior with F γ,in ∝ E −0.18 γ . Even though all these models fit quite well to the observed data below 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 energy cut-off scenario is necessary or not.
In ref. [53], by comparing expansion time scale, interaction time scale of pγ interaction and the high energy proton luminosity to be smaller than the Eddington luminosity in the inner jet region R f 3 × 10 15 cm (size of the outer jet is R b 0.7 × 10 16 cm), the range of 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 as 1.3 × 10 10 cm −3 < 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 M eV ≤ γ ≤ 23.6 M eV . The pair production cross section for γ ≥ 0.35 M eV 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 I. 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].

and
A γ = 28.0 for EBL-I which correspond to very soft spectrum and the intrinsic spectrum is also soft (between -0.68 to -0.58). Even though these parameters fit well to the observed multi-TeV spectrum, in the low energy limit the spectrum shoots up very high as shown in Fig. 4 which is certainly not observed by HESS telescopes. So we can ignore these soft power-law fit to the observed flare data and to look for α < 3. This soft power-law problem arises because β = 0.48 is small. So we can use the leptonic model which have β > 0.48 and will 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  Table   I.
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 behavior of these models have to be modified to accommodate these radio data. The SSC 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.7 γ . Comparison of F γ,in of both 2004 and 2010 multi-TeV flaring shows that different spectral shape 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.

IV. SUMMARY
We used the photohadronic scenario complemented by two EBL models (EBL-D and comparison of the attenuation factors in Fig. 1. The photohadronic scenario with both the 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 attributes 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 observation of the tail region of the SSC SED with TeV observation would be most useful to constrain the photohadronic model.