# Electromagnetic cascade masquerade:

a way to mimic -axion-like particle mixing effects in blazar spectra

###### Key Words.:

astroparticle physics – radiation mechanisms: non-thermal – methods: numerical – BL Lacertae objects: individual – cosmic background radiation###### Abstract

Context:Most of the studies on extragalactic -ray propagation performed up to now only accounted for primary gamma-ray absorption and adiabatic losses (“absorption-only model”). However, there is growing evidence that this model is oversimplified and must be modified in some way. In particular, it was found that the intensity extrapolated from the optically-thin energy range of some blazar spectra is insufficient to explain the optically-thick part of these spectra. This effect was interpreted as an indication for -axion-like particle (ALP) oscillation. On the other hand, there are many hints that a secondary component from electromagnetic cascades initiated by primary -rays or nuclei may be observed in the spectra of some blazars.

Aims:We study the impact of electromagnetic cascades from primary -rays or protons on the physical interpretation of blazar spectra obtained with imaging Cherenkov telescopes.

Methods:We use the publicly-available code ELMAG to compute observable spectra of electromagnetic cascades from primary -rays. For the case of primary proton, we develop a simple, fast and reasonably accurate hybrid method to calculate the observable spectrum. We perform the fitting of the observed spectral energy distributions (SEDs) with various physical models: the absorption-only model, the “electromagnetic cascade model” (for the case of primary -rays), and several versions of the hadronic cascade model (for the case of primary proton). We distinguish the following species of hadronic cascade models: 1) “basic hadronic model”, where it is assumed that the proton beam travels undisturbed by extragalactic magnetic field and that all observable -rays are produced by primary protons through photohadronic processes with subsequent development of electromagnetic cascades 2) “intermediate hadronic model”, the same as the basic hadronic model, but the primary beam is terminated at some redshift 3) “modified hadronic model” that includes the contribution from primary (redshifted and partially absorbed) -rays.

Results:Electromagnetic cascades show at least two very distinct regimes labeled by the energy of the primary -ray (): the one-generation regime for the case of ¡10 and the universal regime for ¿100 and redshift to the source ¿0.02. Spectral signatures of the observable spectrum for the case of the basic hadronic model, = 0.186 and low energy (¡200 ) are nearly the same as for purely electromagnetic cascade, but for ¿200 the spectrum is much harder for the case of the basic hadronic model. In the framework of the intermediate hadronic model, the observable spectrum depends only slightly on the primary proton energy, but it strongly depends on at ¿500 . As a rule, both electromagnetic and hadronic cascade models provide acceptable fits to the observed SEDs. We show that the best-fit model intensity in the multi- region of the spectrum in the framework of the electromagnetic cascade model is typically greater than the one for the case of the absorption-only model. Finally, for the case of blazar 1ES 0229+200 we provide strong constraints on the intermediate hadronic model assuming the blazar emission model of Tavecchio (2013) and the model of magnetic field around the source according to Meyer et al. (2013).

Conclusions:

## 1 Introduction

The development of ground-based -ray astronomy with imaging Cherenkov telescopes has been very fast during the last two decades (for a review see, e.g., Hillas (2012)). Indeed, only 7 years after the detection of photons from the active galactic nucleus (AGN) Mkn 421 in Punch et al. (1992), the first detailed study of a blazar (AGN with the jet presumably pointed towards the observer) spectrum was made by Aharonian et al. (1999). Almost immediately, these observations were utilized to put some constraints on the intensity of extragalactic background light (EBL) (in Stecker & de Jager (1993), de Jager et al. (1994) for the former observations, and in Aharonian et al. (1999) itself for the latter).

Indeed, primary very high energy (VHE, ¿100 ) -rays are absorbed on EBL photons by means of the process (Nikishov (1962), Gould & Shroeder (1967)), and for the case of primary energy and higher — on the cosmic microwave background (CMB) as well (Jelley (1966)). More recently, the signatures of this absorption process were observed with the Fermi LAT instrument (Ackermann et al. (2012)) and the H.E.S.S. Cherenkov telescope (Abramowski et al. (The H.E.S.S. Collaboration) (2013)) with high statistical significance ( and 8.8 , respectively).

However, Horns & Meyer (2012) found that the strength of absorption at high optical depth (¿2, hereafter simply ) appears to be lower than expected. This result was obtained on a sample of blazar spectra measured with imaging Cherenkov telescopes by comparing the distributions of the flux points for the regions 1¡¡2 and ¿2 around the intensity extrapolated from the optically-thin regime ¡1. The statistical significance of this effect, according to Horns & Meyer (2012), is 4.2 . While this result was not confirmed by Biteau & Williams (2015), very recently Horns (2016) again found an indication for this anomaly with another analysis method. Such an anomaly, in fact, closely resembles the so-called “-IR crisis” (Protheroe & Meyer (2000)) that was derived from the already mentioned observations of Mkn 501 (Aharonian et al. (1999)). The “-IR crisis”, however, was later found to be less severe after the development of more advanced EBL models. On the contrary, the “new” anomaly of Horns & Meyer (2012), Horns (2016) persists for most of these contemporary models of EBL intensity.

The authors of Horns & Meyer (2012) interpreted their result as an indication for the existence of some non-conventional physical effect, for instance, the process of oscillations of -rays into axion-like particles (ALPs) and back into photons in magnetic field on the way from the source to the observer (). Indeed, a part of photons reconverted from ALPs near the observer can significantly enhance the observed intensity in the ¿2 region. Moreover, once the anomaly is well established, it is possible to put constraints on the gamma-ALP mixing parameters (the mass of ALP) and (the photon-ALP coupling constant). In Meyer et al. (2013) a lower limit on the was found, depending on the value. For any fixed in the range considered in Meyer et al. (2013), some values of greater than the lower limit could explain the observed anomaly. Together with the upper limits from the CAST experiment (Andriamonje et al. (CAST Collaboration) (2007)), this result allowed to construct a confidence interval for .

Besides the anomaly at ¿2, there exists another signature of oscillation, namely, a step-like irregularity that is situated at the energy lower than the starting energy of the VHE anomaly (e.g, Sanchez-Conde et al. (2009)). The drop of intensity associated with this spectral feature is usually about 1/3 as photons (two polarization states) attain equipartition with ALPs (one polarization state). A very recent analysis of Fermi LAT data (Ajello et al. (2016)) (observations of the NGC 1275 Seyfert galaxy were used), however, did not find this signature. Moreover, other bounds (Ayala et al. (2014), Abramowski et al. (H.E.S.S. Collaboration) (2013), Payez et al. (2015), Wouters & Brun (2013)) on the parameters of -ALP mixing allowed to strongly constrain the scenario considered in Meyer et al. (2013). Therefore, the hypothesis that the VHE anomaly is explained by oscillation appears to be less attractive than before; one needs to search for some other physical mechanism of the anomaly.

In this respect, we note that most of extragalactic -ray propagation studies were performed with account of only two elementary processes: the absorption of primary photons (by means of pair-production) and their adiabatic losses. This model (hereafter “the absorption-only model”) rests on the assumption that the secondary electrons and positrons (hereafter simply “electrons” unless otherwise stated) are deflected and delayed by the extragalactic magnetic field (EGMF), thus the cascade photons produced by these electrons by means of the inverse Compton (IC) scattering do not contribute to the observed spectrum.

However, the EGMF strength in voids of the Large Scale Structure (LSS) may be small enough to violate this assumption. The existing constraints (Blasi et al. (1999), Pshirkov et al. (2016), Dolag et al. (2005), Neronov & Vovk (2010), Dermer et al. (2011), Taylor et al. (2011), Vovk et al. (2012), Abramowski et al. (H.E.S.S. Collaboration) (2014), Takahashi et al. (2012), Takahashi et al. (2013)) on the EGMF strength on characteristic coherence scale 1 (Akahori & Ryu (2010)) (some of them were summarized, e.g., in Fig. 4 of Dzhatdoev (2015a)) do not exclude the probability that cascade emission may contribute to the observed spectrum.

The recent work Finke et al. (2015), using a more conservative method than the above-mentioned references, excluded , again on characteristic coherence scale 1 , with statistical significance . Arlen et al. (2014) did not find any reason to reject the null hypothesis of at all. Tashiro & Vachaspati (2015) found . This study was made using the angular correlation pattern of Fermi LAT (Atwood et al. (2009)) diffuse -rays; however, it is not clear how much this result may be affected by the existence of comparatively strong magnetic field in galaxy clusters.

