On the narrow spectral feature at $\sim$3 TeV in the MAGIC spectrum of Mrk 501

The broadband spectral energy distributions (SEDs) of blazars can be characterized by a double-peak structure extending from radio frequencies to the $\gamma$-ray energies (Abdo et al., 2010a; Giommi et al., 2012; Acciari et al., 2021b). The low-energy peak located in the infrared to X-ray range is interpreted as synchrotron emission from non-thermal relativistic electrons accelerated in jets. Based on the position of the synchrotron peak ($\nu_{\rm pk}$), blazars are divided into high-synchrotron-peaked (HSP; $\nu_{\rm pk}>10^{15}$ Hz), intermediate-synchrotron-peaked (ISP; $10^{14}<\nu_{\rm pk}<10^{15}$ Hz), and low-synchrotron-peaked (LSP; $\nu_{\rm pk}<10^{14}$ Hz) blazars (Abdo et al., 2010b). The high-energy peak, located at $\gamma$-ray energies, can be well explained by leptonic models. $\gamma$-ray photons from HSPs are usually generated via inverse-Compton (IC) scattering the synchrotron photons by the same electron population that produces the synchrotron emission (SSC) (Maraschi, Ghisellini & Celotti, 1992; Yan, Zeng & Zhang, 2014).

Blazars exhibit fast and large-amplitude variability across the entire electromagnetic spectrum (e.g., Hayashida et al., 2012; Paliya, Sahayanathan & Stalin, 2015; Bartoli et al., 2016; Yan et al., 2018). However, the physical origin of the variability and the associated electron acceleration mechanism are currently not well understood. Diffusive shock acceleration (e.g., Drury, 1983; Kirk, Rieger & Mastichiadis, 1998; Summerlin & Baring, 2012) and stochastic acceleration (e.g., Dermer, Miller & Li, 1996; Katarzyński et al., 2006; Yan, Zeng & Zhang, 2012) are widely considered in blazar jets.

Mrk 501 is a famous TeV HSP (Quinn et al., 1996), which is intensively monitored by various astronomical detectors, from radio frequencies to TeV $\gamma$ rays (e.g., Abdo et al., 2011; Neronov, Semikoz & Taylor, 2012; Shukla et al., 2015; Aleksić et al., 2015; Ahnen et al., 2018). Very recently, an extreme X-ray activity in July 2014 was reported by Acciari et al. (2020). During the flaring activity, the X-ray and TeV emissions have been found to be correlated, and the amplitude of the TeV variability is twice as large as that of the X-ray variability. This indicates that the leptonic model is most likely favored. Interestingly, a narrow feature at $\sim 3$ TeV was observed by MAGIC stereoscopic telescope on 2014 July 19 (MJD 56857.98), and the significance of this peak-feature is estimated to be $\sim 4\sigma$ (Acciari et al., 2020). This kind of spectral feature is rare, which will provide insight into the underlying emission mechanism in the jet. Acciari et al. (2020) proposed several explanations for this narrow spectral feature, including pileup in the electron energy distribution due to stochastic acceleration and IC pair cascade induced by electrons accelerated in a magnetospheric vacuum gap. Wendel et al. (2021) provided details for IC pair cascade scenario.

Interestingly, Acciari et al. (2021a) reported that the X-ray flux measured by Swift-BAT in the 15-50 keV band is significantly higher than that expected from the simple extrapolation of the Swift-XRT spectral data for HSP Mrk 421 during MJD 57422-57429 (2016 February 4-11). This suggests an additional emission component beyond the single synchrotron emission (Acciari et al., 2021a). This BAT excess may be related to the presence of the additional (and narrow) component in the TeV spectrum of Mrk 501 (Acciari et al., 2021a). Actually, the results of NuSTAR data analysis have revealed possible excess in hard X-ray emission for Mrk 421 (Kataoka & Stawarz, 2016) and PKS 2155-304 (Madejski et al., 2016). These results may imply that the appearance of the occasional TeV narrow spectral feature is not unique in Mrk 501. The future improved TeV $\gamma$-ray detectors, such as Cherenkov Telescope Array (CTA; e.g., Actis et al., 2011; Acharya et al., 2013; Sol et al., 2013), would detect such a narrow sharp feature with high significance.

Our goal is to interpret the narrow feature at $\sim 3$ TeV of Mrk 501 during the extreme X-ray activity in 2014 July 19 in the framework of a one-zone leptonic model. This paper is structured as follows: we describe our model in Section 2 and results in Section 3, followed by the discussion and conclusions in Section 4.

