Nonthermal Radiation of the Extreme TeV Blazar 1ES 0229+200 from Electromagnetic Cascades on Infrared Photon Field

Extreme TeV blazars (ETBs) are active galactic nuclei with jets presumably pointing towards the observer having their intrinsic spectral energy distributions (SEDs) peaked at an energy in excess of 1 TeV. These sources typically reveal relatively weak and slow variability as well as an extremely high frequency of the low-energy SED peak compared to other classes of blazars. It proved to be exceedingly hard to incorporate all these peculiar properties of ETBs into the framework of a reasonable $\gamma$-ray emission model. ETB physics have recently attracted great attention in the astrophysical community, underlying the importance of the development of self-consistent ETB emission model(s). We propose a new scenario for the formation of X-ray and $\gamma$-ray spectra of ETBs assuming that electromagnetic cascades develop in the infrared photon field surrounding the central blazar engine. This scenario does not invoke compact fast-moving sources of radiation (so-called"blobs"), in agreement with the apparent absence of fast and strong variability of ETBs. For the case of the extreme TeV blazar 1ES 0229+200 we propose a specific emission model in the framework of the considered scenario. We demonstrate that this model allows to obtain a good fit to the measured SED of 1ES 0229+200.


Introduction
Blazars are active galactic nuclei with their jets pointing in the close direction to the observer's line of sight. The observable spectral energy distribution (SED= E 2 dN/dE, where E is the observable photon energy and N is the number of photon counts) of a typical blazar reveals two prominent peaks. As a rule, the lower-energy SED peak occupies the frequency range spanning from infrared to X-rays; the higher-energy SED peak falls into the γ-ray domain. Strong and fast variability, sometimes on timescales as short as minutes [1,2], is typical for most classes of blazars. This flaring emission is usually interpreted in a theoretical framework assuming particle acceleration and subsequent γ-ray production inside relativistic plasmoids [3] (usually called "blobs" [4,5]) propagating along the jets.
In 2006, the H.E.S.S. collaboration reported the discovery of γ-ray emission from the blazar 1ES 1101-232 at redshift z = 0.186. The SED of this source was measured up to the energy of about 3 TeV [6]. At such energies, γ rays are subject to strong absorption on extragalactic background light (EBL) photons due to the γγ → e + e − pair production (PP) process [7,8] 1 . Neglecting the effects arising from intergalactic electromagnetic cascades 2 , the intrinsic spectrum of 1ES 1101-232 (i.e. the spectrum of γ rays escaping into the intergalactic medium) was reconstructed. Remarkably, this intrinsic spectrum does not reveal a cutoff even at the highest observable energy (for instance, see Fig. 17 (top-left) of [9]).
Such blazars with intrinsic SEDs peaked at an energy E ph > 1 TeV are called "extreme TeV blazars" (ETBs). ETBs have recently came into the spotlight of the topical research in astrophysics [10]. By the end of 2019 about ten ETBs were discovered [10]. However, it proved to be notoriously difficult to arrive at a satisfactory theoretical understanding of these sources. In the present paper we inquire into the nature of high energy (E > 100 MeV) γ rays and X-rays (E = 0.3 − 100 keV) from extreme TeV blazars. In Section 2 we summarize the basic peculiar properties of ETBs which any model seeking to explain their emission must be able to reproduce. In Section 3 we briefly discuss the difficulties of the existing models.
In the rest of the paper we propose and discuss an emission model for the blazar 1ES 0229+200 (z = 0.14 [11]). This source is convenient for the aims of the present study because its broadband spectrum was well measured from the infrared to the very high energy (VHE, E > 100 GeV) domain and, moreover, because the spectrum was made publicly available in Supplementary Information of [10]. In Section 4 we examine the low-energy part of the spectrum and speculate about the possible spectrum of thermal radiation inside 1ES 0229+200. In Section 5 we fit the X-ray and γ-ray part of the spectrum with semi-analytic templates for electromagnetic cascade spectra, including both the inverse Compton and the synchrotron components. This paper is meant to be the first in series of works aimed at the understanding the nature of extreme TeV blazars. Therefore, we leave some subjects for future study (see Section 6). Finally, we conclude in Section 7.

