Modeling Gamma-Ray SEDs and Angular Extensions of Extreme TeV Blazars from Intergalactic Proton-Initiated Cascades in Contemporary Astrophysical EGMF Models

The intrinsic spectra of some distant blazars known as “extreme TeV blazars” have shown a hint at an anomalous hardening in the TeV energy region. Several extragalactic propagation models have been proposed to explain this possible excess transparency of the Universe to gamma-rays starting from a model which assumes the existence of so-called axion-like particles (ALPs) and the new process of gamma-ALP oscillations. Alternative models suppose that some of the observable gamma-rays are produced in the intergalactic cascades. This work focuses on investigating the spectral and angular features of one of the cascade models, the Intergalactic Hadronic Cascade Model (IHCM) in the contemporary astrophysical models of Extragalactic Magnetic Field (EGMF). For IHCM, EGMF largely determines the deflection of primary cosmic rays and electrons of intergalactic cascades and, thus, is of vital importance. Contemporary Hackstein models are considered in this paper and compared to the model of Dolag. The models assumed are based on simulations of the local part of large-scale structure of the Universe and differ in the assumptions for the seed field. This work provides spectral energy distributions (SEDs) and angular extensions of two extreme TeV blazars, 1ES 0229+200 and 1ES 0414+009. It is demonstrated that observable SEDs inside a typical point spread function of imaging atmospheric Cherenkov telescopes (IACTs) for IHCM would exhibit a characteristic high-energy attenuation compared to the ones obtained in hadronic models that do not consider EGMF, which makes it possible to distinguish among these models. At the same time, the spectra for IHCM models would have longer high energy tails than some available spectra for the ALP models and the universal spectra for the Electromagnetic Cascade Model (ECM). The analysis of the IHCM observable angular extensions shows that the sources would likely be identified by most IACTs not as point sources but rather as extended ones. These spectra could later be compared with future observation data of such instruments as Cherenkov Telescope Array (CTA) and LHAASO.


Introduction
Besides having many features that attract scientific interest and prompt their continuous study, blazars are also used as standard probes to determine the opacity of the Universe to gamma-rays. Gamma-rays from distant blazars are subject to absorption on the photons of Extragalactic Background Light (EBL) through the pair production process. This means that the observable Spectral Energy Distributions (SEDs) of gamma-rays from such sources would be attenuated and that one can obtain the intrinsic SEDs of the sources by correcting for the gamma-ray absorption. In this work extreme TeV blazars are defined as blazars, which have intrinsic SEDs peaking at TeV energies. By now, approximately 10 extreme TeV blazars have been registered [1]. Extreme TeV blazars, as compared to such blazars as Mkn 501 and Mkn 421, typically show weak and slow variability for the high energy and very high energy bands, which makes them suitable for studying the effects of intergalactic magnetic and photon fields.
Authors of [2] have obtained intrinsic SEDs for blazars of various classes (some of which were extreme TeV blazars) and found that for a number of sources there is an apparent excess of gamma-rays in the optically thick region of the spectra. The authors of [3] have also conducted a search for this anomaly, favoring its presence with a significance of ∼ 2σ.
Several gamma-ray propagation models have been developed to account for the possible hardening of the intrinsic spectra. One of such models introduces the so-called axion-like particles (ALPs) which are assumed to be pseudoscalar bosons with extremely low masses that oscillate into photons and vice versa in a similar way to how different flavors of neutrinos oscillate [4][5][6][7][8][9][10][11]. The theory of gamma-ALP mixing was first described in [12] and later developed in [13][14][15][16][17][18]. The authors of [2] credited ALPs as one of the possible explanations to their high energy anomaly. Indeed, if gamma-rays convert to ALPs in the presence of Extragalactic Magnetic Field (EGMF) and propagate an appreciable distance as ALPs (which are not absorbed on EBL photons) before converting back into gamma-rays near the observer, these gamma-rays experience a lesser degree of absorption on EBL photons, which can explain the hardening of the intrinsic blazar spectra. ALPs are characterized by two parameters: their mass m ALP and the gamma-ALP coupling constant g γ−ALP . The latest existing constraints on these parameters are described in [19] (see also [20][21][22][23]).
Alternatively, the excess of photons could be accounted for in the so-called cascade models. These models presuppose that a part of gamma-rays is produced not inside the source but rather in the intergalactic cascades as secondary particles [24][25][26][27]. In the Electromagnetic Cascade Model (ECM) [28,29] gamma-rays are considered to be the primary particles, and the cascades are initiated by the inverse Compton scattering of electrons produced in primary gamma-ray interactions with EBL photons on the photons of Cosmic Microwave Background (CMB). The other type of cascade models, Intergalactic Hadronic Cascade Model (IHCM) [30,31], assumes Ultra-High Energy Cosmic Rays (UHECR) as primary particles that initiate cascades via Bethe-Heitler pair production and photohadronic processes on various photon fields (CMB, EBL, Universal Radio Background (URB)).
For cascade models the account of state-of-the-art models of EGMF is vital because in these models stronger magnetic fields deflect primary CRs and electrons of intergalactic cascades, which has a considerable impact on both SED and angular extension of the source. Without the account of EGMF, the observable SEDs in the IHCM (in this case called "basic hadronic cascade model") inside a typical IACT PSF would have longer high energy tails because in this case primary protons and cascade electrons would not be deflected by the EGMF. An extensive research of basic hadronic cascade model has been done by the author of this paper and his co-authors in [32]. For this work I use and compare several of the latest EGMF models based on MHD simulations following the large-scale structure of the Universe [33,34] in the IHCM. For ALP models stronger magnetic fields mean smaller values of the mixing constant g γ−ALP and, thus, higher sensitivities to ALPs. Thus, the use of astrophysically accurate EGMF models could help constrain the parameter space of g γ−ALP and facilitate the search for gamma-ALP oscillations. This work does not consider the ALP models, although a similar simulation of these models would be beneficial.
In this work I provide model SEDs and angular extensions of observable gamma-rays from two extreme TeV blazars in the IHCM using several models of EGMF. Since the possible excess of observable gamma-rays in the TeV region of blazar spectra is one of the problems that will be investigated by the soon-to-be-operational Cherenkov Telescope Array (CTA) [35], these results would provide a valuable model comparison to the future CTA observations of such sources. The collaboration of the LHAASO particle detector [36] has also included some extreme TeV blazars in its scientific program.
The paper is organized as follows. In Section 2 the contemporary models of EGMF are discussed and the choice of models used for this work is specified. In Section 3 the simulation methods are described, including the choice of modeling codes, photon field models, and sources. In Section 4 the main results of this paper are presented and discussed. Finally, in Section 5 the main conclusions of this paper are drawn, its limitations are outlined, and prospects for future research are discussed.

