EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 685 (May 2024) A&A, 685 (2024) A110 Full HTML
Open Access
Issue
A&A
Volume 685, May 2024
Article Number A110
Number of page(s) 9
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202347809
Published online 15 May 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Blazars are a rare class of active galactic nuclei (AGN) with relativistic jets that are powered by accretion onto a central supermassive black hole (Begelman et al. 1984) and closely aligned with our line of sight (Urry & Padovani 1995). They are the most powerful persistent astrophysical sources of nonthermal electromagnetic radiation in the Universe, and promising candidate sources of other cosmic messengers like high-energy neutrinos (for a recent review, see Murase & Stecker 2023).

The spectral energy distribution (SED) of blazars spans ∼15 decades in energy, starting from radio frequencies and extending to high-energy γ-rays. The flux across the electromagnetic spectrum is variable on different timescales, ranging from minutes to months, with more varying amplitudes found usually in X-rays and γ-rays (for a review, see Böttcher 2019). The blazar SED has a characteristic double-humped shape, with the low-energy component peaking between infrared and X-ray energies and the high-energy component peaking in γ-rays (e.g., Ulrich et al. 1997; Fossati et al. 1998). The low-energy hump is usually explained by electron synchrotron radiation, whereas the origin of the high-energy component is less clear (for a recent review, see Cerruti 2020). According to the most widely accepted scenario for high-energy blazar emission, γ-ray photons are produced via inverse Compton scattering between relativistic electrons and their own synchrotron photons (synchrotron self-Compton, SSC, see e.g., Maraschi et al. 1992; Bloom & Marscher 1996; Mastichiadis & Kirk 1997) or low-energy external radiation fields (external Compton, see e.g., Dermer et al. 1992; Sikora et al. 1994; Ghisellini & Madau 1996).

The likely association of at least some of the high-energy neutrinos detected by IceCube with AGN (e.g., IceCube Collaboration 2018a,b, 2022; Aartsen et al. 2020) has revived interest in so-called hadronic emission models. According to these, the high-energy SED component of a jetted AGN is produced directly via synchrotron radiation of protons (Aharonian 2000; Mücke & Protheroe 2001) or indirectly via a hadronic-initiated electromagnetic cascade (e.g., Mannheim et al. 1991; Mannheim 1993). Detailed SED modeling of the blazar TXS 0506+056, the first astrophysical source associated with high-energy neutrinos in time and space at the ∼3 − 3.5σ level (IceCube Collaboration 2018a,b), has revealed that the electromagnetic imprint of the proton population responsible for the neutrino signal should be subdominant, with the emission of the primary leptonic population producing the observed SED (Keivani et al. 2018; Cerruti et al. 2019; Gao et al. 2019; Petropoulou et al. 2020). These constraints are insensitive to the details of the models, because they basically rely on energy conservation and on the fact that electromagnetic cascades lead to broad spectra that extend in the X-ray and γ-ray ranges (Murase et al. 2018). This new class of hybrid lepto-hadronic models has been applied so far to TXS 0506+056 and a few other blazars that were spatially associated with high-energy neutrinos (Oikonomou et al. 2021; Sahakyan et al. 2023) and they have been used to make neutrino predictions from bright Fermi-detected blazars (Oikonomou et al. 2019; Rodrigues et al. 2024).

Recently, MAGIC Collaboration (2020) reported an unusual spectral feature in the very high-energy (VHE, E > 100 GeV) spectra of the blazar Mrk 501 on the day of the highest X-ray activity above 2 keV measured by Swift/XRT in the last 17 years. This feature at ∼3 TeV was inconsistent (at a 3σ − 4σ confidence level) with the standard analytic functions used to describe VHE spectra (e.g., a power law or a log-parabola). Instead, it was better described by a narrow Gaussian-like component or an exponential log-parabolic component. This rare feature, if true, cannot be explained by the standard one-zone SSC scenario (for alternative leptonic scenarios, see MAGIC Collaboration 2020; Wendel et al. 2021). This raises the question of whether there are conditions under which hadronic-related spectral features can appear at VHE during X-ray flares in the context of the usual SSC scenario for blazar emission.

In this paper, we study the emergence of narrow features in VHE γ-ray spectra of flaring blazars produced via the decay of neutral pions (henceforth, π0 bumps). We focus on high-synchrotron peaked (HSP) blazars (for blazar spectral classes see e.g., Abdo et al. 2010), as these are often detected at VHE energies (Chang et al. 2019). In fact, HSPs are the majority of the ∼90 extragalactic sources detected at VHE by imaging atmospheric Cherenkov telescopes (see TeVCat catalog1). We perform numerical calculations of HSP SEDs during X-ray flares, reaching peak synchrotron energies ∼100 keV, for various model parameters. Our goal is to search for conditions leading to the emergence of π0 bumps in the VHE spectra of flaring blazars and apply our model to Mrk 501.

This paper is structured as follows. In Sect. 2 we outline the model and the numerical code used. In Sect. 3 we present the results of our numerical investigation. We then apply our model to the observations of Mrk 501 in Sect. 4. We discuss the detectability of spectral features related to π0 decay with the Cherenkov Telescope Array (CTA) in Sect. 5, and conclude in Sect. 6 with a discussion of our results. We adopt a flat ΛCDM cosmology with H0 = 70 km s−1 Mpc−1 and Ωm = 0.3.

2. The SSC+π0 model

We modeled the emission region as a spherical blob of radius R′ containing a tangled magnetic field of strength B′. The blob moves with a bulk Lorentz factor, Γ= ( 1− β Γ 2 ) −1/2 $ \Gamma=\left(1-\beta_\Gamma^2 \right)^{-1/2} $, at an angle, θ, with respect to the observer’s line of sight. The Doppler factor was then defined as δ = Γ−1(1 − βΓ cos θ)−1. Henceforth, primes will be used to denote quantities in the frame co-moving with the blob, and unprimed quantities will be used for the observer’s frame.

We assumed that, prior to the onset of an X-ray flare, the blazar SED is fully explained by the SSC radiation of a population of relativistic electrons. These are injected continuously in the emitting region at a constant rate given by Q e ( γ e ′ )= Q e0 γ e ′ − s e H( γ e ′ ; γ e,min ′ , γ e,max ′ ) $ Q_{\rm e}(\gamma^\prime_{\rm e})=Q_{\rm e0} \gamma_{\rm e}^\prime-s_{\rm e} H(\gamma^\prime_{\rm e};\gamma^\prime_{\rm e,min}, \gamma^\prime_{\rm e, max}) $, where H(x; x1, x2) = 1 for x1 ≤ x ≤ x2 and 0 otherwise. The normalization, Qe0, can be expressed in terms of the bolometric injection electron luminosity, L e ′ $ L^\prime_{\rm e} $, as Q e0 = L e ′ /( m e c 2 F( s e , γ e,min ′ , γ e,max ′ )) $ Q_{\rm e0} = L^\prime_{\rm e}/(m_{\rm e} c^2 F(s_{\rm e}, \gamma^\prime_{\rm e, min}, \gamma^\prime_{\rm e, max})) $, where F=ln( γ e,max ′ / γ e,min ′ ) $ F=\ln(\gamma^\prime_{\rm e, max}/\gamma^\prime_{\rm e, min}) $ for se = 2 and F=( γ e,max ′ − s e +2− γ e,min ′ − s e +2)/(2− s e ) $ F=(\gamma_{\rm e, max}^\prime-s_{\rm e}+2-\gamma_{\rm e, min}^\prime -s_{\rm e}+2)/(2-s_{\rm e}) $ otherwise.