Peculiar properties of extreme TeV blazars
Below we list some observational properties of extreme TeV blazars. ETB emission model building could be guided by the requirement that a successful model must reproduce these properties.

1.
By definition, the intrinsic SED of any ETB is peaked at an energy E ph > 1 TeV (see Section 1). In particular, the reconstructed intrinsic spectrum of the blazar 1ES 0229+200 [12] does not show an evidence for a break or cutoff even at the highest accessible energy of 10 TeV [9]. 2.
ETBs do not reveal fast (day-scale or even week-scale) or strong (flux change by an order of magnitude or more) γ-ray variability 3 .

3.
At the same time, relatively weak and slow (if compared to other classes of blazars) γ-ray variability of ETBs was indeed detected. In particular, the VERITAS Collaboration finds that the blazar 1ES 0229+200 is variable on yearly timescales [12]. In the present paper we strive to explain the variability quasi-period as small as several months, at the same time explaining the absence of faster variability. 4.
The energy of the low-energy SED peak E pl for ETBs is usually relatively high compared to other classes of blazars and falls into the X-ray domain [10]; this could incur significant difficulties in some models as discussed in Section 3.

Models of extreme TeV blazar emission and their difficulties
Below we list some models that were proposed to describe the multiwavelength emission of ETBs. We note that most of these models assume the existence of blobs inside the jets of ETBs, at odds with the absence of fast and strong variability of these sources. As well, many models require that the magnetic field energy density (u B ) is well below the energy density of particles (u p ): The K u ≫ 1 condition is very hard to meet. Indeed, in this case the accelerating particles quickly escape from the blob, stopping acceleration immediately (for a discussion of the K u ∼ 1 asumption see e.g. [14] and [15], p.7). In addition, such systems are usually characterised by fast magnetic field amplification driven by the particle current on a timescale comparable to the particle acceleration timescale [16]. This process would lead to a shift of the low-energy SED component towards higher energies, impairing the fit to the observed spectrum.

1.
In synchrotron self-Compton (SSC) models γ rays are produced as a result of inverse Compton (IC, e − γ → e − ′ γ ′ or e + γ → e + ′ γ ′ ) scattering of accelerated electrons on synchrotron photons radiated by the same electron population. For the case of ETBs (assuming γ-ray production in blobs) these models require K u ≫ 1 as well as high values of minimal Lorentz factor of the primary electrons (γ e−min > 10 4 ) (e.g. [12]). The last condition requires specific acceleration mechanism such as the Blandford-Znajek mechanism [17] that operates only in the immediate vicinity of the event horizon of the central black hole in the blazar. This mechanism does not chime well with the concept of particle acceleration inside a fast-moving blob that quickly escapes the central engine.

2.
A leptonic model presented in [18] assumes that a high value of γ e−min is due to a specific preacceleration process, namely, the transfer of energy from protons to electrons. However, it is not clear why the same process does not operate in jets of other (non-extreme) blazars. This model, as well, implies K u ≫ 1.

3.
Another leptonic model assuming the IC process on cosmic microwave background (CMB) photons was proposed in [19]. It still suffers from the γ e−min problem outlined above. In addition, the CMB energy density typically dominates very far from the central black hole, at the multi-parsec or even kiloparsec scale. Thus appears the problem of ultrarelativistic electron transport to such great distances without appreciable energy losses and/or the problem of accelerating such electrons very far from the central engine.

