Very hard gamma-ray emission from a flare of Mrk 501
Aims: We investigate the peculiar properties of a large TeV -ray flare from Mrk 501 detected during the 2009 multiwavelength campaign.
Methods: We identify the counterpart of the flare in the Fermi LAT telescope data and study its spectral and timing characteristics.
Results: A strong order-of-magnitude increase of the very-high-energy -ray flux during the flare was not accompanied by an increase in the X-ray flux, so that the flare was one of the “orphan”-type TeV flares observed in BL Lacs. The flare lasted about 1 month at energies above 10 GeV. The flaring source spectrum in the 10-200 GeV range was very hard, with a photon index , harder than that observed in any other blazar in the -ray band. No simultaneous flaring activity was detected below 10 GeV. Different variability properties of the emission below and above 10 GeV indicate the existence of two separate components in the spectrum. We investigate possible explanations of the very hard flaring component. We consider, among others, the possibility that the flare is produced by an electromagnetic cascade initiated by very-high-energy -rays in the intergalactic medium. Within such an interpretation, peculiar spectral and temporal characteristics of the flare could be explained if the magnetic field in the intergalactic medium is of the order of G.
Variable -ray emission from BL Lacs is most commonly interpreted as Doppler-boosted inverse Compton (IC) emission by high-energy electrons in relativistic jets aligned with the line of sight. In the Synchrotron-self-Compton (SSC) model, the seed photons for the IC scattering are provided by the lower energy (radio-to-X-ray) synchrotron emission produced by the same electrons. Since both the synchrotron and IC emissions are produced by the same electrons, the variability of the -ray emission is expected to be accompanied by a similar variability in the X-ray emission. However, this is not always the case. Strong “orphan” -ray flares in the TeV energy band are sometimes observed without accompanying X-ray synchrotron flares (Krawczynski et al., 2004; Daniel et al., 2005). The nature of these orphan flares is not clear. Different mechanisms have been proposed to explain these orphan flares in the context of leptonic (Kusunose & Takahara, 2006) or hadronic (Böttcher, 2005) models of blazar emission.
In the framework of the SSC model for the -ray emission, the emission spectrum is expected to be characterized by a photon index ( is defined as the slope of the differential spectrum ) in the energy range where electron cooling via synchrotron and/or IC energy losses is efficient. The limiting value of is considered as a lower bound on in the analysis of VHE -ray flux attenuation through interactions with Extragalactic Background Light (EBL) (Aharonian et al., 2006, 2007a, 2007b). At the same time, harder than spectra could, in principle, be produced in the VHE -ray band as a result of absorption of -rays in pair production inside the source (Aharonian et al., 2008), due to the presence of low-energy cut-offs or hardenings in the electron distribution (Katarzynski et al., 2006a; Böttcher et al., 2008) or the formation of narrow energy distribution of emitting particles in the source (Aharonian et al., 2002). The observation of blazars with intrinsically hard -ray spectra, harder than , would have important consequences for EBL studies and for the physical models describing the -ray emission.
In this letter we consider the spectral and timing characteristics of a large -ray flare of Mrk 501 observed during a multi-wavelength campaign in 2009 (Abdo et al., 2011). The flare was remarkable both because it was an “orphan” flare (the factor of 10 increase in the VHE -ray flux was not accompanied by a similar simultaneous increase in the X-ray flux) and because the spectrum of the -ray emission during the flare was much harder than the spectrum of the quiescent emission from the source. We show that the spectral and timing characteristics of the source during the flare suggest the existence of two distinct components in the -ray emission: a non-flaring one, characterized by a softer photon index , close to the time-averaged photon index of the source, and a flaring one, with an extremely hard spectrum . We discuss the possible physical origins of this hard flaring component.
2 Fermi data selection and data analysis
For our analysis we used the publicly available data collected by the Large Area Telescope (LAT) on board the Fermi satellite (Atwood et al., 2009) between August 4, 2008 and March 15, 2011. We processed the “PASS 6” data using the Fermi Science Tools111http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/. We filtered the entire data set with gtselect and gtmktime tools following the recommendations of the Fermi team222http://fermi.gsfc.nasa.gov/ssc/data/analysis/ and retained only events belonging to the class 4, which are most confidently identified with -rays, as described by Abdo et al. (2010a). To estimate the flux from the photon counts in 30-400 GeV band we used the gtexposure tool. The spectral analysis was done using the unbinned likelihood method (Mattox et al., 1996), taking into account all the sources from the 1-st year Fermi catalog (Abdo et al., 2010b) within from Mrk 501 as well as an additional bright source at the position of a quasar B1638+3952, not listed in the catalog.
Timing analysis. The long-term lightcurves of Mrk 501 for the whole period of observations are shown in Fig. 1 for three different energy bands. For the highest energy band 30-400 GeV, the lightcurve was produced using the aperture photometry method, using time bins with an equal signal-to-noise ratio . To estimate the diffuse background level in this energy bin we used a nearby source-free region of radius centered on Ra, Dec. In the lower energy bands the binning was homogeneous in time, with time bins of 30 days, produced using the likelihood analysis.
A flare with a flux increase by a factor of was readily identified in the lightcurve in the 30-400 GeV energy band. The flare period includes the Veritas flare reported by Abdo et al. (2011). The moment of onset of the Veritas flare is marked in Fig. 1 by the vertical dashed line. The flare could not be identified in the lower energy lightcurves shown in Fig. 1. In principle, if the flux enhancement during the flare was moderate and the flare duration was much shorter than 1 month, the flare signal would have been “diluted” in the lightcurves binned using the 30 day binning. We have verified that reducing the size of the time bins down to several days (the time scale on which the source could be significantly detected at lower energies) does not reveal the flare in the lower energy lightcurves.
The duration of the flare in the VHE band was not constrained by the Veritas observations (Abdo et al., 2011), since their observations stopped 3 days after the onset of the flare, when the source was still in the flaring state. In spite of the low statistics of the VHE signal in LAT, the LAT data could be used to estimate the duration of the flare. Fig. 2 shows the time delay between subsequent photons in the 100-400 GeV band. During the flare the time delay between subsequent photons dropped below 10 days for about one month. During the rest of the Fermi exposure the GeV photon count rate from the source was about 1 photon per 100 days. To estimate the duration of the flare we divide all the 100-400 GeV photons into “flare signal” and “background”. The “background” count rate is ph/day. The time interval d following MJD 54953 has 5 photons, four of which we ascribe to the “flare signal” (the probability for a background photon to come in the d interval is 0.5). The d time scale can be considered as an estimate of the flare duration. The uncertainty of this estimate is by a factor of . Indeed, if the flare duration was d, instead of d, the probability for all four “flare signal” photons to come in the first half of the time interval would be , so that the d flare duration is excluded at confidence level. The waiting time between subsequent photons drops down to d at the peak of the flare. Just before the peak no photons were detected during the previous d. This constrains the rise time of the flare. Assuming that the flux level comparable to the peak flux was sustained over a rise time interval , one could estimate an upper limit on . The probability that no photons are detected during a time interval d before the peak is %, so that d is ruled out at % confidence level.
Spectral analysis. The timing behaviour of the source below and above 10 GeV is different: the flare appears most prominently at the highest energies ( GeV) and is not detectable at low energies. Taking this into account we have analyzed the low- and high-energy part of the LAT energy band separately. We have extracted the spectra in the 0.3-10 GeV and 10-400 GeV bands for the time period of the flare, between and d, corresponding to the highest count rate of GeV photons, see Fig. 2. We assume that in each energy band the spectrum is well described by a simple powerlaw .
The results of the likelihood analysis of the spectrum in the two energy bands are shown in Fig. 3. One can see that the slopes of the spectrum below and above 10 GeV are significantly different. Below 10 GeV the likelihood analysis gives the slope , while above 10 GeV the spectrum is much harder, with . We note that the highest energy photon from the source detected by LAT during the flare period had an energy GeV. Taking this into account we show the powerlaw model for the source spectrum only up to 200 GeV in Fig. 3.
The fact that most of the highest energy photons from the source come in a narrow time window around the time of the flare was already noticed by Abdo et al. (2011). This peculiarity of the highest energy photon timing was found to lead to the hardest source spectrum during the 30 d flare period, with an average slope of in the 0.2-400 GeV band. Our analysis of the spectral variability details during the flare suggests that the flaring activity appears only at the highest energies. Taking this into account, we fit the spectrum extracted from narrow energy bins (black thick data points in Fig. 3) with a broken powerlaw, instead of the simple powerlaw model. We find that the use of the broken powerlaw model improves the quality of the fit, with an F-test probability of the chance fit improvement at the level of 8%. The best fit value of the break energy is found fitting the data by the broken powerlaw model is GeV. The quality of the data is not sufficient for limiting the break energy from above, only a lower bound GeV could be found (Fig. 4). If the break energy is left free during the fitting, the range of allowed values of is wider, so that the lower bound on .
The difference in the spectral slopes and in the timing behaviour at low and high energies point to the presence of two separate components in the emission from the source. The “soft” component which appears below GeV is steady while the “hard” component, visible mostly above 30 GeV, is variable during the flare. Two additional arguments which support the hypothesis for the existence of a slow-variable soft component are the closeness of the slope of the spectrum in the 0.3-10 GeV band to the time-average slope (Abdo et al., 2010b) and the consistency of a powerlaw extrapolation of the low-energy spectrum to the VHE band with the quiescent spectrum of the source measured by Veritas and MAGIC (shown by the dashed data points in Fig. 3).
3 Discussion: origin of the hard -ray emission
The photon index of the high-energy component is extremely hard. In particular, it is harder than , often adopted as a lower bound on the photon index, based on simple SSC models for the broad band emission from blazars. Such a limit on the photon index of blazars underlies the derivation of constraints on the EBL density (see e.g. Aharonian et al. (2006, 2007a, 2007b)). The detection of a hard spectral component in the flare of Mrk 501 has, therefore, important implications for the modeling of blazar spectra and the EBL. In the following sections we discuss possible mechanism(s) for the production of a very hard spectral component.
3.1 Mechanisms intrinsic to the source
Several models which could explain very hard intrinsic blazar spectra in the -ray band have already been proposed.
3.1.1 Internal absorption origin
Aharonian et al. (2008) considered the possibility for the formation of a hard -ray blazar spectra through the absorption of -rays in the dense soft photon background, with a narrow (e.g. thermal) spectrum. -rays with energies are most efficiently absorbed through interactions with soft photons of energy eV. Optical depths would be sufficient to produce a hard spectrum above the energy GeV. The observation of a hard spectral component in Mrk 501 would imply, therefore, the existence of a dense field of soft photons with energies eV in a compact -ray emission region. If the soft photon field has a thermal spectrum, its temperature should be K, much higher than the typical temperature of accretion disks in Active Galactic Nuclei (AGN) accreting in a radiatively-efficient way. At the same time, this temperature range might not be unrealistic for radiatively-inefficient accretion flows, like those found in the Fanaroff-Riley I radio galaxies such as the nearby radio galaxy M 87 (Neronov & Aharonian, 2007). A potential problem for the “absorption feature” explanation of the hard spectrum is that if the soft photon field has a narrow energy distribution, the suppression of the -ray signal should become small at energies much lower than GeV. This means that the flaring component of the spectrum should reappear at the low-energy end of the LAT energy range, which is not the case. This difficulty could, however, be overcome if the soft photon distribution has a broad spectrum extending to energies much higher than eV. This might also be typical for the radiatively inefficient accretion flows where the emission in the 0.1-1 keV range is supposed to be produced via inverse Compton (IC) by the accretion flow electrons with a quasi-thermal distribution (Neronov et al., 2008).
3.1.2 IC from hard spectrum electron origin
Katarzynski et al. (2006a) and Böttcher et al. (2008) proposed that a hard spectra could be produced via the IC mechanism if the spectrum of IC emitting electrons is hard enough. In particular, Katarzynski et al. (2006a) consider the possibility of a low energy cut-off in the electron spectrum. In this case the spectrum of IC emission in the SSC model could have a photon index as low as . Böttcher et al. (2008) consider the possibility of a hard spectrum for the case when the seed photons for IC scattering are provided by the Cosmic Microwave background (CMB) and the electron spectrum has a slope with . Each of these models carries with it a potential problem with relation to source variability.
On the one hand, in the model of Katarzynski et al. (2006a), an immediate effect of cooling is that an type spectrum forms at energies below the low-energy cut-off, giving rise to a photon spectra. Furthermore, if the synchrotron cooling of the highest energy flux occurs in the Klein-Nishina regime, even softer, , photon spectra are obtained. The cooling times of electrons emitting above 10 GeV are days. Thus only fields of G would allow the hard spectrum IC emission to live for long enough to explain the persistence of the hard spectrum during the 30 day flare. Furthermore, the acceleration time of 10 GeV electrons in such fields (Rieger et al., 2007), days, leaves the question as to the origin of the rise time of the flare unanswered. A possible way out of this difficulty is that the flare rise time is determined not by the acceleration time scale, but by the dynamical time scale of the system (this would imply, however, that there is a pre-existing population of accelerated electrons). Cooling times shorter than the synchrotron cooling time could also be achieved if an additional cooling mechanism is present, such as adiabatic cooling (Finke et al., 2008).
Similarly, the problem of the cooling time for the model of Böttcher et al. (2008) limits its applicability to short flares. For this model, the cooling time of electrons emitting at GeV energies is very long, yr, leaving this mechanism inapplicable for the description of a short flaring episode.
SSC type models embedded within the stochastic acceleration framework, however, do not suffer so badly from these timescale problems (Katarzynski et al., 2006b). The steady state spectra in this scenario, obtained when stochastic acceleration is balanced by radiative losses, can indeed be hard enough to explain the observed hard flare of Mrk 501. Furthermore, for these models, the longer acceleration time for 10 GeV electrons in 0.1 G magnetic fields, (100 km s days, can more naturally explain the rise-time of the flaring episode. The Alfven speeds motivated by this scenario, however, require rather dense astrophysical environments with cm (O’Sullivan et al., 2009).
3.1.3 Hadronic models
Hard intrinsic spectra in the -ray range can also arise in hadronic models (Rachen&Meszaros, 1998; Aharonian, 2000; Mücke & Protheore, 2001; Reimer et al., 2004; Zacharopoulou et al., 2011). Variability times as short as s can originate from the synchrotron cooling time of protons in strong G magnetic fields, if the observed -ray emission has proton and/or muon synchrotron origin. A hard spectrum, within this framework, can form as the result of intrinsic absorption of the -ray flux on soft photon fields in the source (Zacharopoulou et al., 2011), in a way similar to the leptonic models with intrinsic absorption. Alternaviely, the spectrum of accelerated protons might be intrinsically hard and have a pileup at the highest energies (Aharonian, 2000). Finally, if significant contribution to the source flux is produced through muon synchrotron emission, the hard spectrum could arise by muons with an energy below a certain critical energy decaying before they lose their energy through synchrotron emission (Rachen&Meszaros, 1998; Reimer et al., 2004).
3.1.4 Vacuum gap acceleration origin
The problem of the absence of cooling of electrons encountered in the models of Katarzynski et al. (2006a) and Böttcher et al. (2008) disappears if one assumes that the observed -ray emission comes directly from the electron acceleration region. In this case the hard spectrum of electrons is maintained by the balance between the energy gain and energy loss rates. Such a possibility is considered e.g. in models describing particle acceleration in black hole magnetospheres (Beskin et al., 1992; Levinson, 2000; Neronov et al., 2005; Aharonian & Neronov, 2005; Neronov & Aharonian, 2007; Neronov et al., 2009) in which an almost monochromatic distribution of particles is maintained at the vacuum gap in the magnetosphere. The spectrum of synchrotron and/or curvature emission from electrons and/or protons accelerated in the vacuum gap can be very hard. For example, in the regime when the acceleration of particles of mass and charge in a magnetosphere with magnetic field is balanced by the synchrotron loss rate ( is acceleration efficiency), typical particle energies are . The synchrotron radiation from such particles is sharply peaked at an energy independently of the magnetic field strength. The spectrum of emission below this peak energy is hard (see e.g. Vincent & LeBohec (2010) for an example of the spectrum calculation). In principle, this hardness could be close to the limiting case for the synchrotron emission spectrum from a monochromatic electron distribution, which is characterized by a photon index . If the energy loss mechanism limiting the maximal proton energies is curvature radiation, the maximal photon energy depends on the magnetic field and on the particle energy as , where is the black hole mass which determines the size of the acceleration region. The -ray emission from the vacuum gap is highly anisotropic (see e.g. (Neronov & Aharonian, 2007)). Most of the hard spectrum luminosity is directed along the AGN jet so that the -ray flux from the vacuum gap can be comparable or higher than the -ray emission from electrons accelerated in the jet. Since the energy deposited in the accelerated particles is released in the form of the -rays (Neronov & Aharonian, 2007), the -ray production efficiency in the vacuum gap is extremely high.
3.1.5 IC from electrons within a cold relativistic flow origin
Alternatively, a narrow particle energy distribution could be maintained in a “cold” relativistic wind with very large bulk Lorenz factor, similar to the case of pulsar winds, if such winds are ejected by the central supermassive black hole. Hard spectrum -ray emission could then originate through IC emission, produced by the wind propagating through an external radiation field, as suggested by Aharonian et al. (2002).
3.2 Emission from -ray induced cascades in the intergalactic medium
Very hard -ray emission spectra can also form through the propagation of multi-TeV -rays from their source to Earth. The absorption of VHE -rays through their interactions with the EBL initiates the development of electromagnetic cascades. If the energy of the primary -rays is high enough, most of the primary source power is transferred into an electromagnetic cascade. The energy in the electromagnetic cascade is released in the form of detectable secondary cascade -ray emission with energies below the energy at which the mean free path of -rays is comparable to the distance to the source (Aharonian et al., 1994; Neronov&Semikoz, 2007). This results in the possibile formation of a hard -ray emission spectrum, resulting from the development of a cascade in the intergalactic medium, as was suggested by Aharonian et al. (2002). An example calculation of the cascade spectrum in the case of Mrk 501 is shown in Fig. 3. In this calculation primary -rays with energy 100 TeV were injected at the source. The broad band type spectrum in the energy range 0.1-10 GeV is formed during the propagation of the -ray beam towards the Earth. The observed high-energy cut-off at TeV sits at the energy for which the photon mean free path is about the distance to the source.
The low energy suppression of the cascade spectrum arises due to the influence of intergalactic magnetic fields on the geometrical structure of the cascade. If the intervening intergalactic magnetic field (IGMF) strength were negligible, secondary pairs, as well as cascade -rays would propagate in the same direction as the primary -ray beam. Small deflections of the trajectories of electrons and positrons by weak, but non-negligible magnetic fields lead to the time delay of the cascade signal (Plaga, 1995; Murase et al., 2008; Neronov & Semikoz, 2009). At energies for which the time delay is much longer than the flare duration, the cascade emission is suppressed. This time delay increases with the decrease of the energy of the cascade photons (Neronov & Semikoz, 2009) The time delay suppression of the cascade flux becomes stronger at lower energies. This is illustrated in Fig. 3 where the cascade spectrum calculated for magnetic field values of G, G and G are shown. The low-energy suppression of the spectrum is observed below energies of GeV. Monte Carlo calculations which take into account the production spectrum of the pairs and of the secondary cascade photons show that at energies for which the time delay suppression is important, the spectrum of the cascade emission is hard, with a photon index , close to the observed photon index of emission above 10 GeV in the Mrk 501 flare. Details of the Monte-Carlo code used for the calculation can be found in Taylor et al. (2011).
The range of IGMF strengths required in the cascade model is consistent with the lower bounds on the IGMF which have been recently derived from the non-observation of the -ray emission from electromagnetic cascades emission from several hard spectrum TeV blazars Neronov& Vovk (2011); Taylor et al. (2011); Tavecchio et al. (2010); Dermer et al. (2011). The exact level of this lower bound depends on the nature of suppression of the cascade emission signal (time delay or extended nature of the cascade signal, see Taylor et al. (2011) for details) and is in the range of G if the IGMF correlation length is Mpc. For shorter correlation lengths, the bound strengthens to G. This cascade explanation of the observed hard spectrum of Mrk 501 is valid only if the real IGMF strength is close to the previously derived lower bounds.
The cascade model provides an alternative possible explanation of the extremely hard spectrum in the 10-200 GeV band. The interpretation of the data in terms of this model provides a measurement of the magnetic field in the intergalactic medium along the line of sight toward Mrk 501. The duration of the flare was month. If the magnetic field were weaker than G, the characteristic time delay in the energy range GeV would be much shorter than the flare duration and no suppression of the -ray flux would be observed (dotted curve in Fig. 3). On the other hand, if the magnetic field were stronger than G, the type spectrum formed by the time delay suppression effect would continue to energies much higher than GeV. This would place heavy demands on the source TeV luminosity, requiring a TeV flare luminosity of more than an order of magnitude higher than that observed in the 100 GeV energy band. Although this possibility can not be fully ruled out (no systematic monitoring of the source during the flare was done by the ground-based -ray telescopes), such a scenario would need a source flux larger than erg/(cm s), in excess of the highest historical flux detected from the source during the large 1997 flare (the spectrum of the 1997 flare measured by HEGRA telescope (Aharonian et al., 1996) is shown in grey in Fig. 3).
An alternative possibility for the cascade contribution to the Mrk 501 spectrum was considered by Takahashi et al. (2012). Assuming that the observed TeV band flux is intrinsic to the source, rather than resulting from the development of an electromagnetic cascade, the cascade emission is still expected to appear in the GeV energy band. Takahashi et al. (2012) made an attempt to constrain the strength of IGMF along the line of sight towards Mrk 501 by searching for time delayed GeV band “pair echo” emission following the hard flare discussed in this paper. We note, however, that even if the peak flare flux in the TeV band was at the level measured by VERITAS, the Fermi data are of insufficient quality to provide sensible constraints on the IGMF. The reason for this is the high level of the steady-state flux from the source in the GeV energy band. From Fig. 1 one can see that no variations in the flux level are detected in the 0.3-30 GeV band, neither during the flare, nor in the period following the flare. The source flux in the GeV band is constantly at the level of erg/cms, which is above that expected for the pair echo flux at any moment of time (see Fig. 1 of Takahashi et al. (2012)). Thus, the pair echo signal is not expected to be detectable, no matter how weak the magnetic field is. At the same time, if the peak flux of the flare was at the level indicated by the Fermi measurement in the GeV energy band (see Fig. 3), the strength of the pair echo signal considered by Takahashi et al. (2012) would be larger by a factor 3 – 5, so that the cascade emission in the GeV band would be above the level of the steady state emission observed from the source.
The observation of a very hard spectrum of -ray emission from Mrk 501 during an “orphan” -ray flare challenges theoretical models of -ray emission from blazars. Several mechanisms for the formation of the very hard spectra are able to provide explanations for the observed source behaviour. To clarify the nature of the hard flare(s) a systematic analysis of the frequency of occurrence and properties of the hard and/or “orphan” -ray flares in the TeV -ray loud BL Lacs is needed.
Acknowledgements. The work of AN and AMT is supported by the Swiss National Science Foundation grant PP00P2_123426.
- Abdo et al. (2010a) Abdo A.A. et al., 2010a, Phys.Rev.Lett, 104, 101101
- Abdo et al. (2010b) Abdo A.A., et al., 2010b, Ap.J.Supp., 188, 405
- Abdo et al. (2011) Abdo A., et al., 2011, Ap.J., 727, 129.
- Atwood et al. (2009) Atwood B., et al., 2009, Ap.J., 697, 1071
- Aharonian et al. (1996) Aharonian F.A., et al., 1999, A&A, 349, 11
- Aharonian (2000) Aharonian F.A., 2000, New Astronomy, 5, 377
- Aharonian et al. (2002) Aharonian F.A., Timokhin A.N., Plyashechnikov A.V., 2002, A&A, 384, 834
- Aharonian et al. (2006) Aharonian F.A., et al., 2006, Nature, 440, 1018
- Aharonian et al. (2007a) Aharonian F.A. et al., 2007a, A&A, 475, L9
- Aharonian et al. (2007b) Aharonian F.A., et al., 2007b, A&A, 473, L25
- Aharonian & Neronov (2005) Aharonian F.A., Neronov A., 2005, Ap.J., 619, 306
- Aharonian et al. (2008) Aharonian F.A., Khanulyan D., Costamante L., 2008, MNRAS, 387, 1206
- Aharonian et al. (1994) Aharonian F.A., Coppi P.S., Volk H.J., 1994, Ap. J., 423, L5
- Beskin et al. (1992) Beskin V. S., Istomin Y. N., Parev V. I., 1992, SvA, 36, 642
- Böttcher et al. (2008) Böttcher M., Dermer C.D., Finke J.D., 2008, Ap.J., 679, L9
- Böttcher (2005) Böttcher M., 2005, Ap.J., 621, 176
- Daniel et al. (2005) Daniel M.K., et al., 2005, Ap.J., 621, 181
- Dermer et al. (2011) Dermer C., et al., 2011, Ap.J., 733, L21.
- Finke et al. (2008) Finke J.D., Dermer C.D., Bottcher M., Ap.J., 2008, 686, 181.
- Katarzynski et al. (2006a) Katarzynski K., Ghisellini G., Tavecchio F., Gracia J., Maraschi L., 2006, MNRAS, 368, L52
- Katarzynski et al. (2006b) Katarzynski K., Ghisellini G., Mastichiadis A., Tavecchio F., and Maraschi L., 2006, A&A 453, 47
- Krawczynski et al. (2004) Krawczynski H., et al., 2004, Ap.J., 601, 151
- Kusunose & Takahara (2006) Kusunose M., Takahara F., 2006, Ap.J., 651, 113
- Levinson (2000) Levinson A., 2000, Phys. Rev. Lett., 85, 912
- Mattox et al. (1996) Mattox J.R. et al. 1996, Ap.J., 461, 396
- Mücke & Protheore (2001) Mücke A., Protheroe R.J., 2001, Astropart.Phys., 15, 121
- Murase et al. (2008) Murase K., Takahashi K., Inoue S., Ichiki K., Nagataki S., 2008, Ap.J., 686, L67
- Neronov et al. (2005) Neronov A., Tinyakov P., Tkachev I., 2005, JETP Lett., 100, 656
- Neronov&Semikoz (2007) Neronov A., Semikoz D.V., 2007, JETP Lett., 85, 473
- Neronov & Aharonian (2007) Neronov A., Aharonian F.A., 2007, Ap.J., 671, 85
- Neronov et al. (2009) Neronov A., Semikoz D.V., Tkachev I., 2009, New J.Phys., 11, 065015
- Neronov et al. (2008) Neronov A., Semikoz D., Sibiryakov S., 2008, MNRAS, 91, 949
- Neronov & Semikoz (2009) Neronov A., Semikoz D.V., 2009, Phys.Rev. D 80, 123012
- Neronov& Vovk (2011) Neronov A., Vovk. Ie., 2010, Science, 328, 73.
- O’Sullivan et al. (2009) O’Sullivan S., Reville B., Taylor A.M., 2009, MNRAS, 400, 248.
- Plaga (1995) Plaga R., 1995, Nature, 374, 430
- Rachen&Meszaros (1998) Rachen J.P., Meszaros P., 1998, Phys.Rev.D, 58, 123005
- Reimer et al. (2004) Reimer A., Protheroe R.J., Donea A.C., 2004, A&A, 419, 89.
- Rieger et al. (2007) Rieger F.M., Bosch-Ramon V., Duffy P., 2007, Ap.SS, 309, 119.
- Takahashi et al. (2012) Takahashi K.; Mori M.; Ichiki K., Inoue S., 2012, Ap.J., 744, L7.
- Tavecchio et al. (2010) Tavecchio, F., et al., 2010, MNRAS, 406, L70.
- Taylor et al. (2011) Taylor A.M., Vovk Ie., Neronov, A., 2011, A&A, 529, A144.
- Vincent & LeBohec (2010) Vincent S., LeBohec S., 2010, MNRAS, 409, 1183
- Zacharopoulou et al. (2011) Zacharopoulou O., Khangulyan D., Aharonian F.A., Costamante L., 2011, Ap.J., 738, 157.