For hard electron distributions (se < 2), like those typically used in SSC models of HSP blazars (e.g., Konopelko et al. 2003; Cerruti et al. 2015), the peak synchrotron energy, εs, is determined by the maximum electron Lorentz factor,

γ e , max ′ ≃ 10 5 ( 1 + z ) 1 / 2 ( ε s 3 keV ) 1 / 2 ( B ′ 1 G ) − 1 / 2 ( δ 20 ) − 1 / 2 . $$ \begin{aligned} \gamma ^{\prime }_{\rm e, \max } \simeq 10^5 \ (1+z)^{1/2} \left(\frac{\varepsilon _{\rm s}}{3~\mathrm{keV}} \right)^{1/2} \left( \frac{B^{\prime }}{1\,\mathrm{G}}\right)^{-1/2} \left(\frac{\delta }{20} \right)^{-1/2}. \end{aligned} $$(1)

Therefore, the most straightforward way to simulate a flare where the synchrotron peak shifts to higher energies (henceforth, hard X-ray flares) is to assume that the maximum electron energy also increases, while the other parameters remain constant. The injection electron luminosity could be the same as or higher than its pre-flare value, but for simplicity we adopted the former option. Motivated by observations of hard X-ray flares in Mrk 501 (e.g., Tavecchio et al. 2001), we assumed that the power law slope of the electron distribution does not change.

The same process that accelerates electrons can also accelerate protons to relativistic energies. To minimize the free model parameters, we assumed that the acceleration mechanism at work produces power law distributions for both particle species with sp = se and γ p,max ′ ∼( m e / m p ) γ e,max ′ $ \gamma^\prime_{\rm p,\max} \sim (m_{\rm e}/m_{\rm p}) \, \gamma^\prime_{\rm e, max} $. Similarly to electrons, protons are injected into the radiation region (blob) with a rate given by Q p ( γ p ′ )= Q p0 γ p ′ − s p H( γ p ′ ; γ p,min ′ , γ p,max ′ ) $ Q_{\rm p}(\gamma^\prime_{\rm p})=Q_{\rm p0} \gamma^\prime_{\rm p}{}^{-s_{\rm p}} H(\gamma^\prime_{\rm p};\gamma^\prime_{\rm p,min}, \gamma^\prime_{\rm p, max}) $, where Q p0 = L p ′ /( m p c 2 F( s p , γ p,min ′ , γ p,max ′ )) $ Q_{\rm p0} = L^\prime_{\rm p}/(m_{\rm p} c^2 F(s_{\rm p}, \gamma^\prime_{\rm p, min}, \gamma^\prime_{\rm p, max})) $ and L p ′ $ L^\prime_{\rm p} $ is the bolometric injection proton luminosity.

To produce TeV γ-ray photons (in the observer’s frame) from π0 decay, the parent proton energy in the comoving frame should be

ε p ′ ≈ 10 ε γ δ ( 1 + z ) ≃ 0.5 TeV ( 1 + z ) ( ε γ 1 TeV ) ( 20 δ ) . $$ \begin{aligned} \varepsilon ^{\prime }_{\rm p}\approx 10 \frac{\varepsilon _\gamma }{\delta }(1+z) \simeq 0.5\,\mathrm{TeV} \, (1+z) \left(\frac{\varepsilon _{\gamma }}{1\,\mathrm{TeV}} \right) \left(\frac{20}{\delta }\right). \end{aligned} $$(2)

In contrast to other hadronic scenarios for blazar γ-ray emission (e.g., Cerruti et al. 2015; Petropoulou et al. 2015), the required proton energies are low. The target photons used in photomeson interactions with protons of energy ε p ′ $ \varepsilon^\prime_{\rm p} $ should have comoving energy, ε t ′ ≳ ε ¯ th m p c 2 / 2 ε p ′ $ \varepsilon\prime_{\mathrm{t}}\gtrsim \bar{\varepsilon}_{\mathrm{th}} m_{\mathrm{p}} c^2/ 2 \varepsilon\prime_{\mathrm{p}} $, where ε ¯ th = m π 0 c 2 ( 1 + m π 0 / 2 m p ) ≃ 0.145 $ \bar{\varepsilon}_{\mathrm{th}}=m_{\pi^0} c^2 (1+m_{\pi^0}/2 m_{\mathrm{p}}) \simeq 0.145 $ GeV. The threshold energy condition for interactions with jet photons can then be written as

( ε t 0.5 MeV ) ( ε γ 1 TeV ) ≳ 11.6 ( 1 + z ) − 2 ( δ 20 ) 2 · $$ \begin{aligned} \left(\frac{\varepsilon _{\rm t}}{0.5\,\mathrm{MeV}}\right) \left(\frac{\varepsilon _{\gamma }}{1\,\mathrm{TeV}}\right) \gtrsim 11.6 \, (1+z)^{-2} {\left( \frac{\delta }{20}\right)}^2\cdot \end{aligned} $$(3)

Therefore, the target photons should be very energetic (hard X-rays or soft γ-rays) to enable the production of pions that will decay into TeV γ-rays (in the observer’s frame).

In non-flaring states, the maximum proton energy is even lower than the one needed for the production of pionic TeV γ-rays, namely ε p,max ′ ∼0.05−0.15TeV $ \varepsilon^\prime_{\rm p, max} \sim 0.05{-}0.15\,{\rm TeV} $ for B′∼0.1 − 1 G and δ ∼ 20 (see Eqs. (1) and (2)). During hard X-ray flares with peak synchrotron energies > 100 keV, however, the condition can be met by the most energetic protons in the blob. In this case, the bolometric luminosity of the pionic γ-rays can be estimated as Lπ0 → 2γ ≈ fpπLp, where fpπ is the photopion production efficiency that depends on the number density of synchrotron photons and the Doppler factor (e.g., Murase et al. 2014). Therefore, Lπ0 → 2γ will be a function of L e ′ , L p ′ , B ′ , R ′ $ L^\prime_{\rm e}, L^\prime_{\rm p}, B^\prime, R^\prime $, and δ. The competition between the pionic γ-ray flux and the SSC flux in the relevant energy range, which depends on the aforementioned parameters except L p ′ $ L^\prime_{\rm p} $, will determine the final shape of the VHE spectrum in the SSC+π0 model.