EGMF Models
There have been several works that model the EGMF with different assumptions for the seed magnetic field [33,34,[37][38][39][40][41]. Most of these models are assume that the seed magnetic fields were generated during the formation of the large-scale structure (i.e., astrophysical models) and use magneto-hydrodynamical simulations to model the evolution of these fields.
Three of these models, Refs. [33,34,39] are being frequently used for various astrophysical research, such as extragalactic CR propagation and deflection studies. An important difference between the models of Sigl et al. [39] and Dolag et al. [34] (hereafter named D05 for simplicity) is that in [39] the seed magnetic field is generated from a zero field through the Biermann battery, while in [34] it is assumed as a uniform field with strength B seed ∼ 10 −9 G. Moreover, Ref. [34] is a model of the local part of the large-scale structure with Milky Way in the center constrained by the astrophysical data.
The models of Hackstein et al. [33] (hereafter named H18 for simplicity) improve on the D05 by using more advanced numerical methods and initial conditions of [42], which are based on more contemporary observational constraints. H18 include six different models of EGMF: three primordial ones (i.e., models in which the seed magnetic field is introduced in the early Universe, e.g., during the inflation or phase transition epochs) and three astrophysical ones. The astrophysical models assume a seed magnetic field of B seed ∼ 10 −11 nG and impulsive thermal and magnetic feedback in haloes as the mechanism of EGMF generation. These models differ in the energy budget assumptions with astrophysical model assuming the release of 5 × 10 58 erg per feedback episode starting from redshift z = 4, astrophysical R model-10 60 erg per feedback episode from z = 4, and astrophysical 1R model-an energy budget changing from 10 60 to 5 × 10 58 per feedback episode starting from z = 1. A recent work by [43] has revealed that astrophysical R and astrophysical 1R show almost the same dependence of mean deflections of CRs on their rigidity. For these models the deflections are strongest among all considered astrophysical models expect [39], followed closely by D05. The astrophysical model demonstrates the lowest mean deflection values.
This work considers three astrophysical EGMF models: D05, and astrophysical and astrophysical R of H18.

Simulation Methods
The simulation was done in several steps. Firstly, the publicly available CRPropa3 simulation framework [44] was used to model the propagation of primary protons and their deflection on the filaments of EGMF. CR-Propa3 allows to load the local cubic volumes of D05 and H18 magnetic fields with edge lengths of 132 Mpc and 249 Mpc respectively. For the simulations 10 5 protons with primary energies 30 EeV were isotropically emitted from random points throughout the simulated EGMF volume. The replication of the EGMF volume allowed for propagating these protons to large distances (up to 1200 Mpc) and tracking their energies, positions, and directions every 5 Mpc of the way. Throughout the proton propagation the code accounted for adiabatic as well as interaction losses. The interactions included Bethe-Heitler pair production on CMB and EBL (using the model by Gilmore et al. [45]) photon fields and pion production on CMB, EBL, and URB (using the model by Protheroe et al. [46]) photon fields. As a result, the deflections of 10 5 protons were calculated for every 5 Mpc of their way.
The next step of the simulation was to calculate the observable spectra for gammarays produced in the intergalactic cascades. For that the underlying assumption of the so-called weak universality was used as introduced in [47] and previously used by the author of this paper in [32]. As explained in [32], weak universality means that "the spectrum of observable gamma-rays is practically independent of the energy and the type, photon/electron, of the primary particle but depends on the redshift of the source". Thus, observable cascade gamma-ray spectra could be calculated for a range of redshifts, and an SED of a source with redshift z s could be obtained by summing the cascade SEDs from all the intermediary redshifts from 0 to z s .
To this end the publicly available ELMAG 2.03 code [48] was used to compute a database of the universal EM cascade spectra with primary photon energy 1 PeV and redshift z distributed randomly and uniformly varying from 0 to 0.30. The EBL model [45] was used for these calculations. This code was used in the previous works of this author, and its applicability and sensitivity of the results to the choice of EBL models are discussed in [32]. In this work the deflections of secondary particles in voids were neglected. If they were taken into account, it would broaden the observable angular extensions of the sources as compared to the results below.
The final step of the simulation utilized an improved version of the hybrid code developed by the author of this paper and his colleagues and used previously in [32]. This code propagated the protons from the source to the observer using semi-analytic formulas that calculate adiabatic losses of protons, their pair production losses on CMB [49], and pion production losses on both CMB and EBL [50] and combining them with Monte Carlo simulations of precomputed cascade spectra and proton deflections in EGMF described above. As an output, this code produced a file with coordinates, directions, energies, and observable cascade spectra of 10 4 protons on every redshift step dz = 10 −4 . This file was then used to produce the averaged SEDs and model angular extensions of the sources shown in the next section.
Two extreme TeV blazars were chosen as sources for this work because of their low variability. They were 1ES 0229+200 with redshift z = 0.14 and 1ES 0414+009 with redshift z = 0.287. The spectrum of the first source was measured by both the H.E.S.S. Collaboration [51] (with threshold energy E thr = 580 GeV, total exposure time 41.8 h, and mean zenith angle Z mean = 46 • ) and the VERITAS Collaboration [52] (with threshold energy E thr = 300 GeV and total exposure time 54.3 h). The spectrum of the second source was measured by H.E.S.S. [53] (with threshold energy E thr = 200 GeV, total exposure time 73.7 h, and zenith angles varying between 22 • and 41 • with a mean value of Z mean = 26 • ).

Results
In order to see how the observable SEDs in the framework of the IHCM would look like inside the Point Spread function (PSF) of observing gamma-ray telescopes the following cut was applied to the value of the observable angle (i.e., the angle between the direction of an incident cascade gamma-ray and the line source-observer): θ obs < 0.1 • . The value 0.1 • represents a typical average value of PSFs of IACTs like H.E.S.S. Figure 1 shows these observable SEDs (denoted as blue curves) in comparison with the integral observable SEDs (black curves), for which no cuts on the observable angles were applied (these SEDs correspond to the basic hadronic cascade model, which assumes a zero EGMF), and the universal SEDs of ECM, where gamma-rays were the primary particles (red curves). One can see that the model SEDs from D05 and astrophysical R EGMF models demonstrated a significant attenuation in comparison with the integral spectra, especially towards the higher energies. This attenuation was the strongest for the D05 model and the weakest for the astrophysical model. Figure 2 demonstrates a comparison between IHCM model SEDs for the source 1ES 0229+200 (z = 0.14) and all considered EGMF models and ALP model SEDs obtained in [14] for a similar z = 0.116, ALP masses m ALP = 10 −10 eV and two EBL options: Primack et al. [54] (P05) and Kneiske et al. [55] (K04). All IHCM models were normalized at 1 TeV. It can be seen that IHCM SEDs in this case had longer high energy tails than ALP SEDs. However, this result was strongly dependent on the assumed parameters of the ALPs. A more thorough study examining different ALP models with varied parameters is underway and will be published elsewhere. Black curves denote the integral SEDs or SEDs for the basic hadronic cascade model (i.e., SEDs, for which no cuts on the observable angles have been applied). Red curves denote the universal spectra in the framework of the Electromagnetic Cascade Model (ECM). (a,c,e) demonstrate the results for the source 1ES 0229+200 (z = 0.14) and D05, astrophysical (H18), and astrophysicalR (H18) models of Extragalactic Magnetic Field (EGMF) respectively; (b,d,f)-for the source 1ES 0414+009 (z = 0.287) and D05, astrophysical (H18), and astrophysical R (H18) models of EGMF respectively.  . Comparison of IHCM model SEDs of the blazar 1ES 0229+200 (z = 0.14) with ALP SEDs from [14] (z = 0.116) and IACT measurements. Solid lines denote IHCM model SEDs: black line corresponds to astrophysical R (H18) model of EGMF, magenta line-to astrophysical (H18), and olive line-to D05. Green dashed curve corresponds to ALP model SEDs assuming P05 EBL model [54], blue dashed curve-to ALP model SEDs assuming K04 EBL model [55]. Dotted lines represent primary spectra in ALP models. Finally, symbols with error bars denote IACT measurements of the observable SEDs of 1ES 0229+200: red dots correspond to VERITAS measurements, red triangles-H.E.S.S. measurements. Figure 3 demonstrates the observable angular extensions of sources in the framework of the IHCM in the form of containment angle dependence on observable energy. Different curves denote different quantiles or containment angle vs. observable energy dependences for different containment values (e.g., black solid curve denotes an angle, inside which 5% of the observable spectrum was contained). One can see two distinct patterns in these graphs. For the D05 model, all the values of containment angles increased substantially with energy for E > 10 12 eV so that 68% of all observable angles were greater than 0.73 • at E = 10 TeV for the source 1ES 0229+200 (z = 0.14) and 2.35 • for the source 1ES 0414+009 (z = 0.287) and the same energy. Containment angle (f) Figure 3. Containment angle vs. observable energy dependences for different containment values in the framework of the IHCM. Black solid curve denotes a 5% containment angle, red solid curve-10% containment angle, green solid curve-20%, blue solid curve-40%, black dashed curve-68%, red dashed curve-80%, green dashed curve-90%, and blue dashed curve-95%. (a,c,e) demonstrate the results for the source 1ES 0229+200 and D05, astrophysical (H18), and astrophysicalR (H18) models of EGMF respectively. (b,d,f)-for the source 1ES 0414+009 and D05, astrophysical (H18), and astrophysical R (H18) models of EGMF respectively. A different pattern was noticeable for H18 models. While the upper quantiles of containment angles rose rather quickly with energy, the lower ones (less than 40% for the astrophysical model and less than 10% for the astrophysical R model) showed only a very slight increase in the values of containment angles. This could mean that in these models a sizeable part of protons did not encounter any filaments on their way and, thus, did not experience strong deflections. So, the deflection effect of filaments for the H18 models was stronger, while at the same time the filaments themselves could be thinner or more far between in comparison with the D05 model.
It can also be said that the characteristic angular extension of these sources (for instance, the value of θ 68 ) mostly fell outside the PSFs of the gamma-ray telescopes, which is clearly seen in Figure 4 (showing PSFs of Fermi-LAT, CTA, H.E.S.S., MAGIC, and LHAASO). For this figure angular resolutions for Fermi-LAT were taken from [56], for H.E.S.S., MAGIC and CTA-from [35], for LHAASO-from [57]. The angular resolution for VERITAS [58,59] was similar to that of H.E.S.S. and MAGIC. The results show that these telescopes would see these sources as extended ones as compared to point sources. Only for the astrophysical model of EGMF and the source 1ES 0229+200 (z = 0.14) the value of θ 68 is lower than the PSFs of these telescopes for energies E < 0.7 TeV. Given that the source 1ES 0229+200 has been reported as a point source based on IACT observations [51], the angular extension predicted by the IHCM diminished this model's plausibility. Further statistical analysis is required to constrain or possibly disfavor the intergalactic hadronic cascade model.  These SEDs and angular extensions could later be compared with future observation data by imaging atmospheric Cherenkov telescopes (IACTs) such as CTA or particle detectors such as LHAASO. For example, Figure 5 shows the SEDs measured by H.E.S.S. and VERITAS with the model SEDs of the three models of EGMF considered in this paper inside the PSF of the LHAASO particle detector (normalized to fit the data using the same normalization constant) as well as the 1-year and 5-year sensitivity plots of LHAASO obtained in [36] with statistical significance of 5 standard deviations. The PSF of the LHAASO detector was taken from [57]. Since the PSF was available only for energies E > 10 TeV, an integral SED was used for lower energies. It can be seen that LHAASO particle detector had an outstanding sensitivity at energies >10 TeV allowing it to register the hard tails in the SEDs of certain blazars like 1ES 0229+200, which were included in its scientific program. . Comparison of IHCM model SEDs of the blazar 1ES 0229+200 inside the PSF of LHAASO experiment (the PSF was taken from [57] for energies E > 10 TeV; for lower energies an integral SED is used) with SEDs measured by H.E.S.S. [51] and VERITAS [52] as well as LHAASO sensitivity curves. Solid lines denote model SEDs: blue line corresponds to astrophysical R (H18) model of EGMF, cyan line-to astrophysical (H18), and black line-to D05. Dashed lines denote 1-year (green line) and 5-year (magenta line) sensitivities of the LHAASO experiment particle detector. Finally, symbols with error bars denote IACT measurements of the observable SEDs of 1ES 0229+200: red dots correspond to VERITAS measurements, red triangles-H.E.S.S. measurements.

Discussion and Conclusions
This work presents SEDs and angular extensions of observable gamma-rays from two extreme TeV blazars, 1ES 0229+200 (z = 0.14) and 1ES 0414+009 (z = 0.287), in the framework of the Intergalactic Hadronic Cascade Model for several contemporary models of EGMF (astrophysical and astrophysical R of H18) and D05.
There is a number of features that the observable SEDs and angular extensions demonstrate, which could be used to distinguish between different propagation models. In the IHCM there is a visible attenuation of observable SEDs inside the PSF of observing instruments revealing itself at higher energies (as compared to hadronic models neglecting the EGMF). This attenuation is more pronounced for D05, where at E = 100 TeV the difference between the integral SEDs and the SEDs inside a PSF of 0.1 • has a factor of ∼80 for the source 1ES 0229+200 and ∼220 for the source 1ES 0414+009. The ALP models exhibit one important feature: a step-like irregularity at lower energies corresponding to a sizeable drop in the intensity of photons (which could be up to 1/3 of its value) related to two photon polarization states attaining equipartition with one polarization state of ALP above a certain critical energy [14]. As of yet, this irregularity has not been observed in the spectra of several sources [20,21], which led to some of the constraints on the mixing parameters. In comparison with the ALP model SEDs from [14] the IHCM SEDs showed longer high energy tails. Further comparative analysis of these models including ALP models with varied parameter values is necessary for determining the discrimination criteria between the SEDs in IHCM and ALP models.
The analysis of angular extensions of extreme TeV blazars in the framework of IHCM demonstrated clear distinctions between D05 and H18 models with H18 models having fewer or thinner but at the same time more deflective filaments than D05. It has also shown that these sources would be mostly identified by the gamma-ray telescopes as extended sources, not the point ones. The only sources that could in principle be identified as point sources are near sources in the astrophysical model of EGMF. Since this kind of angular extension has not been supported by the IACT observational data, the author of this paper concludes that a thorough statistical analysis of angular extension is crucial for constraining or disfavoring the IHCM.
These SEDs and angular extensions could be compared with future IACT and particle detector observations. For example, the characteristics of CTA allow this telescope to investigate the gamma-ray signal from distant extreme TeV blazars with its angular resolution improving from ∼ 0.3 • at 20 GeV to ∼ 0.03 • at 100 GeV. Moreover, LHAASO's outstanding sensitivity at energies > 10 TeV allows it to observe the high-energy part in the blazar SEDs as well.
A part of cascade energy can be lost due to plasma beam instabilities [60]. There have been many works modeling these collective plasma losses using Particle-In-Cell simulation methods (they are listed and reviewed in relation to cascade pair beams in Section 4.3 of [61]). However, as of yet there is no clear consensus on the extent of these losses in the context of extragalactic gamma-ray propagation from distant blazars. Such losses are not accounted for in this work; if they are present and appreciable, they could further attenuate the observable gamma-ray spectra.
The consideration of other EBL models, primary proton energies, source redshifts, primary proton angular distributions as well as the account of electron deflection and "additional" processes such as triplet pair production and double pair production could also change the observable SEDs and angular extensions. The effects of these changes are discussed in detail in [62].