EBL effect on the observation of multi-TeV flaring of 2009 from Markarian 501

Markarian 501 is a high-peaked BL Lacertae object and has undergone many major outburst since its discovery in 1996. As a part of the multiwavelength campaign, in the year 2009 this blazar was observed for 4.5 months from March 9 to August 1 and during the period April 17 to May 5 it was observed by both space and ground based observatories covering the entire electromagnetic spectrum. A very strong high energy $\gamma$-ray flare was observed on May 1 by Whipple telescope in the energy range 317 GeV to 5 TeV and the flux was about 10 times higher than the baseline flux. We use the photohadronic model complimented by the extragalactic background radiation (EBL) correction to this very high state flare and have shown that the EBL plays an important role in attenuating the very high energy flux even though Markarian 501 is in the local Universe.


Introduction
Blazars are a subclass of AGN and the dominant extra galactic population in gamma rays [1]. These objects show a rapid variability in the entire electromagnetic spectrum and have non-thermal spectra which implies that the observed photons originate within the highly relativistic jets oriented very close to the observers line of sight [2]. Due to the small a e-mail: sarira@nucleares.unam.mx b e-mail: vladimir@ciencias.unam.mx c e-mail: luis.miranda@correo.nucleares.unam.mx d e-mail: albertoros4@ciencias.unam.mx viewing angle of the jet, it is possible to observe the strong relativistic effects, such as the boosting of the emitted power and a shortening of the characteristic time scales, as short as minutes [3,4]. Thus these objects are important to study the energy extraction mechanisms from the central supermassive black hole, physical properties of the astrophysical jets, acceleration mechanisms of the charged particles in the jet and production of ultra high energy cosmic rays, very high energy γ -rays and neutrinos.
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 self-produced synchrotron photons [5,6]. Depending on the location of the first peak, blazars are often sub-classified into low energy peaked blazars (LBLs) and high energy peaked blazars (HBLs) [7]. In LBLs, the first peak is in the nearinfrared/optical energy range and the second peak is around GeV energy range. For HBLs, the first peak is in the UV or X-rays range and the second peak is in the GeV-TeV energy range. The above scenario is called leptonic model and is very successful in explaining the multi wavelength emission from blazars and FR I galaxies [8][9][10][11].
Flaring seems to be the major activity of the blazars, which is unpredictable and switches between quiescent and active states involving different time scales. While in some blazars a strong temporal correlation between X-ray and multi-TeV γray has been observed, outbursts in some others have no low energy counterparts (orphan flaring) [12,13] and the explanation of such extreme activity needs to be addressed through different mechanisms. It is also very important to have simultaneous multiwavelength observations of the flaring period to constrain different theoretical models of emission in different energy regimes.
The TeV photons of the flare can interact with the background soft photons in the jet to produce e + e − pairs. However, production of the lepton pair within the jet depends on the size of the emitting region and the photon density in it. Also the required target soft photon threshold energy ε γ ≥ 2m 2 e /E γ is needed. It has been observed that the jet medium is transparent to pair production where the optical depth is very small [14,15]. Also the TeV photons on their way to Earth can interact with the extragalactic background light (EBL) to produce the lepton pair [16][17][18][19][20].