Numerical approach. To study the impact of the π0 component in the VHE spectra of hard X-ray flares, we used the numerical code ATHEνA (Mastichiadis & Kirk 1995; Dimitrakoudis et al. 2012), which computes the multiwavelength photon and all-flavour neutrino spectra as a function of time by solving the kinetic equations for relativistic protons, secondary electrons, and positrons, photons, neutrons, and neutrinos. The physical processes that are included in the code and couple the various particle species are: electron and proton synchrotron emission, synchrotron self-absorption, electron inverse Compton scattering, γγ pair production, photopair (Bethe-Heitler) production, and photomeson production. The latter is modeled based on the results of the Monte Carlo event generator SOPHIA (Mücke et al. 2000), while photopair production is modeled using the Monte Carlo results of Protheroe & Johnson (1996); see also Mastichiadis et al. (2005). For the modeling of other processes we refer the reader to Mastichiadis & Kirk (1995) and Dimitrakoudis et al. (2012). With the adopted numerical scheme, energy is conserved in a self-consistent way, since all the energy gained by one particle species has to come from an equal amount of energy lost by another particle species. The adopted numerical scheme is ideal for studying the development of in-source electromagnetic cascades in various regimes (e.g., Petropoulou 2014; Petropoulou et al. 2019).

We assumed that the non-flaring SED, which is described by a SSC model, has a peak synchrotron frequency, ϵs = 3 keV, which increases to ∼100 keV during the flare. This was achieved by suitably increasing γ e,max ′ $ \gamma^\prime_{\rm e, max} $ and for simplicity we used a step function in time. At the same time we assumed that γ p,max ′ $ \gamma^\prime_{\rm p, max} $ undergoes a similar change, and evolved the system til it reached a new steady state characterized by higher X-ray flux and energy. A steady state was reached within ∼2 light crossing times for all parameter sets we used. Our approach is sufficient for studying the flaring spectra. If, however, one would like to investigate the shape of flares at different wavelengths, a Lorentzian or Gaussian function describing the evolution of the maximum particle energy would be more appropriate (e.g., Moraitis & Mastichiadis 2011).

We computed multi-wavelength flaring spectra for different values of B′, R′, and particle injection luminosities (see Table 1). These were expressed in a dimensionless form through their compactnesses, which were defined as

ℓ e , p = σ T L e , p ′ 4 π R ′ m e , p c 3 . $$ \begin{aligned} \ell _{\rm e,p} = \frac{\sigma _{\rm T} L^{\prime }_{\rm e,p}}{4 \pi R^{\prime } m_{\rm e, p} c^3}. \end{aligned} $$(4)

Table 1.

Input parameters for the simulation of flares with the SSC+π0 model.

In all cases, we used se = sp = 1.7, γ e,min ′ = γ p,min ′ =1 $ \gamma^\prime_{\rm e, min}=\gamma^\prime_{\rm p, min}=1 $, and δ = 20 (unless stated otherwise). For the calculation of the fluxes we considered a blazar at the redshift of Mrk 501, z = 0.034, which corresponds to a luminosity distance, dL = 149.4 Mpc.

3. Results

Figure 1 shows the temporal evolution of the spectrum from a pre-flare steady state (dotted black line) to a new steady state corresponding to a hard X-ray flare (dashed blue line) for different choices of the electron and proton injection compactness. For comparison purposes, we also show the SSC spectrum of the flare (dashed orange line). In both cases, the VHE flare spectrum exceeds the SSC prediction due to the injection of energetic γ-rays from π0 decay. The spectral signature of the latter is clearly visible on the right-hand-side panel of Fig. 1, where the ratio ℓp/ℓe is largest. The π0 bump is also narrow, suggesting that only a small part of the proton energy distribution (here, the highest energy protons) satisfy the pion production threshold with the X-ray synchrotron photons. Even in a less extreme case (in terms of the proton-to-electron luminosity ratio), like the one depicted in the left-hand-side panel, there is a visible hardening of the VHE spectrum that cannot be explained by the SSC emission alone.

thumbnail Fig. 1.

SED snapshots (solid lines) showing the evolution from a pre-flare steady state (dotted black line) to a new steady state that corresponds to the peak of a hard X-ray flare when protons are injected with γ p,max ′ $ \gamma^\prime_{\rm p, max} $ = 103.5 (dashed blue line). For comparison, the pure SSC flare spectrum is overplotted (dotted orange line). The time interval between successive snapshots is 0.1 t cr ′ $ 0.1 \, t^\prime_{\rm cr} $ and the total duration in the observer’s frame is 10 t cr ′ ≈4.6 $ 10 \, t^\prime_{\rm cr}\approx 4.6 $ h. The left and right panels show the results for ℓe = 10−4, ℓp = 1 and ℓe = 10−6, ℓp = 10−1, respectively. Other parameters used are B′ = 0.86 G, R′ = 1015 cm, se = sp = 1.7, γ e,max ′ $ \gamma^\prime_{\rm e, max} $ = 105 (pre-flare) and 106.5 (flare), γ p,max ′ ∼( m e / m p ) γ e,max ′ $ \gamma^\prime_{\rm p,\max} \sim (m_{\rm e}/ m_{\rm p})\gamma^\prime_{\rm e, \max} $, and δ = 20. The attenuation of VHE γ-rays due to the extragalactic background light is not included here.

The effects of the electron injection compactness, ℓe, on the VHE spectra are better illustrated in Fig. 2. Here, we plot spectra of hard X-ray flares (solid lines) computed for three values of ℓe, while keeping all other parameters fixed. For comparison reasons, the pure SSC spectrum of the flare is also shown (dotted lines). As ℓe increases, we recover well-known scaling relations of the synchrotron and Compton fluxes with electron luminosity (e.g., Bloom & Marscher 1996). For a fixed proton luminosity, the π0 bump becomes less prominent as the electron luminosity increases. This might seem counterintuitive at first, but it can be understood as follows. On the one hand, the number density of target photons for pion production increases almost linearly with the electron luminosity, resulting in Lπ0 → 2γ ∝ ℓe. On the other hand, the SSC flux scales roughly quadratically with ℓe, and thus leads to L π 0 → 2 γ / L SSC ∝ ℓ e − 1 $ L_{\pi^0\rightarrow 2 \gamma}/L_{\mathrm{SSC}} \propto \ell_{\mathrm{e}}^{-1} $.

thumbnail Fig. 2.

Non-flaring and X-ray flaring SEDs computed for three values of the electron compactness, ℓe, as is indicated in the plot and ℓp = 1. All other parameters are the same as in Fig. 1.

The dependence of this ratio on the electron compactness is better illustrated in Fig. 3, which summarizes the findings from a suite of numerical calculations performed for different parameter values (see Table 1). Here, we plot the ratio of the VHE γ-ray luminosities (integrated above 1 TeV) for the SSC+π0 and SSC models as a function of ℓe. For a given proton compactness, as ℓe increases, the VHE spectrum becomes leptonically dominated. For sufficiently high ℓe values, the luminosity ratio should asymptotically approach unity. The contribution of the pionic γ-rays becomes evident for lower values of the electron compactness where ( L SSC + π 0 / L SSC ) > 1 TeV ∝ ℓ e − 1 $ \left(L_{\mathrm{SSC+\pi^0}}/L_{\mathrm{SSC}} \right)_{ > 1 \rm TeV} \propto \ell_{\mathrm{e}}^{-1} $. In this regime, the luminosity ratio is proportional to ℓp, as is shown by comparing symbols of the same color. Changes in the radius of the emission region have a lesser impact on the luminosity ratio (for all other parameters fixed), as is evidenced by comparing markers of the same shape and different colors.

thumbnail Fig. 3.