5.
The proton sychrotron model (e.g. [20]) requires the acceleration of primary protons up to ∼ 10 20 eV. The maximum energy of synchrotron photon produced by the protons is E sp−max ∼ 100 GeV, and for electrons co-accelerated with the protons E se−max ∼ 50 MeV [21]. To produce hard (spectral index < 2) observable spectrum up to the energy of E ph = 10 TeV, the bulk motion of a blob or a jet corresponding to the Doppler factor of D > E ph /E sp−max = 100 is required. In this case, the typical observable energy of synchrotron photons emitted by electrons would be ∼ E se−max D ≈ 5 GeV, and not ∼1-10 keV typical for extreme blazars. Therefore, it is difficult to account for the low-energy SED component self-consistently.
For these reasons, a radically different approach to the ones proposed in the abovediscussed works is justified when trying to understand the emission of extreme TeV blazars. As a first step of this understanding process, let us consider the nature of photon field on which the observed γ rays could be produced.

The low-energy part of the spectrum
The low-energy part of the observed SED for the blazar 1ES 0229+200 is plotted in Figure 1. All data are taken from tables provided as Supplementary Material for [10]. Xray data are from the Neil Gehrels Swift Observatory [22], X-ray Telescope (Swift-XRT) and the NuSTAR Observatory [23], according to the analysis of [24] (see red circles in their Figure 2). The data in the energy range of 0.05-10 eV are from the the Swift Ultra-violet Optical Telescope (Swift-UVOT) (for two epochs, see Supplementary Information of [10] and [24] for more details), the Wide-field Infrared Survey Explorer (WISE) satellite [25], the ASI Space Science Data Center (ASI-SSDC) (again, see [24] for more details), and [26] (see legend in Figure 1). Below the hydrogen ionisation edge, the SED consists of three distinct bumps. The SED of the central bump is well fitted with the giant elliptical galaxy template of [27] (hereafter S98); indeed, the host galaxies for most ETBs are believed to be giant ellipticals [24]. However, in the ultraviolet (UV) and infrared (IR) ranges some residuals with respect to this template are evident; these residuals are well fitted with two blackbody components that have the temperatures of 175 K and 2.1 · 10 4 K and the luminosities of 4 · 10 43 erg/s and 1.2 · 10 44 erg/s, respectively.
The ultraviolet component may represent a stellar phenomenon [28]. In principle, some contribution to the UV component from an accretion flow is also possible. 4 It is believed that the inner part of the accretion flow in Bl Lacs is radiatively inefficient, with high temperature beyond the UV energy range. However, as discussed in [30], the outer part of the accretion flow could represent a standard Shakura-Sunyaev accretion disk [31]. Now let us discuss the infrared part of the spectrum in more details (see Figure 2). The template of S98 predicts a bump peaked at the energy of ∼ 10 −2 eV 5 , while observations reveal a more energetic component peaked at ≈ 6 · 10 −2 eV. The discrepancy in the peak energy and intensity of the IR component between the template of S98 and the observations may be in part due to the underestimation of the typical temperature of the dust and/or of the dust content in the star-forming regions.
However it is also possible that a part of the observed IR photons was produced near the central engine, at the spatial scale between ∼ 0.01 pc and 100 pc. 6 This IR photon field could serve as a target for high-energy electrons that could produce VHE γ rays via the IC process. The same electrons could produce the X-ray component of the spectrum via the synchrotron process. Let us see what shape of the observable spectrum we could expect 4 this was the case of the elliptical galaxy NGC 4552 [29] 5 this component is probably due to radiation of dust 6 the intensity of photon fields in the energy range ≤ 5 · 10 −2 eV in 1ES 0229+200 is poorly constrained with available observations in the framework of this scenario. For simplicity we assume that the IR photon field is isotropic.

The measured SED
The high-energy part of the observed SED for the blazar 1ES 0229+200, together with the X-ray component, is presented in Figure 3. The spectra obtained with imaging atmospheric Cherenkov telescopes were taken from [12] (VERITAS, red circles) and [32] (H.E.S.S., red triangles).
We have performed a Fermi-LAT [33] data analysis for 1ES 0229+200 over the time period of 2008-08-04 -2020-09-20 and the energy range of 1 GeV -500 GeV following the standard recommendations for a point-like source spectral analysis issued by the Fermi-LAT collaboration, using the 4FGL catalog [34], the galactic background model gll_iem_v07, and the isotropic background model iso_P8R3_SOURCE_V2_v1. The Fermitools software (version 1.2.23) 7 as well as the fermiPy package [35] (version 0.19.0) 8 were utilized for this analysis. For details concerning the H.E.S.S. and VERITAS data analysis the reader is referred to [32] and [12], respectively; for details concerning all datasets except the H.E.S.S., VERITAS, and Fermi-LAT ones -to [10]. We note that the observations with different instruments were not strictly simultaneous. However, as the source under study does not display a very strong or fast variability, in what follows we use the datasets resulting from these observations. The Fermi-LAT SED of 1ES 0229+200 obtained by us is shown in Figure 3 as red squares. The last and the third from the last energy bins yielded only upper limits which are not shown in the figure.