Moreover, there are some hints that the cascade component does contribute to the observed spectrum of blazars at energies . Neronov et al. (2012), using observations of Mkn 501 in a flaring state by the Fermi LAT instrument and the VERITAS Cherenkov telescope (Abdo et al. (2011)) during the 2009 multiwavelength campaign, obtained the energy spectrum ( is the number of photons per energy bin ) from 300 up to 5 . It was found that the Fermi LAT spectrum in the 10–200 energy range had a power-law () index , while in the 300 – 5 . It is interesting that the Fermi LAT lightcurves in the 0.3-3 and 3-30 energy bins do not show any evidence for strong, fast variability, while in the 30-300 energy bin the flare is readily identified on the lightcurve. Neronov et al. (2012) found that such a behaviour of the spectral and timing characteristics is typical for intergalactic cascade, assuming on the maximum spatial scale 1 .

Another result supporting the incompleteness of the absorption-only model was also obtained with the Fermi LAT telescope. Namely, Furniss et al. (2015) found a correlation between 10-500 energy flux of blazars with relatively hard () observed spectra above 10 and the fraction along the line of sight occupied by voids in the LSS. This effect may be explained by the same physical mechanism: the cascade emission that is likely to be angle-broadened by the EGMF at some energy below 200-300 . Finally, Chen et al. (2015), again with Fermi LAT data, found the evidence for the existence of “pair halos” (extended emission) around the positions of various blazars, thus giving additional support to the hypothesis that the cascade component may contribute to the observed spectra.

These works motivated us to study how the inclusion of the cascade component into the fitting of the VHE spectra would influence the data interpretation, especially in the optically thick (¿1) region. Our paper is by far not the first study of intergalactic cascade spectra; besides the already mentioned work (Neronov et al. (2012)), there are many others that included more or less detailed duscussions of these, namely: Aharonian et al. (1999), Aharonian et al. (2002), d’Avesac et al. (2007), Murase et al. (2012), Takami et al. (2013).

However, very recently in a conference paper Dzhatdoev (2015b) it was shown that if the primary spectrum is hard enough () and the cutoff in the spectral energy distribution (SED) of the source is situated at a sufficiently high energy ( ), then the intersection of a low-energy cascade component and a high-energy primary (absorbed) component forms a kind of an “ankle” that to some extent may mimic the mixing effect mentioned above. One of the main aims of the present paper is to discuss this effect in more depth. We will see that the last model is qualitatively different from other “electromagnetic cascade models”, i.e. the models that involve intergalactic cascades initiated by primary gamma-rays.

There exists another class of intergalactic cascade models of blazar spectra dealing with primary protons or nuclei that initiate secondary photons by means of photopion losses and pair production with subsequent development of electromagnetic cascades (“hadronic cascade models”) (e.g. Uryson (1998), Essey & Kusenko (2010a), Essey et al. (2010b), Essey et al. (2011), Murase et al. (2012), Takami et al. (2013), Essey & Kusenko (2014), Zheng et al. (2016)). We will compare hadronic models with the electromagnetic cascade model of Dzhatdoev (2015b) to reveal their advantages and diffuculties.

The paper is organized as follows. In Sec. 2 we discuss some general properties of cascades from primary photon or proton developing on EBL/CMB and give a detailed description of our calculation method for the case of primary proton. In Section 3 we define the sample of blazar spectra to be used in the analysis (in this paper we focus on the observations made with Cherenkov telescopes). In Section 4 we present the fits to observed SEDs. Section 5 contains discussion, where our main results are recalled and discussed in the broader context; and, finally, in Section 6 the conclusions to our work are drawn.

## 2 Electromagnetic cascades from primary gamma-rays and primary nuclei

After a primary gamma-ray, travelling through the Universe, is absorbed by an EBL photon, secondary electrons may upscatter CMB/EBL photons and produce the new generation of “cascade” photons by means of the IC process. If the energy of the primary photon is high enough and the source is sufficiently distant, the number of generations in such a cascade may be considerably higher than unity. In this case the spectrum of the cascade takes a universal form. On the contrary, if the energy of the primary photon is sufficiently low (but still high enough to be absorbed on EBL), the cascade takes a “degenerate” form and may be well described by a one-generation approximation (see. e.g., Dermer et al. (2011)). In the next subsection we discuss both regimes, as well as the transition between them.

When the calculations presented in this section were almost finished, we noticed an insightful paper Berezinsky & Kalashev (2016) that discusses the universal regime in depth. We will express our results in terms of this work, when applicable. All results presented in this section are calculated for the case of a cascade particle threshold of 10 , unless otherwise stated, and the redshift to the source = 0.186. We find it convenient as we will see that 3 out of 6 blazars studied in this paper (see Section 3) have redshifts very close to 0.186. In the present paper we use the ROOT (Brun & Rademakers (1997)) analysis framework for drawing figures and performing a part of analysis.

### 2.1 Cascades from primary gamma-rays

Throughout this paper we use the publicly-available code ELMAG (Kachelriess et al., (2012)) version 2.02 to calculate the observable spectra of -rays for the case of primary photon. In what follows, when running the ELMAG code, we assume the Kneiske & Dole (2010) (hereafter KD10) EBL model. The KD10 model was the first to be implemented in the ELMAG code, therefore this code is best tested (and most reliable) with this model. All calculations with ELMAG presented in this paper are full direct MC simulations (with parameter ).

There was a considerable amount of discussion whether pair beams resulting from the development of electromagnetic cascades are subject to plasma instabilities and additional (with respect to the IC process) energy losses. The conclusions of various studies are conflicting: while Broderick et al. (2012), Schlickeiser et al. (2012), Chang et al. (2014), Menzler & Schlickeiser (2015) found that these effects are considerable and must be accounted for, at least in the Fermi LAT region of the spectrum, Miniati & Elyiv (2013), Venters & Pavlidou (2013), Sironi & Giannios (2014), Kempf et al. (2016)) argue that plasma losses are subdominant with respect to IC losses. The last work from this list (Kempf et al. (2016))), running a detailed numerical simulation, found that the growth rate of the instability is likely not sufficient to cause an appreciable effect. Therefore, we do not include such a process in our calculations.

To demonstrate the basic regimes of electromagnetic cascade development on CMB/EBL, Fig. 1 presents observable spectra for different primary energies for the case of monoenergetic primary injection. Two principal components may be readily identified in Fig. 1 — the cascade component itself and the primary (absorbed) component that is seen in the right part of the graph for the case of 10 . The higher the energy, the more the primary component is attenuated due to a rapid rise of the value with increasing energy. A slight shift of the primary energy to the lower values is due to adiabatic losses. For the case of the primary energy above 10 , the primary photons are almost completely absorbed due to very large value of for these energies.

Fig. 1 reveals two very different regimes: for the case of 10 the energy of the maximum in the spectral energy distribution (SED=) of the cascade component is approximately , while, on the contrary, for the case of 100 the cascade spectrum is practically independent of energy. The first case may be well described by the one-generation approximation; the second one is essentially the universal regime. As we will see, the universality of the cascade spectrum is a “weak” one (Berezinsky & Kalashev (2016)), i.e. the spectrum of observable -rays is practically independent of the energy and type (photon/electron) of the primary particle, but depends on . The transition between these regimes appears to be very sharp; the cascade completely changes its mode of propagation on only one decade of the primary energy. The physical reason for such a fast transition is clear: the rapid rise of the interaction rate in the energy range of 10–100 (see, e.g., Kachelriess et al., (2012), Fig. 2, left panel). Up to our knowledge, no study before us has remarked such a fast transition between these completely different regimes of cascade development.

Additional graphs of electromagnetic cascade spectra in the universal regime are presented in Appendix A. They demonstrate that the property of cascade universality is fulfilled in the observable energy range 10 – 30 within 4 decades of the primary energy (from 100 to 1 ) with an accuracy 30 % independently of the primary particle type for 0.02.

To be able to distinquish between two categories of observable -rays: primary (redshifted) (type 1) and cascade (type 2) photons, we modified the output of the ELMAG code and created a database that contains observable spectrum for every primary photon. For every array that contains observable spectrum we performed the following classification procedure. If this array contains only one non-zero entry and the redshifted value of the primary energy fits the energy range of the corresponding bin with the non-zero enrty, this array was classified as a type 1 event, otherwise — as a type 2 event. Appendix B contains several graphs that show primary, absorbed and cascade components for various primary spectra. These figures demonstrate that the contribution of the cascade component to the total intensity is appreciable at 100 (roughly the threshold for Cherenkov telescope observations) and above only if the primary spectrum is hard enough so that and that the cutoff energy of the spectrum 10 .