Ratio of VHE luminosity (> 1 TeV) due to SSC and π0 decay and SSC alone as a function of the electron compactness for various parameters (see inset legend).

The findings presented in this section suggest that synchrotron-dominated flares with peak energies ϵs ≳ 100 keV might be ideal periods for the search of π0−decay bumps in the TeV spectra of HSP blazars. In the next section we discuss the applicability of our model to the flaring spectra of Mrk 501 in 2014 (MAGIC Collaboration 2020).

4. Application to Mrk 501

Mrk 501 is a HSP blazar at z = 0.034 that is well known for its X-ray and VHE flaring activity. It is also characterized as a transient extreme blazar. An example of such behavior was the 1997 flare, during which the peak synchrotron energy exceeded 100 keV, marking a change of more than two orders of magnitude in the peak energy compared to that of the non-flaring state (e.g., Paltani et al. 1998; Tavecchio et al. 2001). However, Mrk 501 has also shown the characteristics of an extreme blazar during non-flaring activity, as was observed throughout the extensive multi-instrument campaign in the year 2012 (Ahnen et al. 2018), indicating that being an extreme blazar may not be a permanent state, but rather a state that changes over time.

MAGIC Collaboration (2020) reported the results of a multi-wavelength campaign performed in July 2014. For a period of about two weeks, Mrk 501 exhibited the highest X-ray activity ever detected by Swift during its operation. VHE flaring activity correlated with the X-rays was also recorded. Moreover, a narrow feature at ∼3 TeV was detected at a significance level of about 4σ in the VHE spectrum measured with the MAGIC telescopes on July 19, 2014 (MJD 56857.98). This was also the day with the highest X-ray flux (> 0.3 keV) measured with Swift. We note that the spectral feature at ∼3 TeV may also have been present the day before and the day after (and hence may have lasted three days, from July 18 until July 20); but the detection of the spectral feature in these two days is not significant, perhaps due to the much shorter MAGIC observation time (half an hour, instead of the 1.5 h from the observation that started on July 19). This feature is inconsistent with the empirical functions used to describe VHE spectra, namely a power law, a log-parabola, or a log-parabola with an exponential cutoff. As a result, the typical one-zone SSC model fails to explain the VHE spectra (MAGIC Collaboration 2020).

Assuming that this feature was real, we applied the model presented in the previous sections to the multi-wavelength data of Mrk 501 on MJD 56857.98. Figure 4 shows the respective single-night SED adopted from MAGIC Collaboration (2020). Most of the data samples were selected from observations as contemporaneous as possible, taken within three hours of each other. We note however that the Fermi-LAT flux points were computed for 4-day and 10-day bins centered on MJD 56857.98. The reason for the longer integration time, in comparison to the X-ray and VHE gamma-ray data, is the limited sensitivity of Fermi-LAT in detecting Mrk 501 over short timescales. Moreover, the hard X-ray flux points (green circles) represent the integrated flux in the 14–195 keV energy range that was derived using two methods (for details see MAGIC Collaboration 2020).

thumbnail Fig. 4.

Left panel: broadband SED of Mrk 501 on MJD 56857.98 compiled using data from MAGIC (orange circles, corrected for EBL attenuation), Fermi-LAT (averaged over four and ten days, blue circles), Swift-XRT (red triangles), Swift-BAT (integrated flux in 14–195 keV from one-day observations, shown with green circles), Swift-UVOT (purple circles), KVA (purple square), and Metsähovi (pink square). The SSC+π0 model SED is overplotted (solid line). For comparison, we also show the SED when Bethe-Heitler pair production is neglected (dashed line), and the purely leptonic (SSC) model (dotted line). Right panel: zoom-in of the γ-ray spectrum.

Figure 4 shows a SSC+π0 model SED of the flare (solid gray line), with the SSC contribution indicated by the dotted line (for the parameters used, see Table 2). The VHE spectrum can be explained by the superposition of SSC emission and the π0 bump. It should be noted that the sharpness of the spectral feature is limited by the energy resolution of the numerical code (i.e., five grid points per decade in photon energy). While for this model we used a practically monoenergetic proton distribution (see Table 2), it is also possible to explain the bump with a hard power law proton distribution, with sp ∼ 1.5 − 1.7. In this case, however, one would need to appropriately increase the proton compactness to compensate for the “inactive” lower-energy protons of the distribution that would not contribute to pion production, because they would not satisfy the energy threshold for pion production with the hard X-ray photons of the flare. The fact that a monoenergetic proton distribution is not a necessary condition for the appearance of the pion bump is also supported by the numerical results presented in Fig. 1, where an extended power law distribution was used. The emergence of the pion bump is instead related to the ratio of ℓp/ℓe, as is shown in Fig. 2.

Table 2.

Parameter values used in the one-zone SSC+π0 SED models of Mrk 501 presented in Fig. 4.

For a better comparison with the BAT data, we also indicate the flux of the model integrated in the same energy range (asterisk), which is consistent with the BAT value. For a more detailed spectral comparison (on a single night) between the model and the data in hard X-rays, more sensitive hard X-ray observations would be needed (e.g., observations with NuSTAR). In addition to the π0 bump, the synchro-Compton emission of secondary pairs from proton-photon pair production (Bethe-Heitler) is visible in the low-energy part of the spectrum (below ∼1 eV) and makes a contribution to the LAT energy range too. Because of the hard (< 2) power law slope of the primary electron distribution and the fact that the secondary emission dominates at energies ≲1 eV, the minimum primary electron Lorentz factor is a nuisance model parameter on which only an upper limit can be placed.

In conclusion, our results motivate an empirical spectral model that consists of two log-parabolic components or a log-parabola and an exponential cutoff power law component to describe the composite VHE spectrum, as is detailed in the next section.

5. VHE simulated spectra for CTA

We performed simulations of VHE spectra using the open-source software GammaPy (Deil et al. 2017; Aguasca-Cabot et al. 2023) to test the detectability of π0 bumps with the CTA (CTA, Cherenkov Telescope Array Consortium 2019). For our simulations and fitting, we followed the steps below. First, we defined the source position in the sky (using the sky coordinates of Mrk 501). Secondly, we assigned the SSC+π0 model shown in Fig. 4 to the source, and derived its theoretical photon flux (in units of s−1 cm−2 TeV−1), after multiplying the theoretical spectrum with e−τEBL(z, E), where the optical depth, τEBL, was adopted from Franceschini et al. (2008). Then, we simulated the observed γ-ray spectrum using the instrument response functions (IRF) for the northern site of CTA (currently including four large-sized telescopes and nine medium-sized telescopes), for a 30-min exposure time. In all of the simulations we used a 0.11 degree extraction region centered on the source location, which for a 30-min exposure yields ∼3850 source counts and ∼130 background counts between 0.1 TeV and 10 TeV.

Motivated by the theoretical spectra shown in Fig. 4, we used two empirical spectral models to fit the simulated data, namely a power law with an exponential cutoff (ECPL),

ϕ ECPL ( E ) = ϕ 0 · ( E E 0 ) − Γ exp ( − λ E ) , $$ \begin{aligned} \phi _{\rm ECPL}(E) = \phi _0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma } \exp (- \lambda E) ,\end{aligned} $$(5)

and a dual-component model (ECPL+LP) defined as