Markarian 501
Markarian 501 (Mrk 501) (RA: 251.46 • , DEC: 39.76 • ) is a HBL at a redshift of z = 0.034 (local Universe) and one of the brightest extragalactic sources in X-ray/TeV sky [14]. It is also the second extragalactic object (after Markarian 421) identified as a very high energy (VHE) emitter by the Whipple telescope in 1996. Since its discovery, the multiwavelength correlation of Mrk 501 has been studied intensively and during this period it has undergone many major outbursts on long time scales and rapid flares on short times scales mostly in the X-rays and TeV energies [21][22][23][24][25][26][27][28][29][30][31]. It has been observed that, during these outbursts, both peaks have shifted to higher energies and during the most extreme case the synchrotron peak ∼ keV range has shifted above 200 keV [1]. Due to the low sensitivity of the previous generation instruments, Mrk 501 was primarily observed in VHE band during the outbursts. However, later on it was observed in all the wave bands. In the year 2009, Mrk 501 was observed as a part of large scale multiwavelength campaign covering a period of 4.5 months (from March 9 to August 1, 2009) [32]. The scientific goal of this extended observation was to collect a simultaneous, complete multifrequency data set to test the current theoretical models of broadband blazar emission mechanism. Also this will help to understand the origin of high energy emission from blazars and the physical mechanism responsible for the acceleration of the charged particles in the relativistic jets. Between April 17 to May 5, Mrk 501 was observed by both space and ground based observatories, covering the entire electromagnetic spectrum including even the variation in optical polarization [32]. A very strong VHE flare was detected first by the Whipple telescope on May 1 and 1.5 h later with VERITAS. Both of these telescopes continued simultaneous observation of this VHE flare until the end of the night. The detected flux was enhanced by a factor of ∼10 above the average baseline flux (3.9 × 10 −11 ph cm −2 s −1 ). A dramatic increase in the flux by a factor ∼4 in 25 min and a falling time of ∼50 min was observed. The flux measured at lower energies before and after the VHE flare did not show any significant variation. But Swift-XRT (in X-ray) and UVOT (in optical) did observe a moderate flux variability [32]. Also both Whipple and VERITAS did observe statistically significant variation in VHE band. Using the one-zone SSC model, the average SED of this multiwavelength campaign of Mrk 501 is interpreted satisfactorily.
Our aim here is to use the photohadronic model of Sahu et al. [15,[33][34][35][36] and the EBL model of Dominguez et al. [19] to interpret the observed very strong VHE flare data of May 1 and the long outburst observed by HEGRA telescopes in 1997. We found that both these flares can be well explained with this model.

TeV flaring model
The photohadronic model of Sahu et al. [15,35,36] rely on the standard interpretation of the leptonic model to explain both low and high energy peaks by synchrotron and SSC photons, respectively, as in the case of any other AGNs and blazars. Thereafter, it is proposed that the flaring occurs within a compact and confined volume of radius R f inside the blob of radius R b (R f < R b ) [15] (henceforth implies the jet comoving frame). Both the internal and the external jets are moving with the same bulk Lorentz factor and the Doppler factor D as the blob (for blazars D). In normal situation within the jet, we consider the injected spectrum of the Fermi accelerated charged particles having a power-law spectrum dN /dE ∝ E −α with the power index α ≥ 2. But in the flaring region the injected proton spectrum is a power-law supplemented with an exponential decay factor given as Here the high energy proton has the cut-off energy E p,c . The high energy protons will interact in the flaring region where the comoving photon number density is n γ,f to produce the -resonance. Subsequently the -resonance decays to charged and neutral pions and the further decay of neutral pions to TeV photons gives the multi-TeV SED. The n γ,f is much higher than the rest of the blob n γ (non-flaring) i.e. n γ,f (ε γ ) n γ (ε γ ). There is no direct way to estimate the photon density in the inner jet region as it is hidden. For simplicity we assume the scaling behavior of the photon densities in the inner and the outer jet region as which assumes that the ratio of photon densities at two different background energies ε γ 1 and ε γ 2 in flaring and non-flaring states remains almost the same. While the photon density in the outer region can be calculated from the observed flux, using Eq. (2) we can express the n γ,f in terms of n γ . The π 0 -decay TeV photon energy E γ and the target SSC photon energy ε γ in the observer frame are related through, The observed TeV γ -ray energy and the proton energy E p are related through The optical depth of the -resonance process in the inner jet region is given by where the resonant cross section is σ ∼ 5 × 10 −28 cm 2 .
The efficiency of the pγ process depends on the physical conditions of the interaction region, such as the size, the distance from the base of the jet, the photon density and their distribution in the region of interest.
In the inner region we compare the dynamical time scale t d = R f with the pγ interaction time scale t pγ = (n γ,f σ K pγ ) −1 to constrain the seed photon density so that multi-TeV photons can be produced. For a moderate efficiency of this process, we can assume t pγ > t d and this gives τ pγ < 2, where the inelasticity parameter is assigned with the usual value of K pγ = 0.5. Also by assuming the Eddington luminosity is equally shared by the jet and the counter jet, the luminosity within the inner region for a seed photon energy ε γ will satisfy (4π n γ,f R f ε γ ) L Edd /2. This puts an upper limit on the seed photon density as From Eq. (6) we can estimate the photon density in this region. In terms of SSC photon energy and its luminosity, the photon number density n γ is expressed as where η is the efficiency of SSC process and κ describes whether the jet is continuous (κ = 0) or discrete (κ = 1). In this work we take η = 1 for 100% efficiency. The SSC photon luminosity is expressed in terms of the observed flux ( SSC (ε γ ) = ε 2 γ dN γ /dε γ ) and is given by Using Eqs. (7) and (8) we can simplify the ratio of photon densities given in Eq. (2) to The γ -ray flux from the π 0 decay is deduced to be The exponential factor in the power spectrum in Eq. (1) is responsible for the decay of the VHE flux, and falls faster for Here E c is the γ -ray cut-off energy corresponding to E p,c . The EBL effect also attenuates the VHE flux by a factor of e −τ γ γ , where τ γ γ is the optical depth which depends on the energy of the propagating VHE γ -ray and the redshift z of the source. So there is a competition between the exponential cut-off and the EBL effect. A 6 TeV photon was observed during the 4.5 months campaign and the attenuation factor e −τ γ γ for this photon is about 0.4-0.5 [14]. So attenuation by EBL is significant for multi-TeV γ -rays even for nearby sources [37]. Here we would like to study the effect of EBL on the strongest VHE flare of May 1 and compare with the exponential cut-off scenario.
Including the EBL effect, the relation between observed flux F γ and the intrinsic flux F int is given as Then the EBL corrected observed multi-TeV photon flux from π 0 -decay at two different observed photon energies E γ 1 and E γ 2 can be expressed as where we have used The SSC at different energies are calculated using the leptonic model. Here the multi-TeV flux is proportional to E −α+3 γ and SSC (ε γ ). In the photohadronic process ( pγ ), the multi-TeV photon flux is expressed as Both ε γ and E γ satisfy the condition given in Eq. (3) and the dimensionless constant A γ is given by Comparing Eqs. (11) and (14), the intrinsic flux F int is given as Using Eq. (14), we can calculate the EBL corrected multi-TeV flux where A γ can be fixed from observed flare data. We can calculate the Fermi accelerated high energy proton flux F p from the TeV γ -ray flux through the relation [35] The optical depth τ pγ is given in Eq. (5). For the observed highest energy γ -ray E γ corresponding to a proton energy E p , the proton flux F p (E p ) will be always smaller than the Eddington flux F Edd . This condition puts a lower limit on the optical depth of the process and is given by From the comparison of different times scales and from Eq. (18) we will be able to constrain the seed photon density in the inner jet region.

Results
The average broadband SED of Mrk 501 is modeled using the standard one-zone leptonic model [32]. The emission takes place from a spherical blob of size R b , which moves down the conical jet with a bulk Lorentz factor and a Doppler factor D. The emission region is filled with an isotropic and non-thermal population of electrons and a randomly oriented magnetic field B . To interpret the VHE flare of May 1, 2009, we use the parameters of the one-zone leptonic model which fits reasonably well the average SED and the parameters are shown in Table 1.  Fig. 1 The average SED of Mrk 501 is shown in all the energy bands which are taken from Ref. [32]. The SED of low state (MJD 54936-54951; blue squares) and high state (MJD 54952-55; red circles) of the 3-week period are shown. The leptonic model fit to the low state (blue curve) and high state (red curve) are also shown. The blue dotted curve corresponds to the optical emission from the host galaxy. The black curve is the photohadronic fit to the Whipple very high state data (red circles) The observed VHE flare of May 1, 2009 by the Whipple telescope was in the range ∼ 317 GeV ≤ E γ ≤ 5 TeV. In the context of photohadronic scenario, this range of E γ corresponds to the Fermi accelerated proton energy in the range 3.2 TeV ≤ E p ≤ 50 TeV. So protons in this energy range will interact with the background SSC photons in the energy range 13.6 MeV(3.29 × 10 21 Hz) ≥ ε γ ≥ 0.86 MeV(2.1 × 10 20 Hz) to produce the -resonance and subsequent decay of it will produce both γ -rays and neutrinos through neutral and charged pion decay. Also the above range of ε γ lies in the beginning of the SSC spectrum and in this range of energy the sensitivity of the currently operating instruments are not good enough to detect Mrk 501. However, from the multiwavelength campaign the average SED is fitted very well (Fig. 1) and we use this low energy flux in the photohadronic model to calculate the observed flux. Also to account for the contribution of the EBL on the multi-TeV photons we consider the EBL model by Dominguez et al. The EBL models of Dominguez et al. [19] and Franceschini et al. [20] are widely used to constrain the imprint of EBL on the propagation of VHE γ -rays by Imaging Atmospheric Cherenkov Telescopes (IACTs). The normalization constant A γ given in Eq. (15) can be calculated from the observed flare data.
The multi-TeV flaring from blazars has an exponential fall which is conventionally modeled as shown in Eq. (1). The cut-off energy E c is a free parameter and depends on some unknown mechanism. On the other hand, the diffuse background radiation also attenuates the high energy γ -rays as a consequence of the lepton pair production. Here instead of the additional exponential cut-off, we take into account the effect of EBL to deplete the intrinsic VHE flux. A very good  Fig. 2 The black curve is the hadronic model fit which includes the EBL attenuation using the EBL model of Dominguez et al. [19] to the Whipple very high state flare data (red filled circles) of Mrk 501 and the red continuous curve is the intrinsic flux in the same model. For comparison we have also shown the Whipple fit to the data (dashed curve) and the exponential fit (dashed dotted curve). We have also shown the HEGRA observation of the outburst during 1997: conventional analysis (open circles) [30] and new analysis with improved energy resolution (blue filled squares) [31] fit to the Whipple very high state data of May 1 is obtained for α = 2.4 and A γ = 89 where the EBL corrected flux is considered. We observed that the EBL correction to the VHE γ -ray is small but not insignificant (black curve in Fig.  2); there is a faster fall above 10 TeV. We have also shown the intrinsic flux (red curve in Fig. 2) to demonstrate the difference. For comparison we have fitted the data with an exponential cut-off function (dashed dotted curve) and the best fit is obtained for α = 2.6, E c = 30 TeV and A γ = 66. Also we have shown the Whipple fit (dashed curve) for comparison, where it is fitted by the function dN γ /dE γ = 9.1×10 −7 (E γ /1 TeV) −2.1 ph m −2 s −1 TeV −1 . It is observed that the very high state data of Whipple fits very well with the above three scenarios. However, above 5 TeV, both the EBL corrected fit and the exponential fit differ from the Whipple fit. Again the EBL fit and the exponential fit differ above 10 TeV and the former one falls faster than the latter as can be seen from Fig. 2. Even though all these fit very well with the Whipple very high state data, we observe deviation in the VHE limit. So observation of the VHE flux above 10 TeV will be a good test to constrain the EBL effect on the propagation of VHE γ -rays. In Fig. 1, we also plotted the Whipple very high state data and our model fit (black curve) along with the complete SED.
The high energy protons will be accompanied by high energy electrons and these electrons will emit synchrotron photons in the energy range ∼ 10 19 to ∼ 10 23 Hz when encountering the magnetic field of the jet. This energy range photons lie in between the high energy end of the synchrotron spectrum and the low energy tail of the SSC spectrum, thus may not be observed due to their low flux in this region. These high energy electrons will also emit SSC photons and their energy is given by E IC ∼ γ 2 e ε syn . As discussed before, in the flaring state, in general, the flux of the two opposing jets can be as high as F Edd /2. However, the highest energy protons with E p = 50 TeV must have a flux F p < F Edd /2 0.8 × 10 −7 erg cm −2 s −1 . This constraint translates into τ pγ > 0.04, which corresponds to n γ,f > 1.5 × 10 10 cm −3 in the inner jet. However, the hidden jet lies between R s (Schwarzschild radius) and R b . As one representative value we take R f 5 × 10 15 . From Eq.(6) the seed photon density for ε γ = 0.86 MeV satisfies the inequality n γ,f < 5.1 × 10 10 cm −3 , which translates to the optical depth to be constrained as τ pγ < 0.13. So the optical depth lies in the range 0.04 < τ pγ < 0.13 and this corresponds to the range of photon density in the inner jet region as 1.5× 10 10 cm −3 < n γ,f < 5.1× 10 10 cm −3 , which shows that the photon density in this region is high. Due to the adiabatic expansion of the inner blob, the photon density will be reduced to n γ and also the optical depth τ pγ 1. The energy will dissipate once these photons cross into the bigger outer cone. This will drastically reduce the -resonance production efficiency from the pγ process.
Between March 16 and October 1, 1997, Mrk 501 was in a flaring state which was monitored in TeV γ -rays with the HEGRA stereoscopic system of imaging atmospheric Cherenkov telescopes (IACTs). During this long outburst period (a total exposure time of 110 h) more than 38,000 TeV photons were detected and the energy spectrum of the source was well beyond 10 TeV. Although there was dramatic flux variation during the observation period, the shapes of the daily γ -ray spectra remained essentially stable. A timeaveraged energy spectrum of this observation period was fitted with a power-law accompanied with an exponential cutoff [30]. However, the same data was reanalyzed with an improved energy resolution [31] and it was found that except for the highest energy the two analyses were in very good agreement. At the highest energy the spectrum was found to be much steeper than the conventional analysis. In Fig.  2, along with the 2009 flare spectrum, we have also shown the conventional analysis and the improved energy resolution analysis of the 1997 outburst observed by HEGRA telescopes. Our model fits well with the 1997 flare data and we found a steeper slope of the spectrum beyond 10 TeV which is perfectly compatible with the highest energy point of the reanalysis result of 1997.

Conclusions
The VHE flare of May 1, 2009 observed by the Whipple telescope can be explained very well through photohadronic model supplemented with the EBL correction. Previously, the decay of the VHE flare can be explained through the exponential fall of the flux which introduces an additional free parameter, the cut-off energy. However, here, the EBL corrected VHE flux automatically falls exponentially without any additional free parameter and fits very well with the Whipple very high state data. For comparison we have also shown the Whipple fit as well as the exponential fit. All these three curves fit very well with the VHE flare data. However, we have shown that their behaviors differ in the high energy limit. The HEGRA telescopes observed an outburst in Mrk 501 during 1997 and the energy spectrum was above 10 TeV. We found that the EBL corrected photohadronic model fits very well beyond 10 TeV with the improved energy resolution analysis where the spectrum has a steeper slope.