Qualitative considerations
Primary electrons could be accelerated near the central black hole of the blazar 9 via the Blandford-Znajek process [17] up to at least 100 TeV (e.g. [36]). IC γ rays produced by these electrons could create electron-positrons pairs which, in turn, generate secondary (cascade) γ rays. If the resulting electromagnetic cascade has more than two generations, its spectrum is usually almost independent on both the type of the primary particle (electron/γ ray) and its energy. This property is called "cascade universality" [37].
Under the condition of universality, the spectrum of IC γ rays has a definite shape: it is well-described with a dN/dE ∝ E −1.5 power law at relatively low energies, below a certain energy E x ; with a dN/dE ∝ E −2.0 dependence above E x ; and above the γ-ray horizon energy E a this spectrum reveals a cutoff caused by the PP process. Below we adopt an analytic representation for the IC component based on a smoothly broken power law functional form. For more details, the reader is referred to [37].
The same equations as those presented in [37] for the IC component could be written for the synchrotron component, except that the PP process does not matter for the synchrotron component due to much lower energies of synchrotron photons compared to IC photons. Indeed, both the synchrotron photon energy and the IC γ ray energy in the Thomson regime are proportional to the square of the primary electron energy (e.g. [38]). For the case of photon field assumed below, the Thomson regime is valid up to the energy of ∼ 1 TeV. The deviation from the Thomson regime will be accounted for with the correction factor f KN . 7 https://github.com/fermi-lat/Fermitools-conda 8 https://github.com/fermiPy/fermipy 9 typically at the distance less that several gravitational radii if counted from the event horizon

The model SED
Assuming the cascade universality, we calculate the spectrum of cascade γ rays and synchrotron photons as follows: where s stands for the synchrotron component, c -for the IC component; K i are the normalization factors for each component, γ 1 = 1.5, γ 2 = 2.0 [37]; E xi are the energies of the break of the power-law index value, ε = 5 is the parameter defining the sharpness of the break, E mi are the maximal energies of produced photons and γ rays, and τ int (E) is the value of the intrinsic optical depth for γ rays 10 . K i are determined from fitting the measured SED. For the synchrotron component we estimate E xs and E ms as follows (see e.g. [39]): (2) 10 we note that for the synchrotron component τ int = 0 where B = 0.9 mG is the characteristic magnetic field strength, and E me = 100 TeV is the primary electron spectrum cutoff energy. E xe = E a /2, where E a is the energy of the intrinsic γ-ray horizon: and E b = 6 · 10 −2 eV is the characteristic energy of the IR photons. We set E t = 1 TeV. The intrinsic optical depth for γ rays (more precisely, its dependence on the energy τ int (E)) is usually calculated by integrating the interaction rate over the distance. The definition of the interaction rate could be found e.g. in [40]. The energy of the γ-ray horizon is the energy corresponding to τ int (E) = 1. E t is the effective threshold energy of a background photon interacting with a γ-ray via the PP process. For instance, the threshold for a 1 TeV γ-ray is ∼ 1 eV; E t is inversely proportional to the γ-ray energy.
For the IC component: where γ xe = E xe /m e is the electron Lorentz factor and m e is the electron mass [eV]. Using the approximation of [41] for the secondary γ-ray spectrum resulting from the IC process for the case of thermal photon field, we obtain f KN (E xe ) ≈ 0.15. We use a simple parametrization τ int (E) = E/E a . In this work we are not interested in the precise shape of the intrinsic γγ absorption cutoff, because the intrinsic γ-ray horizon energy appears to be significantly higher than the EBL γ-ray horizon energy. After thus calculating the spectum of X-rays and γ rays resulting from cascades, we apply the effects of redshift and EBL absorption assuming the EBL model of [42], and plot the resulting model curves in Figure 3 (solid green curve for the synchrotron component and solid black curve for the IC component). This fit was obtained without any formal optimization with a dedicated algorithm, but just by choosing reasonable values of input parameters and performing several tries. We note that the blue curve corresponds to the case of E me = 100 TeV.
Finally, we try an option of E me = 300 TeV (dashed green curve for the synchrotron component). The observable spectrum of the IC component is not sensitive to the value of E me because E me significantly exceeds the EBL γ-ray horizon energy.
For demonstration purposes, in Figure 4 we plot the model γ-ray spectrum multiplied by E 1.5 neglecting the absorption effect on the EBL. Furthermore, green dashed curve in this figure does not include the effect of internal absorption either. Vertical black dashed lines correspond to the values of E xc (left line) and E a (right line). Remarkably, the ratio E a /E xc is below ten. In more standard blazar models such a γ-ray spectrum would imply a very narrow spectrum of parent electrons, as was noted in [43]. In our model, however, this feature appears self-consistently.