ϕ ( E ) = ϕ ECPL ( E ) + ϕ LP ( E ) . $$ \begin{aligned} \phi (E) = \phi _{\rm ECPL}(E) + \phi _{\rm LP}(E). \end{aligned} $$(6)

The second function of the composite model is a log-parabola,

ϕ LP ( E ) = ϕ 0 ( E E 0 ) − α − β log ( E E 0 ) , $$ \begin{aligned} \phi _{\rm LP}(E) = \phi _0 \left( \frac{E}{E_0} \right) ^ { - \alpha - \beta \log { \left( \frac{E}{E_0} \right) } }, \end{aligned} $$(7)

which has a peak energy (in a E2ϕLP(E) versus E plot) given by

E peak = E 0 e ( − α + 2 ) / ( 2 β ) . $$ \begin{aligned} E_{\rm peak} =E_0 \ e^{(-\alpha +2)/(2\beta )}. \end{aligned} $$(8)

In both models, the energies are measured in TeV units, and the flux normalization, ϕ0, has units of cm−2 s−1 TeV−1. The pivot energy, E0, in both models was frozen to 1 TeV during the fitting.

For the spectral analysis, we fit each model to the simulated γ-ray spectrum using cash statistics (Cash 1979), following standardized procedures built within GammaPy2. The fitting was performed via log-likelihood functions, so for each model a value for the fit statistics (−2log L) was returned that allowed us to compare the goodness of models.

To explore the significance of any detection based on the SSC+π0 model, we performed three tests. First, we simulated a single spectrum based on the theoretical model and fit it with the single-component (ECPL) and dual-component (ECPL+LP) models (see left panel in Fig. 5). In addition to the improvement in the fit statistics (∼43), the ECPL+LP model also improves the residuals of the fit, defined as the (model-data)/error (see middle and bottom panels of left-hand side Fig. 5).

thumbnail Fig. 5.

Left panel: simulated γ-ray spectrum of Mrk 501 for CTA based on the SSC+π0 model shown in Fig. 4 and a 30-min exposure. For illustrative purposes, the decomposition of the best-fit dual-component (ECPL+LP) model is also shown (colored histograms). The middle and bottom panels show the residuals of the fit, defined as (model-data)/error, for the ECPL and ECLP+LP models, respectively. Right panel: unfolded binned γ-ray spectrum (blue markers) and best-fit dual-component model (solid orange line), with the residuals shown in the bottom panel. The ECPL and LP components of the best-fit spectrum are also plotted for clarity. For comparison, the observed MAGIC spectrum (obtained with a ∼1.5 h exposure) is overplotted (gray squares).

As a second test, we repeated our simulations to establish the statistical significance of our results. We used the Akaike information criterion (AIC, Akaike 1974) for the model comparison. For each model we computed the AIC number by adding to the fit statistics a penalty of 2K, where K is the number of model parameters. To provide a statistical estimate, we repeated the procedure for 10 000 simulated spectra, and computed the difference in AIC values between two models (i.e., ΔAIC = AICModel 1 − AICModel 2) to determine which model was more likely (see discussion by Burnham & Anderson 2004). We found that for a 30-min exposure with CTA the computed SSC+π0 model has similar AIC values to the ECPL+LP model, while the ECPL model alone yields a worse fit, with a median ΔAIC ∼ −34 (see Fig. 6). This indicates that the ECPL+LP is preferred3 by more than 5σ (based on the median ΔAIC values). Moreover, 99.8% of simulated spectra favor the ECPL+LP model with a significance of more than 99%. Only for a small number of simulated spectra (21 out of 10 000) do we find that ΔAIC > −9 (see orange-colored histogram). We also repeated the simulations for 15-min and 60-min exposures and found that the ΔAIC values follow an almost linear trend with the exposure time. In particular, the dual-component model is favored, with a median ΔAIC of ∼16 and ∼66 over the ECPL model for 15-min and 60-min exposures, respectively. Moreover, considering the ΔAIC distribution for the 60-min exposure simulations, the π0 bump would be detected at a median significance of 7.5σ.

thumbnail Fig. 6.

Histograms of the difference in AIC values computed as AICModel 1 − AICModel 2 when (i) Model 1 is the ECPL+LP model and Model 2 is the simulated model (blue), and (ii) Model 1 is the ECPL+LP model and Model 2 is ECPL model (orange). Vertical lines indicate the median values of ΔAIC. More negative values suggest that Model 1 is more probable than Model 2.

As a final test we simulated 10 000 spectra based on a featureless broadband spectrum, using the best fit parameters of the ECPL model from the previous test. We then fit all simulated spectra with the ECPL model or the ECPL+LP model. As was expected, the ΔAIC values between these two models are clustered around zero, with a standard deviation of less than three, while less than 1% of samples have |ΔAIC|> 1. This indicates that there is no statistical preference for the more complex spectral model.

The log-parabola of the dual-component model in this case does not typically describe a narrow feature, but is rather a broad component that contributes to the continuum spectrum. This is better illustrated in Fig. 7, where we show the fitted parameters of the ECPL+LP model on the spectra that were simulated with the featureless ECPL model (light gray color). For comparison, we overplot the parameters obtained from fitting the spectra that were simulated based on the SSC+π0 model (maroon color). We find that the ECPL parameters are consistent between the two sets of simulations. However, the LP parameters derived from the fit to the SSC+π0 simulations are much more clustered around a central value, corresponding to a peak energy of ≃3 TeV, which is consistent with the SSC+π0 model. On the contrary, the gray-colored contours are more dispersed, as they do not represent a sharp feature like the proposed pion bump but are a result of over-fitting the data with two components.

thumbnail Fig. 7.

Corner plot of the parameters obtained from the fitted models to the simulated CTA spectrum. Maroon-colored contours and distributions are based on simulations with the SSC+π0 model fitted with the dual-component (ECPL+LP) model. The gray samples correspond to simulations based on the featureless ECPL model (i.e., no “bump”) that are also fitted with the dual-component model. Contour lines indicate the 11.8%, 39.3%, 67.5%, and 86.4% levels. The parameter Epeak denotes the peak energy of the LP component (in TeV). This is not a free model parameter, but was obtained from α and β – see Eq. (8).

6. Summary and discussion

In this paper we have investigated the production of narrow spectral features in the VHE spectra of flaring blazars through π0 decay. During X-ray flares reaching peak synchrotron energies of ∼100 keV, the number density of hard X-ray photons increases, thus making pion-production feasible for protons of a few TeV in energy. What ultimately determines the appearance of narrow π0 bumps in the VHE spectra is the competition between the π0γ-ray flux and the SSC flux in the relevant energy range. We have shown that flares in HSP blazars with peak energies ≳100 keV and small values of the Compton ratio can be ideal periods for the emergence of π0 bumps at TeV energies.