### 2.2 Cascades from primary proton

For the case of primary proton we developed an original code based on a hybrid approach. First of all, we propagate primary protons from the source to the observer (= 0) with a small step , updating their energy at every step, and compute their mean energy loss according to Berezinsky et al. (2006), equation (8). In this work we are interested in the case of proton energy losses on CMB. Full energy losses of a proton with energy per unit time are:

(1) | |||

(2) | |||

(3) |

where = 67.8 , = 0.308 (Ade et al. (The Planck Collaboration) (2015)), , and — pair production and pion production energy losses at .

The next step is the calculation of energy spectra of particles produced by primary protons. We follow Kelner & Aharonian (2008) to calculate these energy spectra. For the pair production process the spectra of electrons () and positrons are identical and for the case of pair production on CMB (Kelner & Aharonian (2008), equation (67)):

(4) |

where is the Boltzmann constant; (= 2.725 is the CMB temperature at =0); is proton Lorentz factor; , where is the energy of secondary electron in the comoving frame (i.e. in the same frame where is measured); is the energy of the photon in the rest frame of the proton in units of electron rest energy (, are four-momenta of the proton and photon, respectively); , where and are the energy and the momentum of electron in the rest frame of the proton, respectively (); is the double-differential cross section as a function of energy and emission angle of the electron in the rest frame of the proton (Blumenthal (1970), equation (10)) (, where the angle between the momenta of the photon and the electron is denoted as ). Four examples of calculated electron SED for the case of are shown in Fig. 2 for different values of primary proton energy . The maximum of the SED for the case of = 100 was normalized to 1 in this figure.

The case of photopion losses was also covered in Kelner & Aharonian (2008); we use equations from Section II of their work with energy of CMB photon . As we are interested only in primary energies below 100 , we neglect the flux of and (which, according to Kelner & Aharonian (2008), Fig. 6 (right panel), contributes only about 0.1 % of the total flux at = 100 and =0) and calculate the spectra of (,,,,). We have checked that both panels of Kelner & Aharonian (2008) (Fig. 6) are well reproduced for all types of particles for which we perform our calculations. As an example of this calculation we present Fig. 3 where the SEDs of secondary gamma-rays and positrons are presented for the case of = 0 for different values of primary proton energy. The maximum of the SED for the case of -rays and = 100 was normalized to 1 in this figure.

Now we proceed to calculate observable (cascade) -ray spectra. The energies of secondary -rays and electrons in Fig. 2–3 are mainly above 100 , therefore the cascade universality assumption may be applied to reduce computer time requirements; such a calculation is presented in the following Subsubsection 2.2.1. As an alternative, to validate the results obtained in the cascade universality approximation, we performed a full calculation of observable spectra without this assumption (Subsubsection 2.2.2). Synchrotron radiation was neglected in these calculations. In both these Subsubsections we are interested only in the shape of the observable -ray spectrum. The normalization of the spectrum will be discussed in Section 5.

#### 2.2.1 Observable -ray spectra in the universal spectrum approximation

Here we assume “weak universality”, as described above. Using the ELMAG code, we calculated an array of cascade spectra for the case of primary -rays with fixed energy 1 but different , distributed randomly and uniformly from 0 to 0.30. While adiabatic losses do affect the energy of propagating protons, only “active” (pair and photopion) losses give rise to new particles and, therefore, eventually produce observable signal in -rays. The energy transferred to these particles on the step is:

(5) | |||

(6) |

Fig. 2–3 demonstrate that nearly all secondary particles are ultrarelativistic, therefore, nearly all the energy of these particles (except neutrinos) could be transferred to cascade. Several examples of (i.e. normalized to the value at = 0) are shown in Fig. 4 for different values of primary proton energy. For = 50 and especially 100 primary protons quickly lose energy near the source until they reach comparatively low energy so that the energy loss rate becomes much smaller (see Fig. 1 of Berezinsky et al. (2006)). On the other hand, for the case of = 10 and 30 protons keep an appreciable fraction of their energy until they reach , and, as a result, they keep producing secondaries quite effectively even near the observer.

Finally, we calculate the observable spectrum of -rays at =0. For the case of an isotropic source of primary protons:

(7) |

where is the fraction of “active” losses transferred to -rays and electrons, and is the universal spectrum of cascade with starting point at . Several examples of are shown in Fig. 5 for different values of primary proton energy. For the case of 30 , when the dominating source of proton energy loss is the pair-production process, is practically equal to 1.

In this paper we study the case of ultrarelativistic primary protons, 1 , therefore, in the absence of EGMF, the angular distribution of observable -rays is practically coincident with the angular distribution of primary protons. Thus, equation (7) is justified also for the conditions investigated in our work.

The resulting observable SED is shown in Fig. 6 for the case of monoenergetic primary injection with = 30 for two options of cascade -ray threshold energy in ELMAG. The signatures of the spectrum are highlighted with additional power-law approximations at various energy regions. The low-energy signatures of the observable spectrum (i.e. the parameters of a “poly-gonato” power-law approximation at ¡200 ) are practically identical to the case of electromagnetic cascade from primary -ray (see Berezinsky & Kalashev (2016) for a detailed discussion of the latter case). Indeed, for the case of our calculations (i.e. primary proton) these parameters are: the shape of the black power-law segment is = 1.60 (for spectrum) below the break at = 230 , and = 1.85 below another break at = 200 (red segment). The corresponding typical values for electromagnetic cascade from primary -ray are: 1.5, 100 , 1.9, and 200–400 depending on the fitting conventions (see Fig. 3, upper panel, and Fig. 8 of Berezinsky & Kalashev (2016)).

However, for the higher energy (¿200 ) the situation is very different: for the case of primary proton the observable spectrum is far less steep than for the case of primary -ray, = 2.40 between 200 and = 15.2 (green segment), and = 4.00 above 15.2 (blue segment). In this work we are interested only in below 100 ; we leave the 100 energy region of -ray observable spectrum for further stidues.

The values of the power-law indices and the break energies were derived in Berezinsky & Kalashev (2016) as follows (below we use notations slightly different from this paper). A simplified “dichromatic” model of the target photon field (CMB+EBL) was devised, assuming that all CMB photons have energy = , and all EBL photons — = . Then, a quantity was introduced meaning the minimum energy of a primary -ray that undergoes effective absorption on EBL photons; the corresponding energy of electrons of the produced pair are . Secondary -rays produced by electrons with energy typically have energy (Blumenthal & Gould (1970)); putting together the expressions for , , and , .

Denoting the number of cascade electrons that has energy during the whole time of cascade development as and taking into account that the number of electrons in each cascade generation and the energy of these electrons has nearly the same distribution as the energy of the -rays of the same generation, in average for the energy range (as the primary energy ). Most electrons are produced above , therefore in the low-energy range . Finally, the spectrum of cascade photons may be found from the following equation: ; again using the expression for secondary -ray energy , for (corresponding to ) we obtain for the shape of the observable spectrum , which is similar to . In a like manner, for the shape of the observable spectrum is ; numerical simulations performed in Berezinsky & Kalashev (2016) indicate in the same energy region. The simplified analytic model considered above assumes that all primary -rays with energy are completely absorbed and, consequently, the spectrum has an abrupt cutoff at this energy. By detailed numerical calculations presented here we reproduce the more realistic shape of the high-energy cutoff, both for primary -rays and primary protons.

Fig. 7 shows observable SEDs for different values of . For = 10, 30, 50 the shape of the spectra are nearly identical, while for the case of = 100 the spectrum is somewhat steeper due to a pile-up of the dependence near the source of primary protons. Indeed, in the latter case more energy is injected near the source, therefore, more electromagnetic cascades start to develop at comparatively high distance from the observer, thus leading to lower intensity at high values of the observable energy.