Basic and auxiliary parameters
The proposed model of the spectral shape for the blazar 1ES 0229+200 has four basic parameters listed in Table 1. Besides these, the fifth parameter could be introduced in order to ensure proper absolute normalization to the SED measurements. The first three parameters appeared in the previous Subsection. We note that the physical meaning of B is the strength of the turbulent component of magnetic field. The regular magnetic field along which the particles move could have a significantly higher strength than B. The fourth parameter defines the relative normalization of the IC and synchrotron components as follows: where dN c /dE and dN s /dE are the spectra of the IC and synchrotron component, respectively. For the specific case considered in this paper, in the energy regions covered with observational data these spectra are almost independent from the E me value (see the previous Subsection). For these reasons, the shape of the model broadband spectrum is mostly defined by only three parameters: E b , B, and K cs .
In the previous Subsection we presented a fit to the measured SED assuming specific values of these parameters. We note that this solution is likely not unique and other sets of parameters could lead to a reasonable fit to the experimental data.
This fit included some auxiliary parameters that had their values fixed, but these values could in principle be obtained from ab initio calculations. Such direct and detailed calculations are presently underway; the corresponding results will be published elsewhere. The above-mentioned parameters include (γ 1 ,γ 2 ) (the values of these follow from the theory developed in [37]), E t (its value is defined by the physics of the PP process), and ε. We note that these parameters are in fact defined by the values of E b , B, and K cs . Concerning the f KN parameter, we have already performed a calculation of its value, as explained in the previous Subsection.

Explanation of peculiar properties of ETBs
As a summary of the above discussion, we present Table 2 which is called to explain, in a concise form, how the proposed model meets the expectations formulated in Section 2. Namely, when the IR photon field is sufficiently soft (E b ≪ 1 eV), it is possible to obtain the energy of the γ-ray SED peak E ph > 1 TeV, because the PP effect does not lead to a strong cutoff in the cascade γ-ray spectrum below several TeV, and the Klein-Nishina cross section suppression [38] at such energies is still moderate. A sufficiently high energy of the primary electrons (E me ≥ 100 TeV) allows to produce multi-TeV cascade electrons which, in turn, radiate multi-keV synchrotron photons. Our model does not invoke relativistic blobs, naturally explaining the absence of strong or fast γ-ray variability. Finally, slow or weak variability could be explained by a variation in the accretion rate which, in turn, affects the properties of the primary electron spectrum (see the next Subsection for more details) and/or leads to variations in photon field properties which affect the spectrum of cascade γ rays.