We applied the SSC+π0 model to multi-wavelength data of Mrk 501 that were obtained during a period of high X-ray activity that lasted two weeks. Within this 14-day period, a narrow spectral feature at ∼3 TeV was detected with MAGIC on MJD 56857.98 at a ∼3σ − 4σ confidence level, and it could have been present for a maximum of three consecutive nights (MAGIC Collaboration 2020). Here, we focused on single-night data of MJD 56857.98, and searched for a steady-state SSC+π0 model that describes the broadband spectrum well. We showed that the VHE γ-ray spectrum can be explained by a superposition of γ-ray emission by neutral pion decay, inverse Compton scattering emission of secondary leptons produced via proton-photon pair production, and SSC emission of primary electrons. For a good description of the narrow spectral feature, the proton distribution has to be energetically dominated by protons with energy Ep ∼ 10 TeV(ϵγ/1 TeV). This could be realized by either a mono-energetic proton distribution or by a hard (sp < 2) power law distribution extending to lower energies, as was discussed in Sect. 4.

Given that the effective photopion production cross section is ∼10−3σT and that the number density of the X-ray target photons is low, a very high proton luminosity is needed to produce a pionic bump visible in the TeV spectrum. Let us consider, for example, protons with Lorentz factor γ p ′ ∼ 10 3.2 $ \gamma^\prime_{\rm p} \sim 10^{3.2} $. These protons produce pions at a threshold with photons of comoving energy ∼90 keV, whose number density is n X ′ ≃ 10 3 $ n^\prime_X \simeq 10^3 $ cm−3 for the model spectrum shown in Fig. 4. The photopion production efficiency can be estimated as f p γ = t esc ′ / t p γ ′ ≈ R ′ n X ′ σ ̂ p γ ≃ 10 − 9 $ f_{\mathrm{p\gamma}} = t\prime_{\mathrm{esc}}/t\prime_{\mathrm{p}\gamma} \approx R\prime n\prime_X {\hat{\sigma}}_{\mathrm{p}\gamma} \simeq 10^{-9} $, where tpγ is the proton energy loss timescale due to the photopion production process and σ ̂ p γ ≃ 7 × 10 − 29 cm 2 $ \hat{\sigma}_{\mathrm{p\gamma}} \simeq 7 \times 10^{-29}~\rm cm^2 $ is the cross section of the process weighted by the inelasticity (Dermer & Menon 2009). The amount of energy transferred from relativistic protons to secondaries per unit time (in the observer’s frame) can then be estimated as L pγ ≈ f pγ δ 4 L p ′ ≃ 10 44 erg s -1   δ 1 4 L p,49 ′ $ L_{\rm p\gamma} \approx f_{\rm p\gamma} \delta^4 L^\prime_{\rm p} \simeq 10^{44}\,{\rm erg\, s^{-1}}\ \delta_1^4 L^\prime_{\rm p,49} $, where we introduced the notation qm = 10m (in cgs units). Therefore, very high proton luminosities are required to compensate for the extremely low pion-production efficiency in the SSC+π0 model.

Furthermore, the flaring region that produces the TeV pion bump is particle dominated, with u p ′ ≫ u e ′ ≳ u B ′ $ u^\prime_{\rm p} \gg u^\prime_{\rm e} \gtrsim u^\prime_B $. For the parameter values listed in Table 2, the relativistic proton number density is n p ′ ∼ u p ′ /( γ p ′ m p c 2 )∼ l p /(3 σ T R ′ γ p ′ )≃3× 10 4 $ n^\prime_{\rm p} \sim u^\prime_{\rm p}/(\gamma^\prime_{\rm p} m_{\rm p} c^2) \sim \ell_{\rm p} / (3 \sigma_{\rm T} R^\prime \gamma^\prime_{\rm p}) \simeq 3\times 10^4 $ cm−3 for γ p ′ ∼ 10 3.2 $ \gamma^\prime_{\rm p} \sim 10^{3.2} $. The total number density in relativistic protons would be even higher if, instead of an effectively mono-energetic proton distribution, an extended power law starting from γ p,min ′ $ \gamma^\prime_{\rm p, min} $ = 1 were adopted. For instance, we find similar VHE gamma-ray spectra to those displayed in Fig. 4 for n p,PL ′ / n p ′ ∼166 $ n^\prime_{\rm p, PL}/n^\prime_{\rm p} \sim 166 $, where the subscript “PL” indicates an extended power law distribution with index sp = 1.5. The average jet power may be determined by the conditions during non-flaring states, assuming that pionic bumps, like the feature identified in Mrk 501, are rare. If almost all cold protons in the jet region producing the flare become relativistic with γ p ′ ∼ 10 3 $ \gamma^\prime_{\rm p} \sim 10^3 $ during the flare, we can set a lower limit to the kinetic power of the jet in non-flaring states, which is attributed to the cold protons, namely L p , j ( c ) > π R 2 Γ j 2 β j c n p , c ′ m p c 2 ≃ 1.8 × 10 47 $ L_{\mathrm{p,j}}^{(c)} > \pi R^2 \Gamma_j^2 \beta_j c n\prime_{\mathrm{p, c}} m_{\mathrm{p}} c^2 \simeq 1.8 \times 10^{47} $ erg s−1 for n p,c ′ = n p ′ $ n^\prime_{\rm p, c} = n^\prime_{\rm p} $ (or 3.6 × 1049 erg s−1 for n p,c ′ = n p,PL ′ $ n^\prime_{\rm p, c} = n^\prime_{\rm p, PL} $), Γj ∼ δ = 13, and βj = 1. The lower limit on the jet power estimated for non-flaring conditions can readily exceed the Eddington luminosity of the central black hole in Mrk 501, for a typical mass of MBH = 109 M⊙ (inferred by the MBH − σ* relation, Barth et al. 2002). Therefore, the SSC+π0 model for TeV flares requires a proton-dominated flaring region, protons reaching energies not much higher than a few tens of TeV, and super-Eddington average jet powers. These features, except for the lower maximal proton energies, are similar to those found in other leptohadronic models for blazar emission (e.g., Petropoulou et al. 2015; Zdziarski & Bottcher 2015).

In Petropoulou et al. (2023) we considered a different scenario where relativistic protons interact with an external source of X-ray radiation (e.g., from a radiatively inefficient accretion disk), which is not directly observed, to produce the pion bump. In this case, the regions producing the X-ray flare observed with Swift and the VHE emission are decoupled and are described by different physical parameters. In particular, the “SSC zone” has a weaker magnetic field and is more extended than the “pion zone”, which is assumed to be located closer to the base of the jet (for details, see Table 1 in Petropoulou et al. 2023). The former zone can fully account for the optical-X-ray emission of the flare and produce a fraction of the γ-ray flux, while the “pion zone” can dominate the emission at GeV-TeV energies during the flare (see Fig. 2 in Petropoulou et al. 2023). In this scenario, the appearance of the π0 bump during a hard X-ray flare, as was observed in Mrk 501, should be coincidental. Flares from the SSC zone or pion zone would also result in different strengths of correlation between X-ray and TeV γ-ray fluxes. However, this model would require the presence of a persistent “pion zone” to explain a weak X-ray-TeV correlation during non-flaring epochs as was observed, which would also imply a very high (time-average) jet power.

A distinct spectral feature like the π0 bump could be detected by future CTA observations, as we demonstrated in Sect. 5. In particular, we performed simulations based on the theoretical SSC+π0 model shown in Fig. 4 using the “Alpha Configuration” of the CTA Northern Array. We showed that for events collected within a 30-min CTA observation of Mrk 501, spectral fitting would require a multi-component empirical model, and a feature like the 3 TeV bump could be detected at a 5σ significance level.