As well, different components of the observable SEDs are shown — the one formed by -rays and electrons that were produced by primary protons at 0.06, as well as the other contribution for 0.06. The first component dominates at 5 and has almost the same shape of the spectrum for all considered. This effect may be qualitatively understood as follows. The spectrum of the first component is defined by an equation similar to equation (7), but with limits of integration on redshift equal to 0 and 0.06. Primary protons with energy 50 quickly lose energy to pion production (see Berezinsky et al. (2006)), until they reach the energy range where pair production losses dominate. Therefore, at 0.06 even for =100 (see Fig. 5). The dependence is also similar for all considered at 0.06 (see Fig. 4). Again, this is due to the fact that for 50 protons quickly lose energy so that their energy at 0.06 is similar irrespectively of the value. Thus, the shape of the observable SEDs at 20 , where the first component is strongly dominant, is almost the same. In the 200 –10 energy region and 100 the first component has lower relative normalization (with respect to the second component), once again due to the fact that in this case a larger fraction of primary energy is lost at 0.06 than for 50 . The second component for 100 is also somewhat steeper for 100 due to the same effect. Therefore, the overall observable SED in the 200 –10 energy region is steeper for 100 than for the other considered values of . The difference in the slope between the observable SEDs in this energy range may be estimated as follows. The total normalization of the second component is , and of the first component — ; is the ratio of these quantities; the normalization factor is chosen so that at 30 . Direct numerical calculation yields the values at 10 , at 30 , at 50 , and corresponding to 100 ; this translates to the difference in the power-law slope between the case of 100 and other cases.

#### 2.2.2 Observable -ray spectra in the full hybrid approach

To obtain the observable spectrum without the assumption of cascade universality, we calculated an array of cascade spectra from primary -rays and electrons of energies 10 –100 . From 10 to 10 the primary spectrum of these cascades () had a power-law slope 1, and the slope from 10 to 100 was set to 2. As in the universal spectrum approach, we propagated primary protons, but now, instead of the universal spectrum, we utilized our new database of cascade spectra. We calculated a vast two-dimensional array of secondary particle spectra for the case of the pair production process, as described above, as well as for the photopion process for (,,,,), both for the case of 101 logarithmically spaced values of from 1 to 100 , and for 30 linearly spaced values of from 0 to 0.30. These arrays were used to calculate the observable gamma-ray spectrum, that is:

(8) | |||

(9) | |||

(10) |

where and are, as before, the spectra of cascades, but now the shape of these spectra may depend not only on , but also on the primary energy and the type of the primary particle. and are the spectra of secondary electrons and -rays, respectively, that were produced by primary protons.

The results of our calculations of observable SEDs are presented in Figs. 8–11. The comparison graph of spectra calculated with two different approaches for the case of = 30 (black solid line and circles in Fig. 8) shows very good agreement. We also considered a possibility of a highly structured extragalactic magnetic field (for instance, a galaxy cluster) that is situated at some redshift between the source and the observer. We assume that magnetic field strength in this cluster is sufficiently high so that the proton beam traversing the cluster is practically dissolved. Several cases of observable spectra for different are also shown in Fig. 8. In all these cases the agreement of results obtained with two different approaches is good. Finally, Fig. 8 shows that in the case when the observable spectrum is very near to a universal cascade spectrum from primary -rays. The same SEDs, but without relative normalization, are shown in Fig. 9.

N | Source | Observational period | Reference | |
---|---|---|---|---|

1 | H 1426+428 | 0.129 | 1999-2000 | Aharonian et al. (2003) |

2 | H 1426+428 | 0.129 | 1998-2000 | Djannati-Atai et al. (2002) |

3 | H 1426+428 | 0.129 | 2001 | Horan et al. (2002) |

4 | 1ES 0229+200 | 0.140 | 2005-2006 | Aharonian et al. (2007a) |

5 | 1ES 0229+200 | 0.140 | 2010-2012 | Aliu et al. (2014) |

6 | 1ES 1218+304 | 0.182 | 2012-2013 | Madhavan et al. (2013) |

7 | 1ES 1101-232 | 0.186 | 2004-2005 | Aharonian et al. (2007b) |

8 | 1ES 1101-232 | 0.186 | 2004-2005 | Aharonian et al. (2006) |

9 | 1ES 0347-121 | 0.188 | Aug.-Dec. 2006 | Aharonian et al. (2007c) |

10 | 1ES 0414+009 | 0.287 | 2005-2009 | Abramowski A. et al. (2012) |

The same results as in Fig. 8, but for = 100 , are shown in Fig. 10. For 0.05 the agreement between the results obtained with our two methods is rather good, while for 0.10 SEDs resulting from the universal hybrid approach are more steep at comparatively high values of observable energy, ¿500 . This property of observable SEDs is due to a possible hardening of electromagnetic cascade spectra shape for the case of high primary energy, 1 , see Fig. 26. This hardening is the result of the “extreme high-energy cascade regime” that is realized when , where is the energy of EBL/CMB photon, and is the mass of electron. In this regime, for the case of the pair production process, one electron of every produced pair gets almost all of the energy of the primary photon, and, likewise, the secondary photon in the IC process gets almost all of the energy of the primary electron (for comparatively recent discussions see Khangulyan et al. (2008), Aharonian et al. (2012)). As a result, electromagnetic cascade develops much more slowly than in the universal regime.

We note, however, that the ELMAG code does not account neither for a possible impact of the so-called universal radio background (URB) to the cascade development (Protheroe & Biermann (1996), Settimo & De Domenico (2015)), nor for the higher-order quantum electrodynamics (QED) processes such as double pair production (e.g. Brown et al. (1973), Demidov & Kalashev (2009)) or “triplet production” (Mastichiadis et al., (1994)). The results presented in Settimo & De Domenico (2015) (Fig. 2, bottom) show that the mean free path for the IC process is likely to be affected by these effects even for the primary electron energy of 10 , thus making the “extreme high-energy cascade regime” less pronounced. In what follows we mainly use results obtained in the universal cascade regime. Finally, the same SEDs as in Fig. 10, but without relative normalization, are shown in Fig. 11.

## 3 The sample of blazar spectra

In this study we focus on a subsample of AGN, the so-called extreme blazars (Bonnoli et al. (2015)). Blazars are defined as a class of AGN that has peculiar multi-wavelength properties from radio to -ray bands of the spectrum. From the point of view of a gamma-ray astronomer, a blazar may be defined simply as an active galactic nucleus that is a bright source of -ray emission. Extreme blazars are believed to have particularly high peak energy in their SED, 1 (Bonnoli et al. (2015)).

During the last decade, these sources became very important for the studies of absorption on EBL (Aharonian et al. (H.E.S.S. Collaboration) (2006)). Indeed, these sources allow to study this effect for the highest values of available now. Additionally, they may have a very hard primary spectrum of -rays, thus making the contribution of the cascade component significant even in the VHE energy range.

The sample of observations of extreme blazars used in this paper is presented in Table 1. It contains 9 independent observations of 6 sources, performed with imaging Cherenkov telescopes HEGRA, CAT, Whipple, H.E.S.S., and VERITAS. For the whole history of observations of these sources in the VHE energy region up to 2016-01-01 (Wakely & Horan (2016)), 5 out of 6 of them display a rather slow (if any) variability (with characteristic period of months or even years). The only exception is 1ES 1812+304; for this source a rapid flare with full width on half magnitude about 2-3 days was observed in 2008 by the VERITAS telescope (2010).

## 4 Fitting the spectra of extreme blazars

In this section we compare the predictions of different models with observations. Before we proceed to present the fits for various sources, some remarks are in order. The collaborations operating Cherenkov telescopes usually reconstruct the spectrum of primary -rays, and not just the histogram of measured energy values. Therefore, while presenting our results, we will not perform any convolution of model spectra with instrumental energy resolution templates. Of course, the procedure of the primary spectrum deconvolution always results in some additional systematic uncertainty of the reconstructed spectrum. This systematics, however, for the case of observations listed in Table 1 is usually subdominant with respect to statistical uncertainties, especially in the optically thick region, to which we pay most attention in this work.

Another remark concerns the selection of EBL model for the case of absorption-only fits. Meyer et al. (2012) calculated the significance of the VHE anomaly for different EBL models (Franceschini et al. (2008) (hereafter F08), Kneiske & Dole (2010) (already denoted as KD10), and Dominguez et al. (2011) (D11)) and found that may depend non-trivially on the total normalization of the EBL intensity in a certain wavelength band. The lowest was obtained for the case of the F08 model. Therefore, in what follows we present absorption-only fits using the F08 model; for comparison, we also include the fits for G12 model, as well as for the model directly derived from the ELMAG code (see next subsection for more details).

### 4.1 The case of 1ES 0347-121