Our model assumes a single, homogenous emission region of radial size $R^{\prime}$, which moves along the jet with bulk Lorentz factor $\Gamma$. This region is assumed to filled with random magnetic field with the strength of $B^{\prime}$. The accelerated electrons in the region will radiate photons through synchrotron and inverse-Compton (IC) scattering processes. For an observer located at a small viewing angle $\theta_{\rm obs}$ with respect to the jet, the radiative luminosity in the comoving frame is boosted by a factor of $\delta_{\rm D}^{4}$, and the variability timescale is reduced by a factor of $\delta_{\rm D}^{-1}$ (e.g., Urry & Padovani, 1995). Here $\delta_{\rm D}$ is the Doppler factor, which is given by $\delta_{\rm D}=2\Gamma/(1+N_{\Gamma}^{2})$, where $N_{\Gamma}\equiv\Gamma\theta_{\rm obs}$ (Dermer et al., 2014). The emitting region size is constrained by the causality condition: $R^{\prime}\lesssim ct_{\rm var}\delta_{\rm D}/(1+z)$, where $t_{\rm var}$ is the measured variability timescale and $z$ is the blazar redshift. It is usually assumed $N_{\Gamma}=1$, so that $\delta_{\rm D}=\Gamma$ (e.g., Nalewajko, Begelman & Sikora, 2014).

In the model, we take into account the synchrotron self-absorption, as well as the absorption for $\gamma-$ray photons caused by the interaction with the internal synchrotron radiation field. The numerical method to calculate such emission spectra is the same as that in Dermer & Menon (2009). Throughout the paper, quantities in the comoving frame of the emission region are primed, and those in the stationary frame or observer’s frame are unprimed.

In the model, we do not distinguish the acceleration region from the emission region, and we assume that the distribution of relativistic emitting electrons is the result of the evolution of the injected electrons undergoing acceleration, energy losses and escape. To reduce the number of model parameters, the injection of electrons is assumed to be mono-energetic, which is given by the $\delta-$function

where $V^{\prime}=4\pi R^{\prime 3}/3$ is the comoving volume of the emission region, $m_{\rm e}$ is the rest mass of an electron, $L_{\rm inj,i}$, and $\gamma_{\rm inj,i}^{\prime}$ denote the luminosity and the Lorentz factor of the injected electrons respectively, with the start and finish times denoted by $T_{\rm s,i}$ and $T_{\rm f,i}$, respectively. Here, $H(x;a,b)$ denotes the Heaviside function defined by $H=1$ if $a\leq x\leq b$, and $H=0$ everywhere else.

In this section, we explain the formation of the narrow spectral feature observed at $\sim\$3 TeV. We adopt a two-injection scenario.

In our model, the size of the source ($R^{\prime}$), the characteristic acceleration time ($t^{\prime}_{\rm acc}$) and the luminosity of the first injection electrons ($L_{\rm inj,1}^{\prime}$) could be derived from the relationships in Section 2 with the observed variability timescale ($t_{\rm var}$), the synchrotron peak frequency ($\nu_{\rm pk,1}$), and the total flux of the synchrotron emission ($F_{\rm syn,1}$). Therefore, we can describe the emission with the model parameters: $B^{\prime},~{}\delta_{\rm D},~{}a,~{}\beta_{\rm A},~{}t_{\rm var}$, $\nu_{\rm pk,1},~{}F_{\rm syn,1}$ ,$\gamma^{\prime}_{\rm inj,1},~{}T^{\prime}_{\rm s,1},~{}T^{\prime}_{\rm f,1}$, $\gamma_{\rm inj,2}^{\prime},~{}L^{\prime}_{\rm inj,2}$, $T^{\prime}_{\rm s,2}$ and $T^{\prime}_{\rm f,2}$.

The shortest variability found in the MWL data sample is on the order of one day (Acciari et al., 2020), we therefore use $t_{\rm var}=1$ day. Since we consider a very inefficient electron escape from the blob, we fix $\beta_{\rm A}$ to 0.03, which is compatible with the quasi-linear theory (e.g., Becker et al., 2006). We fix $\gamma_{\rm inj,1}^{\prime}=\gamma_{\rm inj,2}^{\prime}=600$. The values of $\gamma_{\rm inj,1}^{\prime}$ and $\gamma_{\rm inj,2}^{\prime}$ have negligible effect on the SSC emission. We set $T^{\prime}_{\rm s,1}=0$ and $T^{\prime}_{\rm f,1}=\Delta t^{\prime}=0.02R^{\prime}/c$ for the impulsive injection. When the first injection electrons are accelerated to be close to half of the equilibrium energy, we start the second injection. In our test, we find that $T^{\prime}_{\rm s,2}$ is not sensitive in modeling the SED. The time interval between the two injections can be estimated by Equation 24. The remaining parameters are free.