Property
The physical reason The absence of relativistic blobs Slow/weak var. present Variations of accretion rate/photon field 6.2. The origin of the primary electrons As dissussed in Subsection 5.2, the primary electrons could have been accelerated in the vacuum gap near the central black hole via the Blandford-Znajek process [17,36,[44][45][46]. The magnetic field in the vacuum gap B gap could be of order of 1 kG; therefore, the acceleration process could be accompanied by intense curvature radiation. For B gap > 1 kG, the typical energy of the curvature photon radiated by ∼100 TeV electron falls in the gap between the X-ray and γ-ray parts of the spectrum of the blazar 1ES 0229+200 [47]. Future hard X-ray and MeV γ-ray observatories could possibly detect a distinct component from curvature radiation in the spectrum of 1ES 0229+200.
Another possible origin of ∼100 TeV primary electrons is the Bethe-Heitler pair production process [48,49] on the IR photon field. In this case, the relevant primary proton energy is ∼100 PeV -1 EeV. These protons could also have been accelerated via the Blandford-Znajek process [36,50,51].
A variation of the accretion rate induces a change of the magnetic field strength in the vacuum gap. The maximum total luminosity of the gap strongly depends on B gap : L gap ∝ B 2 gap [36,51]. Therefore, the change in the accretion rate could induce slow variability on the timescale of months or years. Finally, we note that a faster variability is not excluded altogether: discharges in the black hole magnetospheres are in principle possible leading to effects somewhat similar to lightning [52].
Above, we had neglected the intergalactic cascade contribution to the observable γray spectrum. Indeed, if B EGMF > 10 fG and λ EGMF >100 kpc, then this contribution from strongly beamed sources could be safely neglected [63,64]. Otherwise, the intergalactic cascade component could yield an observable signal. In 2017-2021, a series of studies was conducted in our research group discussing a possible contribution of the intergalactic cascade component to the observable spectra of extreme TeV blazars [9,65,66]. Preliminary estimates show that the proposed model allows a significant contribution of cascade γ rays to the observable SED at low energies (E < 1 TeV). A detailed model of the observable spectrum of the blazar 1ES 0229+200 for the option of B EGMF ≪ 10 fG is currently in development and will be presented elsewhere.

Additional remarks on the model
Some γ-rays are inevitably produced via the SSC mechanism. However, the ratio of the number density for the IR photons to the X-ray photons is K IR−X ∼ 10 3 /(θ 2 jet ), where θ jet is the angular radius of the jet. 11 For θ jet = 0.1 K IR−X ∼ 10 5 ; therefore, the SSC component is subdominant if both photon fields occupy the same volume. 12 Assuming that the UV component of photon field represents a stellar phenomenon on ∼kpc scale (see Section 4), the pair production depth on this UV excess is small (≪1). In the present paper we do not invoke bulk motion of matter in blobs/jet to explain the non-thermal radiation spectrum of 1ES 0229+200. In this respect, our model is similar to the one presented in [36] for the case of the radiogalaxy M87.
The blackbody fits of the UV and IR components presented in Figure 1 are formal ones and do not have direct physical meaning. In particular, the IR component could be in part produced as the result of synchrotron radiation of electrons with relatively low energy that have escaped the jet.

Conclusions
In the present work we proposed a new scenario for the formation of X-ray and and γ-ray spectrum of the blazar 1ES 0229+200 assuming the acceleration of primary electrons up to the energy of ∼100 TeV near the central black hole of the blazar with the subsequent development of electromagnetic cascades on the infrared photon fields surrounding the central engine. This scenario explains all known peculiar properties of extreme TeV blazars. Remarkably, the spectrum of 1ES 0229+200 is very similar to the spectrum of electromagnetic cascade in the universal regime. Moreover, we presented a specific model that fits the measured SED of 1ES 0229+200 well. A somewhat modified version of this model could be applied to the case of other ETBs. Detailed calculations to this end are currently underway and will be published elsewhere.  Data Availability Statement: The Fermi-LAT spectrum of the blazar 1ES0229+200 resulting from our analysis will be shared on reasonable request to the corresponding author.