Let us start our discussion with the case of blazar 1ES 0347-121 and the absorption-only model. The redshift of this source = 0.188 is very near to 0.186, for which all calculations presented in Section 2 are applicable. The shape of the primary spectrum was chosen as:

(11) |

where is the cutoff energy. For this model we neglect adiabatic losses, because they do not change the shape of the spectrum.

The ROOT analysis framework with the integrated minimization system MINUIT (James & Roos (1975)) (namely, the gradient optimization routine MIGRAD) was used to obtain this fit. During the optimization procedure, every experimental bin was divided to 20 small parts; the flux extinction factor was calculated for each of these sub-bins in order to ensure realistic implementation of the -ray absorption process. After that, a histogram of model SED was evaluated, and a minimization was performed with the MIGRAD routine. The output of the routine includes the optimized values of primary spectrum parameters.

The result of fitting for the absorption-only model is shown in Fig. 12, top-left panel. We note that this and the following spectra are presented in the form of histograms with dense energy sampling using narrow model bins rather than somewhat sparse experimental histograms that were in fact computed during the optimization procedure. Together with fits for the case of G12 and F08 EBL models, we present a similar result for the case of EBL model as implemented in ELMAG. To do this, we evaluated optical depth vs. energy for this model as:

(12) |

where is the number of bin in the energy histogram; is the histogram of primary photons, but redshifted in order to account for adiabatic losses; is the histogram of detected photons, that were not absorbed (i.e. corresponds to certain narrow interval of observable energies for both histograms). Appendix C demonstrates that for the case of KD10 model the ELMAG code slightly (about 10-20 %) overestimates in the 1-10 energy range. However, as the actual EBL level is still uncertain, this EBL model that differs from the original KD10 model is still perfectly acceptable for our study.

Top-right panel of Fig. 12 presents a fit in the framework of the electromagnetic cascade model of Dzhatdoev (2015a), Dzhatdoev (2015b). To obtain this fit, we generated a large array of cascades from primary gamma-rays using the ELMAG code with power-law spectrum and primary energies from 100 to 100 . Observable spectra for any other primary spectrum of -rays may be calculated using a re-weighting procedure with a weight defined as:

(13) |

To test this spectrum synthesis procedure we present a comparison graph (see Appendix D, Fig. 30) of spectra calculated by Vovk et al. (2012) using the F08 EBL model and our results obtained with the primary spectrum parameters from Vovk et al. (2012). Notwitstanding slightly different EBL models, the agreement between the total observable model spectra is rather good.

As for the case of the absorption-only model, we generated histograms of observable spectra for a set of parameters and performed optimization over these parameters, this time evaluating on a grid of 1010 cells with exhaustive calculation. The values of that ensure the minimal value of on this grid were found, and the optimization run was repeated on a similar grid, but with a smaller step . This procedure was repeated several times until and . The step decreasing factor was set to 3. We have put considerable efforts into ensuring that the global minimum of is captured by our optimization method.

The best-fit spectrum consists of two distinct components — the primary (absorbed) that dominates at high energy, and the cascade one, that has rather steep spectrum and contributes mainly in the low-energy region of the total spectrum. These two components, when put together, produce an ankle-like spectral feature at an energy around 1 . The same fit, but with different energy and intensity scale, is shown in Fig. 12, middle-left panel. It is clearly seen that the cutoff in the primary spectrum is situated well below 100 .

Fig. 12, lower-left panel was calculated for the case of the “basic hadronic cascade model” of Essey et al. (2010b), where all observable -rays are of secondary nature, i.e. were produced by protons on the way from the source to the observer. Apart from the option without any considerable EGMF on the line-of-sight, we also include several options with a possible magnetic field structure at , as was described in Subsubsection 2.2.2. Such calculations presented in this section were performed in the weak universality approximation with = 30 . The case of power-law spectrum of primary protons from 1 to 100 with is also included to this Figure.

The case of the “modified hadronic cascade model” (Essey & Kusenko (2014)) that includes both primary and secondary components is shown in two lower panels of Fig. 12. Optimization was performed on 4 parameters: the normalization factors of the primary -ray and cascade components, and the parameters of the shape of the primary component (see equation (11)).

Exhaustive search is hardly possible on such a four-dimensional grid; therefore, we have developed a code in the ROOT framework using the Minuit package to perform gradient optimization. The formal best fit obtained with this procedure is presented in Fig. 12, bottom-left. For the case of the primary -ray component we neglected weak cascade emission produced by this component. In this case the spectrum of the primary component appears to be very narrow, with sharp lower- and higher-energy cutoffs. In Fig. 12, bottom-right we show, however, that even in the case of a more realistic primary component a good fit to the observed SED is obtainable.

Now let us introduce a quantity called the modification factor (following Berezinsky et al. (2006)) also known as the flux boost factor (Sanchez-Conde et al. (2009)) defined as the ratio between the spectra in the electromagnetic cascade model (ECM) and the absorption-only model (AOM):

(14) |

The graph of for blazar 1ES 0347-121 and KD10 EBL model, as implemented in the ELMAG code, is shown in Fig. 13. For calculations of the boost factor, and were interpolated to 1000 bins. is clearly greater than unity for ¿2 ; the ratio of the maximal to the minimal values of is about 3.5. Therefore, electromagnetic cascade model predicts that the intensity in the optically thick region may significantly exceed the one deduced from the fitting of the optically-thin part of experimental SED in the framework of the absorption-only model. Indeed, the cascade component, dominating at low energies, ensures good quality of fit at these energies notwithstanding a very hard primary spectrum. The cascade component, in effect, conceals, or “masks”, the primary component in the optically-thin energy region (hence the name of the present paper). On the other hand, at high energies where ¿2, a very hard primary spectrum provides enough photons to explain observations even after extinction caused by the EBL.

### 4.2 The case of 1ES 0229+200

Blazar 1ES 0229+200 is a frequent source in many discussions of extragalactic -ray propagation due to its hard spectrum, slow variability, and comparatively high total flux. The fits for the case of the absorption-only model, which are similar in every respect to those presented in Fig. 12 for 1ES 0347-121, are shown in top-left panel of Fig. 14.

For this source we also present two fits for a model that includes the oscillation process (top-right panel of Fig. 14), but without account of any secondary (cascade) emission. As was discussed in Introduction, extra photons with respect to the case of the absorption-only model in the optically-thick region of the spectrum may be ascribed to this process. To evaluate the effect of -ALP mixing in the spectrum of 1ES0229+200, we used the results of Sanchez-Conde et al. (2009) that were calculated for the case of similar = 0.116. Sanchez-Conde et al. (2009) presented graphs for boost factor (Fig. 7 of their work) for the case of Primack et al., (2005) EBL model, as well as for the Kneiske et al. (2004) EBL model. In the former case, the values of in the optically-thick region are typically smaller than in the latter case. These two cases (weak and strong flux enhancement, respectively) will serve to qualitatively demonstrate more or less strong -ALP mixing effects. However, in Fig. 12 we actually use fixed EBL model G12.

More detailed study of -ALP mixing is underway and will be reported elsewhere. Given huge uncertainties of EGMF strength and a vast room for the -ALP mixing parameter values, we find it appropriate to use these estimates, even though they were performed for the case of a slightly different redshift and different EBL models.

It is interesting that the case of weak mixing is practically indistinguishable from the absorption-only model fit. The primary spectrum in the case of the exotic model is indeed less hard than without it, but that does not induce any significant observable effect in the 300 — 10 energy range. Moreover, the difference between the spectra in these two cases is clearly insufficient to choose between the models on purely astrophysical grounds. On the other hand, in the low-energy range, together with a spectral irregularity at =20-30 that was qualitatively discussed in Introduction, another feature is present: namely, the spectral slope for the case of exotic model is steeper than the one for the absorption-only fit, reflecting the shape of the primary spectrum. Therefore, it appears that for the case of weak mixing the observable spectrum is the same or even steeper at all energies from the above-mentioned irregularity up to 10 , and never harder. Up to our knowledge, this highly unexpected result was never reported before.

For the case of strong mixing the situation is different — the intensity of the observable model spectrum at 3 is greater than the ones in two other cases, while all other features of the spectrum are qualitatively similar. However, according to Ajello et al. (2016), this strong mixing regime was strongly constrained considering the non-detection of the low-energy spectral irregularity, as was discussed in Introduction.