Fig. 3 shows SSC modeling to the observed SED of Mrk 501 during 2014 July 19 (MJD 56857.98). The values of the parameters are given in Table 1 (Model A), in which we list all input parameters and the physical quantities associated with our model. As shown in the figure, our modeling can reproduce the observed SED well when we select the evolution time $t^{\prime}_{\rm evo}=T^{\prime}_{\rm f,2}=2.74t^{\prime}_{\rm dyn}$. We note that the optical-UV flux produced by the theoretical model is lower than the observed data. However, it may be consistent with the variability behavior of MWL observations, which showed that the variability in the optical-UV bands is much lower than that in X-ray and TeV energies (Acciari et al., 2020). This indicates that the optical-UV emissions have another origin.

We separate contributions of the different segments of the electron spectrum IC scattering different segments of the synchrotron photons. The results indicate that the narrow TeV bump at $\sim$ 3 TeV is the contribution from the pileup bump of electrons IC scattering synchrotron photons produced by the power-law branch of electrons. Moreover, the high flux observed by Swift-BAT above 10 keV is interpreted as the synchrotron emission of the pileup bump of electrons which is resulting from the first injection, while the X-rays below $\sim 1$ keV is dominated by the contribution of the power-law branch of electrons resulting from the second injection. The peak frequency of the synchrotron emission produced by the power-law branch is located at $\nu_{\rm pk,2}\simeq 5\times 10^{16}$ Hz. With the given values of $B^{\prime}$ and $\delta_{\rm D}$ reported in Table 1, one can obtain $\gamma_{\rm pk,2}^{\prime}\simeq 1.7\times 10^{5}$, which is in good agreement with that obtained in the slow cooling model proposed by Acciari et al. (2020). It should be noted that our results are obtained by using a self-consistent physical model. Moreover, it is interesting to note that the shock acceleration rather than the stochastic acceleration is required to be dominant, i.e., a¿1.

To illustrate the impact of duration of the first injection, we also display the result of the SED modeling with $T^{\prime}_{\rm f,1}=T^{\prime}_{\rm s,2}$ in Fig. 4 with the parameters values listed in Table 1 (Model B). We have not found a significant difference between Model A and Model B. The pileup structure in the EED can be formed in both continuous and impulsive injections.

In this work, the emitting electron distribution is determined by numerically solving a Fokker-Planck equation that self-consistently incorporates both shock and stochastic acceleration processes, radiative losses, as well as synchrotron self-absorption and internal $\gamma\gamma$ absorption.

In the model, the parameters $a$, $B^{\prime}$, $\delta_{\rm D}$ and $T^{\prime}_{\rm f,2}$ are four key physical quantities to reproduce the observed SED. In the following, we examine the impact of the four parameters on SED modeling.

In Fig. 5, we show the modeling results with varying $a$. The narrow feature is not visible when the stochastic acceleration dominates over the shock acceleration, i.e., $a<1$. The necessary condition for the formation of the pile-up in electron distribution is that the electrons are well confined within the emission zone (i.e., an inefficient electron escape). Additionally, two injection episodes of electrons are required to reproduce the SED of Mrk 501. Therefore, the narrow feature in the TeV energies disappears naturally if the electron escape effect will be important (i.e., efficient escape scenario), and/or the first population of electrons can be neglected (i.e., a single-injection scenario). In those cases, the broadband SEDs reported in Acciari et al. (2020) without the sharp TeV feature should be reproduced well following the prescription in Section 2.2.2.

A constraint on the Doppler factor $\delta_{D}$ can be obtained from transparency of $\gamma-$rays due to pair production absorption (e.g., Zdziarski & Lightman, 1985; Finke, Dermer, & Böttcher, 2008; Yan, Zeng, & Zhang, 2012). Due to internal photon-photon annihilation, the optical depth for a $\gamma$-ray photon with observed energy $E_{\gamma}$ can be estimated as (e.g., Abdo et al., 2011)

where $F_{0}$ is the observed monochromatic flux energy density at the observed photon energy