Our analysis opens up some interesting possibilities for the interpretation of multi-wavelength blazar spectra during flaring states. We demonstrated that a harder VHE γ-ray spectrum than the usual SSC component or, more occasionally, a distinct neutral pion bump at VHE energies during hard X-ray flares, can be suggestive of a relativistic hadronic component in the jet that otherwise would remain hidden. Further testing of this possibility requires simultaneous and detailed spectral observations in the soft and hard X-rays and TeV γ-rays that can be performed by the current generation of X-ray instruments (e.g., NuSTAR) and Cherenkov Telescopes (namely, H.E.S.S., MAGIC, and VERITAS) or, even better, by more sensitive future instruments such as CTA and X-ray satellites with broadband capabilities, like HEX-P (Madsen et al. 2019).

1

http://tevcat.uchicago.edu

2

See documentation and examples: https://docs.gammapy.org

3

The likelihood of a model over another can be estimated as exp(−ΔAIC/2).

Acknowledgments

The authors would like to thank the anonymous referee for a constructive report. MP acknowledges support from the MERAC Fondation through the project THRILL and from the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “2nd call for H.F.R.I. Research Projects to support Faculty members and Researchers” through the project UNTRAPHOB (Project ID 3013). GV acknowledges support by Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “3rd Call for H.F.R.I. Research Projects to support Postdoctoral Researchers” through the project ASTRAPE (Project ID 7802). JBG acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) through the Spanish State Research Agency, under Severo Ochoa Program 2020-2023 (CEX2019-000920-S). DP acknowledges support from the Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) under Germany’s Excellence Strategy EXC-2094–390783311.

References

  1. Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020, Phys. Rev. Lett., 124, 051103 [Google Scholar]
  2. Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010, ApJ, 716, 30 [NASA ADS] [CrossRef] [Google Scholar]
  3. Aguasca-Cabot, A., Donath, A., Feijen, K., et al. 2023, https://doi.org/10.5281/zenodo.7734804 [Google Scholar]
  4. Aharonian, F. A. 2000, New Astron., 5, 377 [NASA ADS] [CrossRef] [Google Scholar]
  5. Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2018, A&A, 620, A181 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Akaike, H. 1974, IEEE Trans. Autom. Control, 19, 716 [Google Scholar]
  7. Barth, A. J., Ho, L. C., & Sargent, W. L. W. 2002, ApJ, 566, L13 [NASA ADS] [CrossRef] [Google Scholar]
  8. Begelman, M. C., Blandford, R. D., & Rees, M. J. 1984, Rev. Mod. Phys., 56, 255 [Google Scholar]
  9. Bloom, S. D., & Marscher, A. P. 1996, ApJ, 461, 657 [NASA ADS] [CrossRef] [Google Scholar]
  10. Böttcher, M. 2019, Galaxies, 7, 20 [Google Scholar]
  11. Burnham, K. P., & Anderson, D. R. 2004, Sociol. Methods Res., 33, 261 [Google Scholar]
  12. Cash, W. 1979, ApJ, 228, 939 [Google Scholar]
  13. Cerruti, M. 2020, Galaxies, 8, 72 [NASA ADS] [CrossRef] [Google Scholar]
  14. Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, MNRAS, 448, 910 [Google Scholar]
  15. Cerruti, M., Zech, A., Boisson, C., et al. 2019, MNRAS, 483, L12 [NASA ADS] [CrossRef] [Google Scholar]
  16. Chang, Y. L., Arsioli, B., Giommi, P., Padovani, P., & Brandt, C. H. 2019, A&A, 632, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Cherenkov Telescope Array Consortium (Acharya, B. S., et al.) 2019, Science with the Cherenkov Telescope Array (World Scientific Publishing Co Pte. Ltd.) [Google Scholar]
  18. Deil, C., Zanin, R., Lefaucheur, J., et al. 2017, in 35th International Cosmic Ray Conference (ICRC2017), 301, 766 [Google Scholar]
  19. Dermer, C. D., & Menon, G. 2009, High Energy Radiation fromBlack Holes: Gamma Rays, Cosmic Rays, and Neutrinos (Princeton Univerisity Press) [CrossRef] [Google Scholar]
  20. Dermer, C. D., Schlickeiser, R., & Mastichiadis, A. 1992, A&A, 256, L27 [NASA ADS] [Google Scholar]
  21. Dimitrakoudis, S., Mastichiadis, A., Protheroe, R. J., & Reimer, A. 2012, A&A, 546, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini, G. 1998, MNRAS, 299, 433 [Google Scholar]
  23. Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Gao, S., Fedynitch, A., Winter, W., & Pohl, M. 2019, Nat. Astron., 3, 88 [NASA ADS] [CrossRef] [Google Scholar]
  25. Ghisellini, G., & Madau, P. 1996, MNRAS, 280, 67 [CrossRef] [Google Scholar]
  26. IceCube Collaboration 2018a, Science, 361, 147 [NASA ADS] [CrossRef] [Google Scholar]
  27. IceCube Collaboration (Aartsen, M. G., et al.) 2018b, Science, 361, eaat1378 [NASA ADS] [Google Scholar]
  28. IceCube Collaboration (Abbasi, R., et al.) 2022, Science, 378, 538 [CrossRef] [PubMed] [Google Scholar]
  29. Keivani, A., Murase, K., Petropoulou, M., et al. 2018, ApJ, 864, 84 [Google Scholar]
  30. Konopelko, A., Mastichiadis, A., Kirk, J., de Jager, O. C., & Stecker, F. W. 2003, ApJ, 597, 851 [NASA ADS] [CrossRef] [Google Scholar]
  31. Madsen, K., Hickox, R., Bachetti, M., et al. 2019, BAAS, 51, 166 [NASA ADS] [Google Scholar]
  32. MAGIC Collaboration (Acciari, V. A., et al.) 2020, A&A, 637, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Mannheim, K. 1993, A&A, 269, 67 [NASA ADS] [Google Scholar]
  34. Mannheim, K., Biermann, P. L., & Kruells, W. M. 1991, A&A, 251, 723 [NASA ADS] [Google Scholar]
  35. Maraschi, L., Ghisellini, G., & Celotti, A. 1992, ApJ, 397, L5 [CrossRef] [Google Scholar]
  36. Mastichiadis, A., & Kirk, J. G. 1995, A&A, 295, 613 [NASA ADS] [Google Scholar]
  37. Mastichiadis, A., & Kirk, J. G. 1997, A&A, 320, 19 [NASA ADS] [Google Scholar]
  38. Mastichiadis, A., Protheroe, R. J., & Kirk, J. G. 2005, A&A, 433, 765 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Moraitis, K., & Mastichiadis, A. 2011, A&A, 525, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Mücke, A., & Protheroe, R. J. 2001, Astropart. Phys., 15, 121 [Google Scholar]
  41. Mücke, A., Engel, R., Rachen, J. P., Protheroe, R. J., & Stanev, T. 2000, Comput. Phys. Commun., 124, 290 [Google Scholar]
  42. Murase, K., & Stecker, F. W. 2023, High-Energy Neutrinos from Active Galactic Nuclei (World Scientific Publishing Co., Pte. Ltd.), 483 [Google Scholar]
  43. Murase, K., Inoue, Y., & Dermer, C. D. 2014, Phys. Rev. D, 90, 023007 [Google Scholar]
  44. Murase, K., Oikonomou, F., & Petropoulou, M. 2018, ApJ, 865, 124 [Google Scholar]
  45. Oikonomou, F., Murase, K., Padovani, P., Resconi, E., & Mészáros, P. 2019, MNRAS, 489, 4347 [Google Scholar]
  46. Oikonomou, F., Petropoulou, M., Murase, K., et al. 2021, JCAP, 2021, 082 [CrossRef] [Google Scholar]
  47. Paltani, S., Courvoisier, T. J. L., & Walter, R. 1998, A&A, 340, 47 [NASA ADS] [Google Scholar]
  48. Petropoulou, M. 2014, MNRAS, 442, 3026 [CrossRef] [Google Scholar]
  49. Petropoulou, M., Dimitrakoudis, S., Padovani, P., Mastichiadis, A., & Resconi, E. 2015, MNRAS, 448, 2412 [Google Scholar]
  50. Petropoulou, M., Yuan, Y., Chen, A. Y., & Mastichiadis, A. 2019, ApJ, 883, 66 [NASA ADS] [CrossRef] [Google Scholar]
  51. Petropoulou, M., Murase, K., Santander, M., et al. 2020, ApJ, 891, 115 [Google Scholar]
  52. Petropoulou, M., Mastichiadis, A., Paneque, D., & Gonzalez, J. B. 2023, in Proceedings of 7th Heidelberg International Symposium on High-Energy Gamma-Ray Astronomy - PoS(Gamma2022), 417, 081 [CrossRef] [Google Scholar]
  53. Protheroe, R. J., & Johnson, P. A. 1996, Astropart. Phys., 4, 253 [NASA ADS] [CrossRef] [Google Scholar]
  54. Rodrigues, X., Paliya, V. S., Garrappa, S., et al. 2024, A&A, 681, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Sahakyan, N., Giommi, P., Padovani, P., et al. 2023, MNRAS, 519, 1396 [Google Scholar]
  56. Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153 [Google Scholar]
  57. Tavecchio, F., Maraschi, L., Pian, E., et al. 2001, ApJ, 554, 725 [NASA ADS] [CrossRef] [Google Scholar]
  58. Ulrich, M.-H., Maraschi, L., & Urry, C. M. 1997, ARA&A, 35, 445 [NASA ADS] [CrossRef] [Google Scholar]
  59. Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 [NASA ADS] [CrossRef] [Google Scholar]
  60. Wendel, C., Becerra González, J., Paneque, D., & Mannheim, K. 2021, A&A, 646, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Zdziarski, A. A., & Bottcher, M. 2015, MNRAS, 450, L21 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Input parameters for the simulation of flares with the SSC+π0 model.