The fits for the case of electromagnetic cascade model are presented in the middle panels of Fig. 12 for the case of VERITAS (left) and HESS (right) observations. In both cases the ankle feature, which is formed by the intersection of the primary (absorbed) and secondary (cascade) components, is clearly seen. The position of the ankle is similar for both cases, again indicating that this feature, however faint, is not just a statistical fluctuation. For comparison, cascade spectra for the case of monoenergetic primary injection with various energies are also shown. Finally, low panels of Fig. 14 were computed for the case of the basic (left) and modified (right) hadronic cascade models. Both VERITAS and HESS data are presented in the picture, with normalization to the HESS data.

These fits were obtained using exactly the same procedure as those presented in Fig. 12, middle-right and bottom-left. In Appendix D, Fig. 31 we present a comparison graph of spectra calculated for basic hadronic model by Essey et al. (2011), Murase et al. (2012), and by us for the case of =0.14.

For this source we also present several fits for the case of a simple model of structured EGMF. Namely, following Furniss et al. (2015), we use the voidiness parameter defined as the fraction of the (comoving) line-of-sight distance covered by underdense regions of space (voids) to the total line-of-sight distance from the source to the observer. In fact, a spatial region that is comparatively close to the source (typically 10-100 from the source) is the most important for the cascade development, for it contains most of cascade electrons. EGMF strength in voids may be sufficiently low to allow cascade electrons radiate secondary photons before these electrons are strongly deflected from the line-of-sight. On the contrary, denser regions of space, as compared to voids, typically contain comparatively strong magnetic fields ; cascade electrons produced in these regions are strongly deflected and delayed, and secondary photons do not contribute to the observable spectrum. In effect, to account for possible structures with strong magnetic field, we multiply the cascade component to the suppression factor .

Several fits to the spectrum of 1ES 0229+200 for different values of from 0.2 to 1.0 are presented in Fig. 15. The corresponding graphs of boost factor for these fits are shown in Fig. 16. This figure demonstrates that the ankle signature in the total model spectrum does not disappear for , but, on the contrary, only gets more prominent as falls in the range of the considered values.

### 4.3 The case of 1ES 1101-232, 1ES 1812+304, 1ES 0414+009, and H 1426+428

1ES 1101-232 was observed by the HESS Collaboration in 2004-2005, and the result of spectral measurements was published in Aharonian et al. (2006). One year later, a reanalysis of these observations was performed, and the new result on the spectrum was presented in Aharonian et al. (2007b).

In our work we mainly use the latter results, but we also present a fit for the case of the electromagnetic cascade model for the former analysis. Fig. 17, top-left panel contains a fit to the SED of 1ES 1101-232 in the absorption-only model; top-right and middle-left panels of this figure contain the fits for the case of the electromagnetic cascade model for reanalysis and 2006 analysis, respectively.

A basic hadronic model fit is presented in Fig. 17, middle-right panel. The formal best fit for the case of the modified hadronic model is shown in Fig. 17, bottom-left. As for the case of 1ES 0347-121, primary component in this figure is very narrow, therefore, we also present another, more physically plausible fit for the same model in bottom-right panel of Fig. 17.

Similar fits for 1ES 1812+304 and 1ES 0414+009 are shown in Fig. 18. The spectrum of 1ES 1812+304 is not far from the case of universal spectrum of gamma-rays. In the framework of the basic hadronic model this corresponds to the case of only slightly less than of the source.

Finally, fits for the source H 1426+428 are shown in Fig. 19. For this object, the formal best fit for the electromagnetic cascade model is practically coincident with the case of the absorption-only model. The contribution of the cascade component to the total flux in this case is vastly subdominant, only about several percent. In this respect, it is interesting that Furniss et al. (2015) estimated 0 for this source, i.e. that there is small room on the line of sight where electromagnetic cascade could develop. Nevertheless, we also present two fits for the case of dominant cascade contribution at low energy (¡ 1 ) for completeness. The fits for the case of basic and modified hadronic model are presented as well.

## 5 Constraints on hadronic cascade models

As we have seen in the previous section, both electromagnetic and hadronic cascade models provide reasonable fits to the observed shape of the spectra of extreme blazars considered in this paper. Some additional information is required in order to favour one model over another. One of the criteria upon which it could be done is the plausibility of the source emission model in context of total normalization of a certain emission component. Tavecchio (2014) developed such a model (hereafter the T14 model) specifically for the hadronic cascade scenario. In what follows we use the parameters of blazar emission presented in Tavecchio (2014) (Table 2).

In fact, central emitting object of a blazar is circumvented by magnetic fields of an appreciable strength and spatial scale. In what follows we assume the model of the galaxy cluster turbulent magnetic field of Meyer et al. (2013) (equations (14),(15)) with parameters = 200 , = 2/3, = 1/2, coherence scale 10 (the spatial dimensions of all magnetic field domains are equal) and variable parameter (hereafter denoted simply as ).

These magnetic fields, in fact, modify the angular pattern of the proton beam, thus effectively dimming the observable emission. Initial angular pattern of the hadronic beam was set the same as the angular pattern of the leptonic component. The T14 model provided a following constraint to the total power of proton emission: it cannot be greater than the total “magnetic luminosity” . Indeed, in the context of the T14 model, protons during the acceleration are confined to the region of magnetic field, therefore, the energy density of these protons can not be higher than the energy density of magnetic field, otherwise they escape from the region, thus terminating the acceleration process. T14 introduced a parameter of the maximal proton acceleration energy, . We assume the accelerated proton spectrum shape .

It is convenient for us to count the total power of accelerated protons with respect to the power-law spectrum, introducing an additional factor

(15) |

so that . In what follows we introduce additional modification factors of the proton intensity with respect to the power-law spectrum of protons . We set = 1 and = 100 .

We estimate the flux modification factor as:

(16) |

where , is the bulk Lorentz factor of the blazar jet. The deflection angle of protons in the cluster magnetic field is estimated as:

(17) |

where =200 is the number of magnetic field domains and are partial deflections of proton in these domains:

(18) |

This approximation works reasonably well in the small deflection angle limit, 1 . As we will see, the small deflection angle assumption is also well justified.

The interpretation of the equation (16) is as follows: the jet angular pattern is broadened by the deflection of protons in the cluster magnetic field, high energy protons at effectively leave the region of acceleration and only a small number of particles is accelerated far above ; finally, the factor roughly accounts for the effect of the visibility of the second jet (counter-jet) when the deflection angle is large. The modification factor for the case of several values of the parameter is shown in Fig. 20.

After that, we calculate the spectrum of observable -rays for different model parameters. We assume that all energy of secondary electrons and photons was converted to the energy of cascade -rays. As we are interested only in setting an upper limit to the intensity of the observable -rays produced by protons on the way from the source to the observer, this assumption is justified.

In this section we consider the intermediate hadronic cascade model with different values of the parameter. We perform absolute normalization of the spectrum using the T14 model. We consider the source 1ES 0229+200 and two options for model parameters presented in T14 — the one corresponding to = 10, as well as to = 30. 1ES 0229+200 is an archetypal extreme blazar that was considered in Essey & Kusenko (2010b) as a very good candidate for the basic hadronic model. Fig. 21 presents normalized observable gamma-ray spectra for the same values of the parameter as in Fig. 20 and .

Finally, we performed an exhaustive scan over a wide range of parameters on a 200200 grid for two options of the source emission models presented in T14, labelled by the values of =10 and =30. Using a simple goodness-of-fit statistical test, for every version of observable spectrum we calculated a -value according to the prescriptions of Beringer et al. (2012). After that, these -values were converted to -values (i.e. statistical significance) according to the prescriptions of Beringer et al. (2012). Fig. 22 shows -values in color; it was drawn for the case of = 10, and Fig. 23 — for the case of = 30. All values of ¿10 were set to 10 and, thus, are displayed in red. The lower-left parts of Fig. 22–Fig. 23, filled by red color, correspond to the case when the intensity of the observable -rays is much greater than the one required by observations. Other parts of these graphs (violet and blue) denote the range of acceptable models. Finally, the upper-right parts of these graphs display the regions of parameters where the corresponding model is excluded. It is remarkable that all considered configurations with above 100 are excluded with significance ¿7 .