Using $F_{0}\simeq 3\times 10^{-11}~{}\rm erg/cm^{2}/s$ at the observed energy $E_{0}\simeq 11.3$ eV, one has $\tau_{\gamma\gamma}(3\rm TeV)\simeq 0.16$ for the given value of $\delta_{\rm D}=8$. This implies that $\delta_{\rm D}$ used in Model A is close to the lower limit given by the transparency of $\gamma$ rays.

Motivated by the above results, two typical values of $\delta_{\rm D}=10$ and 12 are used, and they are referred as Model C and Model D, respectively. In Fig. 6, we compare the modeling results of Models A, C and D. The parameters are given in Table 1. We find that the location of the narrow feature moves to higher energy for higher value of $\delta_{\rm D}$, and a decrease of $B^{\prime}$ and an increasing of $T^{\prime}_{\rm f,2}$ is required to reproduce the TeV spectrum below 3 TeV. Meanwhile, a smaller value of $a=3.0$ is also needed for matching the observed X-ray spectrum and TeV spectrum below $\sim 3$ TeV. It is worth noting that the duration of injection $\Delta T^{\prime}_{\rm inj,2}=T^{\prime}_{\rm f,2}-T^{\prime}_{\rm s,2}$ can be predicted from the relation

where $\gamma^{\prime}_{\rm pk,2}\simeq\sqrt{\frac{\nu_{\rm pk,2}}{\nu_{0}B^{\prime}\delta_{\rm D}}}$ is the peak frequency of the synchrotron component below 1 keV. Using $\nu_{\rm pk,2}=5\times 10^{16}$ Hz, we obtain that $\Delta T^{\prime}_{\rm inj,2}/t^{\prime}_{\rm dyn}$ is 1.12, 2.14 and 4.89 for Model A, C and D, respectively. These values are in good agreement with that obtained from the SED modeling.

We find that the global properties of the jet including $B^{\prime}$, $\delta_{\rm D}$ and $R^{\prime}$, as well as $\gamma_{\rm eq}^{\prime}$ obtained in Model A are quite similar to that obtained by Acciari et al. (2020). Moreover, it can be seen from Table 1 that the value of the diffusion coefficient $D_{0}$ is $2.4\times 10^{-7}~{}\rm s^{-1}$. We find that similar result was obtained in other HBLs (e.g., Lewis, Becker & Finke, 2016) and Fermi bubbles (Sasaki, Asano & Terasawa, 2015).

In brief, we explain the broadband SED of Mrk 501 on 2014 July 19 by using a self-consistent one-zone leptonic jet model. In this model, the emitting EED is obtained by solving a Fokker-Planck equation including acceleration processes. Two injection episodes of electrons are needed to successfully explain the SED. The sharp and narrow spectral feature observed at $\sim$3 TeV is resulting from the pileup branch of electrons IC scattering synchrotron photons produced by the PL branch of electrons. The extremely high flux observed by Swift-BAT above 10 keV is interpreted as the synchrotron emission from the pileup bump of electrons. The pileup branch of electrons is the accelerated first-injection-electrons, and the PL branch of electrons is the accelerated second-injection-electrons. Moreover, we find that shock acceleration is required to dominate over the stochastic acceleration in modeling the narrow peak feature in TeV spectrum. This is different from the stochastic scenario proposed by Acciari et al. (2020).

The scenario of multiple injections of relativistic electrons has been proposed to provide a explanation for the extreme flux variabilities of blazars (Röken & Schlickeiser, 2009; Röken et al., 2018). Röken et al. (2018) argued that multiple-injection scenario is more realistic than single-injection scenario in blazar jets, as the blazar jets could extend over parsecs to tens of kiloparsecs scales, and thus likely pick up several electron populations from intergalactic media.

In principle, the narrow feature in the VHE spectrum of Mrk 501 is expected to be detected with high significance by CTA, which is designed to reach a sensitivity about an order of magnitude better than that of current imaging air Cherenkov telescopes and water Cherenkov telescopes (e.g., Actis et al., 2011; Acharya et al., 2013; Sol et al., 2013). The precise measurement of the TeV narrow feature by CTA could place strong constraints on the acceleration process working in the blazar jet.

We acknowledge the National Natural Science Foundation of China (NSFC-11803081, NSFC-12065011) and the joint foundation of Department of Science and Technology of Yunnan Province and Yunnan University [2018FY001(-003)]. The work of D. H. Yan is also supported by the CAS Youth Innovation Promotion Association and Basic research Program of Yunnan Province (202001AW070013).

No new data were generated or analysed in support of this research.