In the text
Table 2.

Parameter values used in the one-zone SSC+π0 SED models of Mrk 501 presented in Fig. 4.

In the text

All Figures

thumbnail Fig. 1.

SED snapshots (solid lines) showing the evolution from a pre-flare steady state (dotted black line) to a new steady state that corresponds to the peak of a hard X-ray flare when protons are injected with γ p,max ′ $ \gamma^\prime_{\rm p, max} $ = 103.5 (dashed blue line). For comparison, the pure SSC flare spectrum is overplotted (dotted orange line). The time interval between successive snapshots is 0.1 t cr ′ $ 0.1 \, t^\prime_{\rm cr} $ and the total duration in the observer’s frame is 10 t cr ′ ≈4.6 $ 10 \, t^\prime_{\rm cr}\approx 4.6 $ h. The left and right panels show the results for ℓe = 10−4, ℓp = 1 and ℓe = 10−6, ℓp = 10−1, respectively. Other parameters used are B′ = 0.86 G, R′ = 1015 cm, se = sp = 1.7, γ e,max ′ $ \gamma^\prime_{\rm e, max} $ = 105 (pre-flare) and 106.5 (flare), γ p,max ′ ∼( m e / m p ) γ e,max ′ $ \gamma^\prime_{\rm p,\max} \sim (m_{\rm e}/ m_{\rm p})\gamma^\prime_{\rm e, \max} $, and δ = 20. The attenuation of VHE γ-rays due to the extragalactic background light is not included here.
In the text
thumbnail Fig. 2.

Non-flaring and X-ray flaring SEDs computed for three values of the electron compactness, ℓe, as is indicated in the plot and ℓp = 1. All other parameters are the same as in Fig. 1.
In the text
thumbnail Fig. 3.

Ratio of VHE luminosity (> 1 TeV) due to SSC and π0 decay and SSC alone as a function of the electron compactness for various parameters (see inset legend).
In the text
thumbnail Fig. 4.

Left panel: broadband SED of Mrk 501 on MJD 56857.98 compiled using data from MAGIC (orange circles, corrected for EBL attenuation), Fermi-LAT (averaged over four and ten days, blue circles), Swift-XRT (red triangles), Swift-BAT (integrated flux in 14–195 keV from one-day observations, shown with green circles), Swift-UVOT (purple circles), KVA (purple square), and Metsähovi (pink square). The SSC+π0 model SED is overplotted (solid line). For comparison, we also show the SED when Bethe-Heitler pair production is neglected (dashed line), and the purely leptonic (SSC) model (dotted line). Right panel: zoom-in of the γ-ray spectrum.
In the text
thumbnail Fig. 5.

Left panel: simulated γ-ray spectrum of Mrk 501 for CTA based on the SSC+π0 model shown in Fig. 4 and a 30-min exposure. For illustrative purposes, the decomposition of the best-fit dual-component (ECPL+LP) model is also shown (colored histograms). The middle and bottom panels show the residuals of the fit, defined as (model-data)/error, for the ECPL and ECLP+LP models, respectively. Right panel: unfolded binned γ-ray spectrum (blue markers) and best-fit dual-component model (solid orange line), with the residuals shown in the bottom panel. The ECPL and LP components of the best-fit spectrum are also plotted for clarity. For comparison, the observed MAGIC spectrum (obtained with a ∼1.5 h exposure) is overplotted (gray squares).
In the text
thumbnail Fig. 6.

Histograms of the difference in AIC values computed as AICModel 1 − AICModel 2 when (i) Model 1 is the ECPL+LP model and Model 2 is the simulated model (blue), and (ii) Model 1 is the ECPL+LP model and Model 2 is ECPL model (orange). Vertical lines indicate the median values of ΔAIC. More negative values suggest that Model 1 is more probable than Model 2.
In the text
thumbnail Fig. 7.

Corner plot of the parameters obtained from the fitted models to the simulated CTA spectrum. Maroon-colored contours and distributions are based on simulations with the SSC+π0 model fitted with the dual-component (ECPL+LP) model. The gray samples correspond to simulations based on the featureless ECPL model (i.e., no “bump”) that are also fitted with the dual-component model. Contour lines indicate the 11.8%, 39.3%, 67.5%, and 86.4% levels. The parameter Epeak denotes the peak energy of the LP component (in TeV). This is not a free model parameter, but was obtained from α and β – see Eq. (8).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.