During this calculation, it will be remembered, we have neglected the losses of protons on EBL from the source to the observer. However, there are two main factors that strongly outweigh this one, namely: 1) injection of protons in context of the T14 model is done at an energy much lower than 1 ; therefore, for the case of = 10 and =2 the total energy of accelerated protons is 5 times greater than the energy of UHE protons, thus, the modification factor in this case gains an additional factor of 0.2 2) synchrotron and curvature radiation losses of primary UHE protons will futher diminish their total energy. The last factor is particularly strongly constraining for the case of the Blandford-Znajek emission mechanism (Blandford-Znajek (1977)) due to strong magnetic fields in the acceleration region, see Neronov et al. (2009) and Neronov & Ptitsyna (2016).

## 6 Summary and conclusions

In this work we reviewed the main possibilities of how cascades from primary -rays or protons may influence the data interpretation when testing extragalactic -ray propagation models, mainly dealing with the source redshift range 0.1–0.3. We reviewed the main regimes of electromagnetic cascade development: 1) one-generation regime, which holds when 10 , 2) the universal regime, which is satisfied when 100 3) a possibility of the extreme energy cascade regime at 1 .

We developed a fast, simple and suffuciently precise hybrid code that allows to calculate the observable spectrum of -rays for the case of primary proton. We performed calculations of observable spectra for this case and investigated several versions of the hadronic cascade model — the basic model (all gamma-rays are generated by primary protons), the intermediate model (the same, but the proton beam is dissolved at some , the basic model is the sub-set of intermediate models with = 0), the modified model that includes the primary component. We discussed the signatures of the spectrum for the basic hadronic model for = 0.186. The signatures at a comparatively low observable energy are almost the same as the ones for the case of a purely electromagnetic cascade — that is, the spectrum below and spectrum from to . However, at higher energy (200 ) the spectrum in the hadronic model is much harder than the universal spectrum for the case of a pure electromagnetic cascade.

We also performed the comparison of two hybrid calculations of observable spectra for different values of from 0 up to nearly . The results of these calculations with two different methods are in very good agreement for the case of the proton primary energy 30 , i.e. when the pair-production losses are dominant. For the case of = 100 the agreement between the two methods is also good for the basic hadronic model and the intermediate hadronic model with comparatively low , but not so good if . The reason of this disagreement may be connected to the possibility of the extreme energy cascade regime. New calculations and measurements of the universal radio background are necessary in order to decide which spectrum is more realistic.

While we confirmed the claim that the shape of the observable spectrum in the basic hadronic model depends only slightly on the primary proton energy, it was found that in context of the intermediate hadronic model the shape of the spectrum strongly depends on the parameter in the high energy region.

We performed a series of fits to the spectra of extreme blazars with various extragalactic -ray propagation models, namely: the absorption-only model, the electromagnetic cascade model, the basic hadronic model, the intermediate hadronic model, the modified hadronic model. Most of the fits presented in this paper are formal best fits. It was found that, as a rule, both the electromagnetic cascade model and the hadronic cascade model are able to provide a reasonable fit to the observed spectral shape.

For the case of blazar 1ES 0347-121 and the electromagnetic cascade model we performed the calculation of modification factor with respect to the absorption-only model. We found that ¿1 in the energy region of 1-10 . Therefore, the version of the electromagnetic cascade model that was considered in this work for this source predicts that it is possible that the observable intensity is higher than in the absorption-only model for the same level of EBL without any new physics. For the case of the source 1ES 0229+200 we also performed a fit in the framework of an exotic model with oscillations and a calculation of the modification factor vs. voidiness dependence in the framework of the electromagnetic cascade model. We found that when is lower than 1, the value of does not fall at high energies, but, on the contrary, even grows in the energy range of 1-15 , thus making the effects induced by electromagnetic cascades even more pronounced.

Finally, using the fact that many blazars are situated in galaxy groups or clusters with comparatively strong magnetic fields (Muriel (2016), Oikonomou et al. (2014)), we performed the testing of the emission model of Tavecchio (2014). We found that for ¿100 this emission model is excluded with significance 7 . Our results show that the hadronic cascade model experiences significant difficulties. This conclusion is in line with results obtained by Razzaque et al. (2012). Further testing of the hadronic and electromagnetic cascade models may be performed on an event-by-event level with advanced analysis tools such as GammaLib and ctools package (Knoedlseder et al. (2016)). Such a work is underway in our group.

We conclude that cascades from primary -rays or nuclei may induce an appreciable background for astrophysical searches for oscillation that needs to be taken into account. The electromagnetic and hadronic cascade models deserve further study, also in context of searches for other exotic phenomena, such as Lorentz invariance violation (Tavecchio & Bonnoli (2016)). The knowledge of EGMF would provide a significant insight into the intergalactic cascade phenomena in the Universe. Further measurements with the CTA array of Cherenkov telescopes (Acharya et al. (2013)) will likely clarify these issues. Together with existing instruments Fermi LAT (Atwood et al. (2009)) and AGILE (Rappoldi et al. (2016)), future experiments that have either high sensitivity or good angular resolution, such as GAMMA-400 (Galper et al. (2013)) or the novel balloon-borne emulsion -ray telescope GRAINE (Takahashi et al. (2015)), may also prove to be helpful in this task.

###### Acknowledgements.

This work was supported by RFBR, Grant 16-32-00823. The authors are grateful to Dr. V.V. Kalegaev for permission to use the SINP MSU space monitoring data center computer cluster and to V.O. Barinova, M.D. Nguen, D.A. Parunakyan for technical support.## References

- (2012) Hillas, A.M. 2013, APh, 43, 19
- (1992) Punch, M. et al. 1992, Nature, 358, 477
- (1999) Aharonian, F. et al. 1999, A&A, 349, 11
- (1993) Stecker, F.W. & de Jager, O.C. 1993, ApJ, 415, 71
- (1994) de Jager, O.C., Stecker, F.W. & Salamon, M.H. 1994, Nature, 369, 294
- (1962) Nikishov, A.I. 1962, Sov. Phys. JETP, 14, 393
- (1967) Gould, R.J. & Shreder, G. 1967, Phys. Rev., 155, 1408
- (1966) Jelley, J.V. 1966, Phys. Rev. Lett., 16, 479
- (2012) Ackermann, M. et al. 2012, Science, 338, 11
- (2013) Abramowski, A. et al. (The H.E.S.S. Collaboration). 2013, A&A, 550, A4
- (2012) Horns, D. & Meyer, M. 2012, JCAP, 033
- (2015) Biteau, J. & Williams, D.A. 2015, ApJ, 200, 58
- (2016) Horns, D. 2016, preprint astro-ph/1602.07499
- (2000) Protheroe, R.J. & Meyer, H. 2000, Phys. Lett. B, 493, 1
- (2013) Meyer, M. et al. 2013, Phys. Rev. D, 87, 035027
- (2007) Andriamonje, S. et al. (CAST Collaboration). 2007, JCAP, 04, 010
- (2009) Sanchez-Conde, M.A. et al. 2009, Phys. Rev. D, 79, 123511
- (2016) Ajello, M. et al. 2016, Phys. Rev. Lett., 116, 161101
- (2014) Ayala, A. et al. 2014, Phys. Rev. Lett., 113, 191302
- (2013) Abramowski, A. et al. (H.E.S.S. Collaboration). 2013, Phys. Rev. D, 88, 102003
- (2015) Payez, A. et al. 2015, JCAP 2, 006
- (2013) Wouters, D. & Brun, P. 2013, ApJ, 772, 44
- (1999) Blasi, P. et al. 1999, ApJ, 514, L79
- (2016) Pshirkov, M.S. et al. 2016, Phys. Rev. Lett., 116, 191302
- (2005) Dolag, K. et al. 2005, JCAP, 01, 009
- (2010) Neronov, A. & Vovk, Ie. 2010, Science, 328, 73
- (2011) Dermer, C.D. et al., 2011, ApJ Lett., 733, L21
- (2011) Taylor, A.M. et al. 2011, A&A, 529, A144
- (2012) Vovk, Ie. et al. 2012, ApJ Lett., 747, L14
- (2014) Abramowski, A. et al. (H.E.S.S. Collaboration). 2014, A&A, 562, A145
- (2012) Takahashi, K. et al., ApJ Lett., 744, L7 (2012)
- (2013) Takahashi, K. et al., ApJ Lett., 771, L42 (2013)
- (2010) Akahori, T. & Ryu, D. 2010, ApJ, 723, 476
- (2015a) Dzhatdoev, T.A. 2015, J. Phys. Conf. Ser., 632, 012035
- (2015) Finke, J.D. et al. 2015, ApJ, 814, 20,
- (2014) Arlen, T.C. et al. 2014, ApJ, 796, 18
- (2015) Tashiro, H.& Vachaspati, T. 2015, MNRAS, 448, 299
- (2009) Atwood, W.B. et al. 2009, ApJ, 697, 1071
- (2012) Neronov, A. et al. 2012, A&A, 541, A31
- (2011) Abdo, A. et al. 2011, ApJ, 727, 129
- (2015) Furniss, A. et al. 2015, MNRAS, 446, 2267
- (2015) Chen, W. et al. 2015, Phys. Rev. Lett., 115, 211103
- (2002) Aharonian, F.A. et al. 2002, A&A, 384, 834
- (2007) dâAvezac, P. et al. 2007, A&A, 469, 857
- (2012) Murase, K. et al. 2012, ApJ, 749, 63
- (2013) Takami, H. et al. 2013, ApJ Lett., 771, L32
- (2015b) Dzhatdoev, T. 2015, Izvestiya Rossiiskoi Akademii Nauk, 79, 363 (Preprint astro-ph/1501.00259)
- (1998) Uryson, A., 1998. JETP, 86, 213
- (2010a) Essey, W. & Kusenko, A. 2010, APh, 33, 81
- (2010b) Essey, W. et al. 2010, Phys. Rev. Lett., 104, 141102
- (2011) Essey, W. et al. 2011 ApJ 731 51
- (2014) Essey, W. & Kusenko, A. 2014, APh, 57, 30
- (2016) Zheng, Y.G et al. 2016, A&A, 585, A8
- (2016) Berezinsky, V. & Kalashev, O. 2016, Phys. Rev. D, 94, 023007
- (1997) Brun, R. & Rademakers, F. 1997, NIM A, 389, 81
- (2012) Kachelriess, M. et al. 2012, Comp. Phys. Comm., 183, 1036
- (2010) Kneiske, T.M. & Dole, H. 2010, A&A, 515, A19
- (2012) Broderick, A. E. et al. 2012, ApJ, 752, 22
- (2012) Schlickeiser, R. et al. 2012, ApJ, 758, 102
- (2014) Chang, P. et al. 2014, ApJ, 797, 110
- (2015) Menzler, U. & Schlickeiser, R. 2015, MNRAS, 448, 3405
- (2013) Miniati, F. & Elyiv, A. 2013, ApJ, 770, 54
- (2013) Venters, T.M. & Pavlidou, V. 2013, MNRAS, 432, 3485
- (2014) Sironi, L. & Giannios, D. 2014, ApJ, 787, 49
- (2016) Kempf, A. et al. 2016, A&A, 585, A132
- (2006) Berezinsky, V.S., et al. 2006, Phys. Rev. D, 74, 043005
- (2015) Ade, P.A.R. et al. (Planck Collaboration) 2016, A&A, 594, A13
- (2008) Kelner, S.R. & Aharonian, F.A. 2008, Phys. Rev. D, 78, 034013
- (1970) Blumenthal, G.R. 1970, Phys. Rev. D, 1, 1596
- (1970) Blumenthal, G.R. & Gould, R.J. 1970, Rev. Mod. Phys., 42, 237
- (2008) Khangulyan, D. et al. 2008, MNRAS, 383, 467
- (2012) Aharonian, F. et al. 2012, A&A, 547, A114
- (1996) Protheroe, R.J. & Biermann, P.L. 1996, APh, 6, 45
- (2015) Settimo, M. & De Domenico, M. 2015, APh, 62, 92
- (1973) Brown, R.W. et al. 1973, Phys. Rev. D, 8, 3083
- (2009) Demidov, S.V. & Kalashev, O.E. 2009, JETP, 108, 764
- (1994) Mastichiadis, A. et al. 1994, MNRAS, 266, 910
- (2015) Bonnoli, G. et al. 2015, MNRAS, 451, 611
- (2003) Aharonian, F. et al. 2003, A&A, 403, 523
- (2002) Djannati-Atai, A. et al. 2002, A&A, 391, L25
- (2002) Horan, D. et al. 2002, ApJ, 571, 753
- (2007a) Aharonian, F. et al. 2007, A&A, 475, L9
- (2014) Aliu, E et al. 2014, ApJ, 782, 13
- (2013) Madhavan, A.S. et al. 2013, preprint astro-ph/1307.7051
- (2007b) Aharonian, F. et al. 2007, A&A, 470, 475
- (2006) Aharonian, F.A. et al. (H.E.S.S. Collaboration). 2006, Nature, 440, 1018
- (2007c) Aharonian, F. et al. 2007, A&A, 473, L25
- (2012) Abramowski, A. et al. (The H.E.S.S. Collaboration). 2012, A&A, 538, A103
- (2016) Wakely, S. & Horan, D. http://tevcat.uchicago.edu/
- (2010) Acciari, V.A. et al. (VERITAS Collaboration). 2010, ApJ Lett., 709, L163
- (2012) Meyer, M. et al. 2012, preprint astro-ph/1211.6405
- (2008) Franceschini, A. et al. 2008, A&A, 487, 837
- (2011) Dominguez, A. et al. 2011, MNRAS, 410, 2556
- (2012) Gilmore, R.C. et al. 2012, MNRAS, 422, 3189
- (1975) James, F. & Roos, M. 1975, Comput. Phys. Commun., 10, 343
- (2005) Primack, J. et al. 2005, preprint astro-ph/0502177
- (2004) Kneiske, T.M. et al. 2004, A&A, 413, 807
- (2014) Tavecchio, F. 2014, MNRAS, 438, 3255
- (2012) Beringer, J. et al. (PDG). 2012, Phys. Rev. D 86, 010001
- (1977) Blandford, R.D. & Znajek, R.L. 1977, MNRAS, 179, 433
- (2009) Neronov, A.Yu. et al. 2009, New Journal of Physics, 11, 065015
- (2016) Ptitsyna, K. & Neronov, A. 2016, A&A, 593, A8
- (2016) Muriel, H. 2016, A&A, 591, L4
- (2014) Oikonomou, F. et al. 2014, A&A, 568, A110
- (2012) Razzaque, S. et al. 2012, ApJ, 745, 196
- (2016) Knoedlseder, J. et al. 2016, A&A, 593, A1
- (2016) Tavecchio, F. & Bonnoli, G. 2016, A&A, 585, A25
- (2013) Acharya, B.S. et al. 2013, APh, 43, 3
- (2016) Rappoldi, A. et al. 2016, A&A, 587, A93
- (2013) Galper, A.M. et al. 2013, AIPCP, 1516, 288
- (2015) Takahashi, S. et al. 2015, Prog. Theor. Exp. Phys., 043H01

## Appendix A Electromagnetic cascades in the universal regime

Here we present some additional plots to reveal the borders of the universal regime of EM cascade development. Fig. 24–25 show the spectra for the case of two different redshifts, different primaries (electron and -ray), and two primary energies. We have checked that for the case of greater energies, 10 and 100 , the spectra for both redshifts and both primaries are practically coincident (within the width of the line) with ones shown in these figures for the case of =1 . Fig. 26–26 show the ratio of the spectra to the spectrum for the case of primary -ray and = 1 for different primaries and different , as well for two values of .

## Appendix B Absorbed and cascade components for various primary spectra of -rays

In Fig. 28 we present several other examples of calculations for the case of a purely electromagnetic cascade and the following primary spectrum:

(19) |

where = 1 if and 0 otherwise. Fig. 28 top-left graph is drawn for the case of = 1 and = 100 ; Fig. 28 top-right — for = 2 and = 100 , Fig. 28 middle-left — = 3 and = 100 , Fig. 28 middle-right — = 2 and = 10 , Fig. 28 bottom-left — for = 3 and = 10 , Fig. 28 bottom-right — for = 2 and = 1 .

## Appendix C Optical depth in the ELMAG code

Fig. 29 shows the comparison between several models of optical depth with the version of the KD10 EBL model as implemented in the ELMAG code. Other EBL models include G12, KD10 (original implementation) and F08.

## Appendix D Comparison of observable -ray spectra with other studies

Here we present the comparison of our test calculations with other works. Fig. 30 shows good agreement between our calculations with the ELMAG code (assuming the KD10 EBL model, as implemented in ELMAG) and the result of Vovk et al. (2012), obtained with the F08 EBL model. Fig. 31 presents the comparison of our result for the case of the basic hadronic cascade model with the ones obtained by Essey et al. (2011) and Murase et al. (2012).