# Photohadronic interpretation of the different incarnation of 1ES 2344+514

Sarira Sahu<EMAIL_ADDRESS>Instituto de Ciencias Nucleares,
Universidad Nacional Autónoma de México,
Circuito Exterior, C.U., A. Postal 70-543, 04510 México DF, México  Isabel
Abigail Valadez Polanco<EMAIL_ADDRESS>Facultad de Ingeniería,
Universidad Autónoma de Yucatán,
Industrias No Contaminantes S/N, Sin Nombre de Col 27, Mérida, Yucatán, México
Subhash Rajpoot<EMAIL_ADDRESS>Department of Physics and
Astronomy, California State University,
1250 Bellflower Boulevard, Long Beach, CA 90840, USA

###### Abstract

Since its discovery in 1995, the high-energy peaked blazar 1ES 2344+514 has
undergone several episodes of GeV-TeV flaring and has been observed in the
multiwavelength by several other telescopes. The observed X-ray spectrum of
1996 and the flaring event of 2016 establish that 1ES 2344+514 has a temporary
EHBL-like behavior. Such behavior has also been observed in several nearby
high-energy peaked blazars. We use the photohadronic model to account for the
GeV-TeV flaring observed events of 1995 and 2007. Also, a recently proposed
two-zone photohadronic model, which is successful in explaining the multi-TeV
flaring events of many transient EHBL-like source, is employed to explain the
GeV-TeV flaring spectra of MJD 57611 and MJD 57612. We find that the zone-2
parameters of the two-zone photohadronic model play a central role in
explaining these spectra. Probably this is an indication of a new type of
transient EHBL-like source. We find that our fits to the observed spectra are
comparable or better than the other leptonic and hadronic models employed in
the literature to address the same issue.

High energy astrophysics (739), Blazars (164), Gamma-rays (637), Relativistic
jets (1390), BL Lacertae objects (158)

mainReferences appendixReferences (Appendix)

##  1 Introduction

Blazars are the primary sources of extragalactic $\gamma$-rays. Their spectral
energy distribution (SED) consists of two non-thermal peaks in the
$\nu$-$F_{\nu}$ plane (Dermer & Schlickeiser, 1993). The low energy peak in
the infrared to the X-ray frequency region is due to the synchrotron photons
produced from a population of relativistic electrons in the jet magnetic
field. It is believed that the high energy peak between the X-ray and the
$\gamma$-ray energies is produced either from the Synchrotron Self-Compton
(SSC) scattering of high-energy electrons and the self-produced synchrotron
photons in the jet (Maraschi et al., 1992; Gao et al., 2013) or due to the up-
scattering of low energy external seed photons and the ultra-relativistic
electrons in the jet. The sources of the seed photons could be either the
accretion disk or the broad line regions or the dusty torus (Błażejowski et
al., 2000; Sikora et al., 1994; Dermer & Schlickeiser, 1993). Recently,
different radiation mechanisms related to the leptonic and the hadronic
processes were reviewed in the context of AGN jets by Cerruti (2020). Blazars
are classified according to the position of the synchrotron peak frequency as,
low-energy peaked blazars (LBLs, $\nu_{peak}<10^{14}\,Hz$), intermediate-
energy peaked blazars (IBLs, $10^{14}\,Hz<\nu_{peak}<10^{15}\,Hz$), high
energy-peaked blazars (HBLs, $10^{15}\,Hz<\nu_{peak}\,<10^{17}\,Hz$) (Abdo et
al., 2010) and extreme energy peaked blazar (EHBL, $>10^{17}$ Hz) (Costamante
et al., 2001).

The EHBLs have unusually high synchrotron peak frequency and also the SSC peak
shifts towards the higher energy regime of $\gamma$-rays. Therefore, the shift
of the synchrotron peak and the SSC peak in the multiwavelength (MW) SED
generally translate to a particularly hard X-ray and a very high energy (VHE,
$>\,100$ GeV) $\gamma$-ray spectra. Low luminosity and limited variability are
the characteristic features of EHBLs which make them difficult to detect
(Costamante et al., 2018; Acciari et al., 2020). Thus, the population of the
observed EHBLs are much less than the observed HBLs (Aharonian, 2007;
Aharonian et al., 2007a, b). Different temporal bahavior has been observed
among the EHBLs. For example, 1ES 0229+200, 1ES 0347-232, RGB J0710+591 and
1ES 1101-232 (Aharonian, 2007; Costamante et al., 2018) are the archetypal
EHBLs which constantly exhibit extreme properties. However, some well known
and well studied nearby HBLs such as Markarian 421 (Mrk 421), Markarian 501
(Mrk 501) and 1ES 1959+650 belong to the EHBL family as temporary members
(Ahnen et al., 2018; Foffano et al., 2019). In 1996 and 2016 another
comparatively nearby HBL 1ES 2344+514 had also shown a temporary EHBL-like
behavior and a hard TeV spectra (MAGIC Collaboration et al., 2020a).
Particularly, during the flaring of 2016, both the synchrotron and the SSC
peaks were pushed towards higher energies (Ahnen et al., 2018; MAGIC
Collaboration et al., 2020b; Aleksić et al., 2015a). In March 2010, Mrk 421
was very active and various telescopes made observations in the MW range for
13 consecutive days (Aleksić et al., 2015a). The synchrotron peak position was
above $10^{17}$ Hz and the second peak was also at higher energies during this
observation period, thus exhibiting the signature for extreme HBL-like
behavior. Mrk 501 exhibited EHBL-like behavior during the VHE flaring of May -
July 2005 (Albert et al., 2007a) and June 2012 (Ahnen et al., 2018).
Similarly, during a MW campaign of 1ES 1959+650 between April 29 and November
21 of 2016, the MAGIC telescopes observed multi-TeV flaring events for the
nights of 13th, 14th of June and 1st of July 2016 when the position of the
synchrotron peak was in the EHBL range, and again, exhibiting extreme HBL-like
behavior (MAGIC Collaboration et al., 2020b). Several other authors have also
studied the MW SED of 1ES 1959+650 (Shah et al., 2021; Chandra et al., 2021).

The accumulating evidence from a large number of MW observations so far
suggests that the one-zone SSC model is too simple to account for the SEDs of
these objects. In the model, the emission from optical to VHE $\gamma$-rays
takes place in a single spherical blob. Moreover, in order to explain the
observed spectra, large values of the minimum electron Lorentz factor, the
bulk Lorentz factor (Ahnen et al., 2018; MAGIC Collaboration et al., 2020b),
and unrealistically low values of the magnetic field are used. In order to
avoid the above shortcomings, alternative models have been proposed, such as
the two-zone leptonic model, the inverse Compton scattering of relativistic
electrons with the cosmic microwave background, the spine-layer structured jet
model and different other hadronic models (MAGIC Collaboration et al., 2020a,
c).

Previously, the VHE spectra of many HBLs and EHBLs were successfully explained
using the photohadronic model (Sahu et al., 2019, 2020b). The photohadronic
model is based on the assumption that during the VHE flaring event, the Fermi
accelerated high energy protons in the inner jet region interact with the low
energy tail region of the SSC seed photons through the dominant
$\Delta$-resonance channel $p\gamma\rightarrow\Delta^{+}$ with the multi-pion
channels suppressed. Thus, due to this process, a change in the spectral
behavior of the SSC photons is expected in the VHE spectrum. The
$\Delta$-resonance subsequently decays to $\gamma$-rays via neutral pion
production in the jet environment. These $\gamma$-rays constitute the observed
VHE spectrum. However, given the kinematical conditons and the density of the
background seed photons in the internal jet, the photohadronic process works
well only for $E\gamma\gtrsim 100$ GeV. Below this energy the leptonic
processes have the most important contribution to the multiwavelength SED.
Here, we use the photohadronic process to interpret the VHE spectrum (Sahu,
2019). In several previous studies, we have shown that the flux in the low
energy tail region of the SSC spectrum can be expressed as a power-law which
is taken to be as $\Phi_{SSC}\propto E_{\gamma}^{-\beta}$ where $E_{\gamma}$
is the observed photon energy and $\beta$ is the spectral index. The value of
$\beta$ lies in the range $0<\beta\leq 1$ to explain the low emission state,
the high emission state and the very high emission state of the flaring HBLs
(Sahu et al., 2019). Recently, it was shown that the temporary EHBL-like
behavior of the flaring events of Mrk 421, Mrk 501 and 1ES 1959+650 are
difficult to explain by the photohadronic model (Sahu et al., 2020a, 2021a,
2021b) in its simplest form. This happens because $\Phi_{SSC}$ fails to behave
as a single power-law for the temporary EHBL. During this temporary EHBL-like
behavior, the synchrotron peak and the SSC peak are shifted towards higher
energies. Due to this shift, the standard jet parameters, such as the blob
size, the bulk Lorentz factor, and the strength of the magnetic field, may
differ from the ones used previously in HBL, but the mechanisms to accelerate
particles and the emission mechanisms are taken to be the same (Sahu et al.,
2021b). In this changed environment, we claim that the photohadronic process
should still be the dominant process. In order to explain the temporary EHBL-
like behavior of the VHE flaring events of the above objects, the
photohadronic model was extended to include two different zones for the SSC
flux: the first one with the conventional behavior with $0<\beta\leq 1$ and
the second one with $1<\beta\leq 1.5$ with no other modifications to the
proton spectral index in the jet. Previously, we have shown that in such
scenarios, the photohadronic process is favored over other mechanisms such as
the single zone SSC or the proton synchrotron that employ low magnetic field
and low cutoff energy for the Fermi accelerated protons.

The HBL 1ES 2344+514 has undergone several VHE flaring episodes since its
discovery in 1995. In 1996 the observation in X-rays showed that the
synchrotron peak was above $4\times 10^{18}$ Hz, a signature of EHBL-like
bahavior. Similarly, in August 2016 the ground-based $\gamma$-ray telescope
First G-APD Cherenkov Telescope (FACT) detected 1ES 2344+514 in a high
emission state which triggered multiwavelength observations. During this
observational period, it was observed that the VHE $\gamma$-ray spectrum and
the X-ray spectrum were hard and the frequency shift of the synchrotron peak
was towards $\gtrsim 10^{18}$ Hz, a renewal of extreme-EHBL like behavior of
1ES 2344+514. However, during 1995, 2007-2008 and 2014-2015, the VHE flaring
events were found to be of the standard HBL types which shows that the shift
towards EHBL-like behavior was temporary. Thus, in order to correctly account
for the standard HBL and EHBL we employ the one-zone photohadronic model for
the VHE flaring events for which the synchrotron peak is below $10^{17}$ Hz
and the two-zone photohadronic model to explain the VHE flaring events for
which the synchrotron peak is $\gtrsim\,10^{17}$ Hz. We apply the one-zone
photohadronic analysis to explain the VHE flaring events of December 1995 and
December 2007 of 1ES 2344+514 as the synchrotron peaks were in the HBL limit
($10^{15}-10^{16}$ Hz). We apply the two-zone photohadronic model to the
flaring events of a 11th and 12th of August 2016 as observed by the MAGIC
telescopes in which the source demonstrated EHBL-like behavior.

##  2 Flaring history of 1ES 2344+514

1ES 2344+514 was first detected in the Einstein Slew Survey (Elvis et al.,
1992) at X-ray energies and is classified in the HBL category (Perlman et al.,
1996). It was the third VHE extragalactic source detected after Mrk 421 and
Mrk 501 (Catanese et al., 1998) at a red shift of $z=0.044$ (Perlman et al.,
1996) and a luminosity distance of $d_{L}\sim 200$ Mpc. For the first time,
the VHE $\gamma$-rays from this object were observed during an intense flare
on December 20, 1995 by the Whipple 10 m telescope and the observed flux above
350 GeV was almost 60% of the Crab Nebula flux (Catanese et al., 1998;
Schroedter et al., 2005). In 1996 a flare in X-ray energy was observed with a
large flux variability on timescales of a few hours when the source was at its
brightest (Giommi et al., 2000). During this period the lower limit of the
synchrotron peak frequency was shifted above $3\times 10^{18}$ Hz, a clear
signature of EHBL-like behavior of the source seen for the first time. This
extreme behavior of the source motivated several multiwavelength campaigns to
be undertaken on a regular basis as part of the blazar science programme to
observe and to model the 1ES 2344+514 broad-band SED using the simultaneous
and the quasi-simultaneous data (Albert et al., 2007b; Godambe et al., 2007;
Acciari et al., 2011; Aleksić et al., 2013). Subsequently, it was found that
the source was mostly in a low emission state during the following years with
no signature of any extreme behavior either in the X-ray band or in the
$\gamma$-ray (Catanese et al., 1998; Giommi et al., 2000; Allen et al., 2017).
The MW observations were performed from October 2007 to January 2008 and a
strong VHE flare was observed on December 7, 2007 by VERITAS (Acciari et al.,
2011). The measured flux was about 48% of the Crab Nebula flux. On 10th of
August 2016 the $\gamma$-ray telescope FACT detected a high emission state
which triggered the MW observations from radio to VHE $\gamma$-rays (MAGIC
Collaboration et al., 2020a).

Using the simultaneous observations of the flaring source 1ES 2344+514 from
the MW, the broad-band SED was constructed and studied using both the leptonic
models and the hadronic models. Recently, the VHE spectrum of MJD 57612 (12th
August 2016) is also fitted using the one-zone SSC model, the two component
SSC model and the proton-synchrotron model (MAGIC Collaboration et al., 2020a,
c). During this observation period the VHE $\gamma$-ray flux level above 300
GeV was about 55% of the Crab Nebula flux and the synchrotron peak frequency
was shifted to $\gtrsim 10^{18}$ Hz, signalling an extreme-HBL behavior.
Although these models explain the SED well, the main difference in them lies
in the strength of the magnetic field used, the energy in different components
and the total number of free parameters. The one-zone SSC model employs a
smaller magnetic field than the two-component model. Also, the bulk Lorentz
factor in the one-zone SSC model is much higher than the two-component model.
Low magnetic field corresponds to low magnetic field energy density leading to
the electron energy density dominating over the magnetic one. The two-
component model is based on the co-acceleration of protons and electrons in
the jet and energy is transferred from the protons to the electrons in the
shock transition layer (Zech & Lemoine, 2021). The core component of the model
provides seed photons for the Compton scattering, while X-rays to VHE gamma-
rays are mainly produced in the blob region (MAGIC Collaboration et al.,
2020c). Although the two-component models fit better to the observed data,
they are difficult to constrain due to the large number of free parameters
used. On the other hand, the proton synchrotron model requires ultra-high
energy protons and a very high magnetic field to produce VHE gamma rays. A
major weakness of the proton synchrotron model is that an unrealistic super-
Eddington luminosity in protons is needed to explain the observed VHE
spectrum. However, in the photohadronic scenario this problem is evaded due to
the higher photon density in the inner jet region (Sahu, 2019).

##  3 photohadronic model

We recall the salient features of the photohadronic model (Sahu et al., 2019)
used to explain the VHE flaring events of HBLs. As discussed earlier, when
there is a temporary transition from HBL to EHBL, as in the case of Mrk 421,
Mrk 501 and 1ES 1959+650, the one-zone photohadronic model is inadequate to
explain the spectral behavior. In order to remedy this shortcoming the two-
zone photohadronic model is employed (Sahu et al., 2021a, 2020a, b) which is
an extension of the one-zone photohadronic model. Depending on the nature of
the VHE flaring event of the blazar, the low energy tail region of the SSC
flux can either be expressed as a single power-law or a two component power-
law. We refer to the former as the one-zone photohadronic model and the latter
as the two-zone photohadronic model. These two models are briefly discussed
below.

###  3.1 One-zone photohadronic model

The VHE flaring from a number of HBLs are successfully explained in the
context of the photohadronic model. This model has also successfully explained
the VHE flaring from atypical EHBLs, such as 1ES 0229+200 and 1ES 0347-232
(Sahu et al., 2019). In order to explain the VHE flaring events from HBLs, the
photohadronic model assumes a double jet structure formed along the same axis.
In this scenario a small and compact jet of size $R^{\prime}_{f}$ is enclosed
by a bigger outer jet of size $R^{\prime}_{b}$ with
$R^{\prime}_{f}<R^{\prime}_{b}$. The photon densities in the inner and the out
jet are taken to be as $n^{\prime}_{\gamma,f}$ and $n^{\prime}_{\gamma}$
respectively where $n^{\prime}_{\gamma,f}>n^{\prime}_{\gamma}$ (′ implies in
the comoving frame). The inner jet and the outer jet have bulk Lorentz factors
$\Gamma_{in}$ and $\Gamma_{ext}$ respectively and for simplicity we assume
$\Gamma_{in}\simeq\Gamma_{ext}\simeq\Gamma$ with a common Doppler factor
$\mathcal{D}$ (Sahu, 2019).

During the multi-TeV flaring, the Fermi accelerated high energy protons are
injected into the inner jet region with a power-law spectrum $dN/dE_{p}\propto
E^{-\alpha}_{p}$, where $E_{p}$ is the proton energy and the spectral index
$\alpha\geq 2$. These protons interact with the low energy tail region of the
SSC seed photons to produce $\Delta$-resonances. The $\Delta$-resonance decays
to $\gamma$-rays and neutrinos via the decay of neutral and charged pions
respectively. In this decay process, each neutrino carries about half the
photon energy. From the IceCube point of view, the neutrino energy is very low
for the neutrinos to be detected.

The $\Delta$-resonance is produced from the $p\gamma$ interaction when the
following condition is satisfied

$E_{\gamma}\epsilon_{\gamma}=\frac{0.032\ \Gamma{\mathcal{D}}}{(1+z)^{2}}\
\mathrm{GeV^{2}}.$ (1)

Here $\epsilon_{\gamma}$ and $E_{\gamma}$ are the seed photon energy and the
VHE $\gamma$-ray photon energy in the observer’s frame. The quantities
$\Gamma$, ${\cal D}$ and $z$ are respectively the bulk Lorentz factor, the
Doppler factor and the redshift of the object. In this process about 10% of
the proton energy $E_{p}$ is taken by the VHE photon ($E_{\gamma}\simeq
0.1\,E_{p}$). The observed VHE $\gamma$-ray flux $F_{\gamma}$ is proportional
to the seed photon number density $n^{\prime}_{\gamma,f}$ in the inner jet
region and the incident high energy proton flux $F_{p}\equiv
E^{2}_{p}\,dN/dE_{p}$. This gives $F_{\gamma}\propto
n^{\prime}_{\gamma,f}\,F_{p}$. For all the VHE flaring events observed so far
from the HBLs, the range of $E_{\gamma}$ corresponds to a range of
$\epsilon_{\gamma}$ that can be gotten from Eq. (1), and the latter always
lies in the low energy tail region of the SSC band. The photon density in the
outer jet is given by

$n^{\prime}_{\gamma}(\epsilon_{\gamma})=\eta\left(\frac{d_{L}}{R^{\prime}_{b}}\right)^{2}\frac{1}{(1+z)}\frac{\Phi_{SSC}(\epsilon_{\gamma})}{{\cal
D}^{2+\kappa}\,\epsilon_{\gamma}},$ (2)

where $\eta$ corresponds to the efficiency of the SSC process. We take
$\eta=1$ for 100% efficiency. $d_{L}$ is the luminosity distance to the source
and $\kappa=0(1)$ corresponds to a continuous (discrete) blazar jet. As the
inner jet expands adiabatically into the outer jet, the photon density in the
former decreases. Due to uncertainties in determining the photon density in
the inner region, a simple scaling law is assumed for the photon densities in
the inner jet and the outer jet as follows,

$\frac{n^{\prime}_{\gamma,f}(\epsilon_{\gamma,1})}{n^{\prime}_{\gamma,f}(\epsilon_{\gamma,2})}\simeq\frac{n^{\prime}_{\gamma}(\epsilon_{\gamma,1})}{n^{\prime}_{\gamma}(\epsilon_{\gamma,2})}.$
(3)

This equation implies that the ratio of the photon densities for two different
background energies $\epsilon_{\gamma,1}$ and $\epsilon_{\gamma,2}$ in the
inner and the out jet regions are almost the same. From Eq. (2) and Eq. (3),
$n^{\prime}_{\gamma,f}$ can be expressed in terms of $\Phi_{SSC}$. This gives
$F_{\gamma}\propto n^{\prime}_{\gamma}$ implying that
$F_{\gamma}\propto\Phi_{SSC}\epsilon_{\gamma}^{-1}$. Also, $F_{\gamma}\propto
F_{p}$ implies $F_{\gamma}\propto E_{\gamma}^{-\alpha+2}$. Thus in the low
energy tail region of the SSC band the flux follows a perfect power-law
behavior determined as $\Phi_{SSC}\propto\epsilon^{\beta}_{\gamma}\propto
E_{\gamma}^{-\beta}$. The model based on it is referred to as the one-zone
photohadronic model.

The role of extragalactic background light (EBL) in the propagation of VHE
$\gamma$-rays is crucial as it attenuates the propagating VHE photons by
producing $e^{+}e^{-}$. Thus, in the observation of VHE gamma-rays, EBL
effects have to be properly accounted for. Using different methods several EBL
models have been developed (Salamon & Stecker, 1998; Franceschini et al.,
2008; Domínguez et al., 2011; Gilmore et al., 2012; Domínguez et al., 2013).
With the EBL correction, the observed VHE spectrum in the one-zone
photohadronic model is given by

$F_{\gamma}(E_{\gamma})=F_{0}\left(\frac{E_{\gamma}}{TeV}\right)^{-\delta+3}\,e^{-\tau_{\gamma\gamma}(E_{\gamma},z)}=F_{\gamma,in}(E_{\gamma})\,e^{-\tau_{\gamma\gamma}(E_{\gamma},z)}.$
(4)

Here $F_{0}$ is the normalization constant which can be adjusted from the
observed VHE spectrum. In this model the only free parameter is the spectral
index $\delta=\alpha+\beta$. $F_{\gamma,in}$ and $\tau_{\gamma\gamma}$ are
respectively the intrinsic VHE flux and the optical depth for the lepton pair
production process $\gamma\gamma\rightarrow e^{+}e^{-}$. For the multi-TeV
flaring epochs of the HBLs, the value of $\delta$ lies in the range
$2.5\leq\delta\leq 3.0$ (Sahu et al., 2019).

###  3.2 Two-zone photohadronic model

From recent studies on the EHBL-like behavior of the nearest blazars Mrk 421,
Mrk 501 and 1ES 1959+650, it is shown that the canonical photohadronic model
is not adequate to explain the VHE flaring epochs of these sources. To
overcome this difficulty, it is assumed that the tail region of the SSC band
which is responsible for the production of $\Delta$-resonance has two
different power-laws:

$\Phi_{SSC}\propto\left\\{\begin{array}[]{cr}E^{-\beta_{1}}_{\gamma},&\quad\mathrm{100\,GeV\,\lesssim
E_{\gamma}\lesssim E^{intd}_{\gamma}}\\\
E^{-\beta_{2}}_{\gamma},&\quad\mathrm{E_{\gamma}\gtrsim E^{intd}_{\gamma}}\\\
\end{array}\right..$ (5)

The spectral indices satisfy $\beta_{1}\neq\beta_{2}$ and $E^{intd}_{\gamma}$
is the transition energy between zone-1 and zone-2. The value of
$E^{intd}_{\gamma}$ is different for different spectra. Also, the value of
$E^{intd}_{\gamma}$ provides information on the nature of the flaring event.
Using Eq. (5), we can express the observed VHE spectrum as,

$F_{\gamma}=e^{-\tau_{\gamma\gamma}}\times\begin{cases}F_{1}\,\left(\frac{E_{\gamma}}{TeV}\right)^{-\delta_{1}+3},&\quad\mathrm{100\,GeV\,\lesssim
E_{\gamma}\lesssim E^{intd}_{\gamma}}\,\,(\text{zone-1})\\\
F_{2}\,\left(\frac{E_{\gamma}}{TeV}\right)^{-\delta_{2}+3},&\quad\,\,\,\,\,\,\,\,\,\,\,\,\mathrm{E_{\gamma}\gtrsim
E^{intd}_{\gamma}}\,\,\,(\text{zone-2})\end{cases}.$ (6)

Thus, we now have two different zones and the normalization constants for
zone-1 and zone-2 are defined as $F_{1}$ and $F_{2}$ respectively. As
explained earlier, the photohadronic model is well suited for
$E_{\gamma}\gtrsim 100\,GeV$. Below this energy it is the leptonic model which
contributes to the SED. In the context of the photohadronic model, we only
concentrate on the VHE spectrum. The observed VHE spectrum of the EHBL can be
fitted in zone-1 and zone-2 separately by adjusting the spectral indices
$\delta_{i}=\alpha+\beta_{i}$ ($i=1,2$) respectively. As shown, the VHE events
from HBLs can be explained with $2.5\leq\delta\leq 3.0$ and the temporary
EHBLs with the two-zone model. However, there can be another type of VHE
flaring from blazars. In this situation, the spectrum can be explained by
zone-2 only, i.e.

$F_{\gamma}=e^{-\tau_{\gamma\gamma}}\times
F_{2}\,\left(\frac{E_{\gamma}}{TeV}\right)^{-\delta_{2}+3}.$ (7)

In this case the value of the spectral index will be $3.1\leq\delta_{2}\leq
3.5$.

We take the proton spectral index $\alpha=2$ (Dermer & Schlickeiser, 1993). By
fixing the value of $\alpha$ we automatically fix the value of $\beta_{i}$ for
a given $\delta_{i}$. For the EBL correction to the VHE spectrum the EBL model
of Domínguez et al. (2011) and Franceschini et al. (2008) are widely used. For
the analysis of our results, we use the EBL model of Franceschini et al.
(2008). This model is consistent with a EBL peak flux of $\sim 15\
\mathrm{nW}\mathrm{n^{-2}}\ \mathrm{sr^{-1}}$ at $1.4\,\mathrm{\mu m}$
(Aleksić et al., 2015b; Ahnen et al., 2017) and the results are similar to the
EBL model of Domínguez et al. (2011). In earlier work (Sahu et al., 2019), we
have shown that the value of $\delta_{1}$ (zone-1) allows one to characterize
the flaring state for a given observation and is always in the range
$2.5\leq\delta_{1}\leq 3.0$. The standard VHE flaring events of all the HBLs
lie in this range of $\delta_{1}$. To fit the data in zone-1 (low energy
regime) we take this range of $\delta_{1}$ and fit it with different values of
$F_{1}$. Zone-2 (high energy regime) is similarly fitted for
$3.1\leq\delta_{2}\leq 3.5$ with different values of $F_{2}$. We select the
best fit values of both the zones by calculating their statistical
significance, the point of intersection of both the curves gives the
transition energy $E^{intd}_{\gamma}$.

Figure 1: The VHE flaring event of December 20, 1995 (MJD 50071) observed by
the Whipple telescope (Wh95) (Schroedter et al., 2005) is fitted with the
photohadronic model. Here and in all the subsequent figures the normalization
constants $F_{i}$ ($i=0,2$) are expressed in units of
$\mathrm{erg\,cm^{-2}\,s^{-1}}$. The best fit to the VHE spectrum is shown by
the blue curve. The blue shaded region corresponds to
$F_{0}=8.13^{+6.83}_{-2.97}\times 10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$ and
$\delta=2.9^{+0.67}_{-0.31}$. The dashed curve corresponds to intrinsic flux.
For reference, we have also shown the upper limit of the VHE flux for the
night of 5 December 1996 observed by Whipple telescope.

##  4 results

The first VHE $\gamma$-ray flaring event from the HBL 1ES 2344+514 was
observed on 1995 December 20 by the Whipple 10 m telescope (Wh95) and shown
that it can be fitted with a power-law (Schroedter et al., 2005). During this
period the observed $E_{\gamma}$ was in the energy range $1.15\,TeV\lesssim
E_{\gamma}\lesssim 9.1\,TeV$ and this VHE emission was about 50 times brighter
than the emission during the quiescent phase of this object measured in later
years. This flux level is the historical maximum observed so far from 1ES
2344+514 and the spectrum is fitted well with $\delta=2.9$ and
$F_{0}=8.13\times 10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$ and this
corresponds to a high emission state. $F_{0}$ and $\delta$ are two parameters
in the model and the spectral index is constrained to be in the range
$2.5\leq\delta\leq 3.0$. By using the minimum $\chi^{2}$ fit and varying
simultaneously both $F_{0}$ and $\delta$, the $1\sigma$ confidence interval
defined by $\chi^{2}_{min}+2.3$ is calculated and the errors in these
parameters are estimated for all the flaring epochs. For the above spectrum,
the errors in the parameters are given by $F_{0}=8.13^{+6.83}_{-2.97}\times
10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$ and $\delta=2.9^{+0.67}_{-0.31}$
respectively. In Fig. 1, the blue shaded region corresponds to the errors in
the parameters $F_{0}$ and $\delta$. These are shown in the same figure along
with the best fit values. In this figure we have also shown the intrinsic
spectrum $F_{\gamma,in}$ (dashed curve) which is proportional to
$E^{0.1}_{\gamma}$. The range $1.15\,TeV\lesssim E_{\gamma}\lesssim 9.1\,TeV$
corresponds to the Fermi accelerated proton energy in the range
$11.5\,TeV\lesssim E_{p}\lesssim 91\,TeV$. $\delta=2.9$ corresponds to the
flare in a high emission state and the spectral index of the background SSC
photon is $\beta=0.9$. A very small magnetic field $(B\sim 10^{-4}$ G) is
needed to accelerate the protons up to $\sim 100$ TeV energy in the inner jet
region.

The central black hole (BH) mass of 1ES 2344+514 derived from the stellar
velocity dispersion measurements is $M_{BH}\sim 10^{8.80\pm 0.16}\,M_{\odot}$
(Barth et al., 2003) which corresponds to an Eddington luminosity
$L_{Edd}\sim(0.55-1.15)\times 10^{47}\,\mathrm{erg\,s^{-1}}$. The integrated
flux in the energy range $1.15\,TeV\lesssim E_{\gamma}\lesssim 9.1\,TeV$ is
$F_{\gamma}\sim 8.53\times 10^{-11}\mathrm{erg\,cm^{-2}\,s^{-1}}$,
corresponding to a luminosity of $L_{\gamma}\sim 4.1\times
10^{44}\,\mathrm{erg\,s^{-1}}$.

We can estimate the bulk Lorentz factor $\Gamma$ for this flaring event. As
shown in Fig. 7 of Acciari et al. (2011), the tail region of the SSC starts
around $10^{20}$ Hz. By assuming that the 9.1 TeV observed $\gamma$-ray is
produced from the interaction of 91 TeV proton with the SSC seed photons of
energy $\epsilon_{\gamma}\gtrsim 2\times 10^{20}$ Hz ($\gtrsim 0.83$ MeV), we
obtain the bulk Lorentz factor $\Gamma\gtrsim 16$ which is consistent with the
X-ray variability time scale of about an hour as seen from 1ES 2344+514
(Giommi et al., 2000).

To produce the $\Delta$-resonance within the jet, we take $R^{\prime}_{f}\sim
5\times 10^{15}\,\mathrm{cm}$ and $R^{\prime}_{b}\sim 10^{16}\,\mathrm{cm}$
and the optical depth $\tau_{p\gamma}$ for the $p\gamma\rightarrow\Delta^{+}$
process should satisfy $\tau_{p\gamma}<1$. Also the proton luminosity
$L_{p}=7.5\tau_{p\gamma}^{-1}L_{\gamma}$ should satisfy $L_{p}<L_{Edd}/2$. By
taking $L_{Edd}\sim 10^{47}\,\mathrm{erg\,s^{-1}}$ (Barth et al., 2003), we
obtain $\tau_{p\gamma}>0.06$. By taking $\tau_{p\gamma}\sim 0.2$, the SSC
photon density in the inner jet region is $n^{{}^{\prime}}_{\gamma,f}\sim
8.0\times 10^{10}\,\mathrm{cm^{-3}}$ and $L_{p}\sim 1.5\times
10^{46}\,\mathrm{erg\,s^{-1}}$. Similarly, we can estimate the photon density
and the $L_{p}$ for other flaring epochs also.

In the year 1996 the X-ray spectrum observed was in the EHBL range (especially
the synchrotron peak) and also during this period the Whipple telescope
provided the upper limit to VHE spectrum which we show in Fig. 1 for
comparison. In Fig. 7 of Albert et al. (2007b), the Wh95 spectrum is fitted
with the one-zone SSC model. However, we observe that this spectrum is not
consistent with the original spectrum (Schroedter et al., 2005). Therefore, we
do not compare our results with the one-zone leptonic model fit of Albert et
al. (2007b).

Figure 2: The VHE flaring event of 7th December 2007 (MJD 54441) observed by
VERITAS telescopes is fitted with the photohadronic model. The blue shaded
region corresponding to the errors in the parameters $F_{0}$ and $\delta$ are
shown at $1\sigma$ confidence interval. The dashed curve corresponds to
intrinsic flux. For comparison, we have also shown the one-zone SSC fit
(Acciari et al., 2011)

The VERITAS telescopes observed 1ES 2344+514 for 37 nights during a
multiwavelength observation period between 4th October 2007 and 11th January
2008. On 7th of December 2007, VERITAS observed a strong VHE flare and the
flux above 300 GeV was about 48% that of the Crab Nebula flux. This VHE flare
was observed in the energy range $0.48\,TeV\lesssim E_{\gamma}\lesssim
2.3\,TeV$. Also, during this period, the synchrotron peak was in the HBL
range. In the context of the photohadronic model, the VHE spectrum is fitted
very well with the spectral index $\delta=3.0$ and the normalization factor
$F_{0}=4.17\times 10^{-11}\mathrm{erg\,cm^{-2}\,s^{-1}}$. The errors in the
parameters are given by $F_{0}=4.17^{+0.70}_{-0.79}\times
10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$ and $\delta=3.0^{+0.50}_{-0.24}$
respectively. In Fig. 2, we present the fit to the spectrum, the intrinsic
spectrum and the blue shaded region in the figure corresponds to the $1\sigma$
errors in $F_{0}$ and $\delta$. The spectrum is also fitted with the one-zone
SSC model Acciari et al. (2011). It is observed that our fit is marginally
better than the one-zone SSC fit as shown in Fig. 2. For $\Gamma={\cal D}=16$,
we obtain $3.3\,MeV\lesssim\epsilon_{\gamma}\lesssim 15.7\,MeV$.

Also, the fitted spectrum in the one-zone SSC model falls faster than the
photohadronic fit for $E_{\gamma}\gtrsim 2.3$ TeV. It was observed that by
excluding the flux of the flaring event, the average flux during this
observation period was found to be low. In the photohadronic model, the
integrated flux and the corresponding luminosity in the energy range
$0.48\,TeV\lesssim E_{\gamma}\lesssim 2.3\,TeV$ are respectively
$F_{\gamma}\sim 4.28\times 10^{-11}\mathrm{erg\,cm^{-2}\,s^{-1}}$ and
$L_{\gamma}\sim 2.05\times 10^{44}\,\mathrm{erg\,s^{-1}}$.

Figure 3: The VHE spectrum of MJD 57611 (11th August 2016) observed by MAGIC
is fitted with the photohadronic model ($\delta=3.0$ and $F_{0}=4.76\times
10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$) and with the two-zone photohadronic
model ($\delta_{2}=3.1$ and $F_{2}=4.31\times
10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$ corresponds to zone-2). Here zone-1
is absent. The dashed curve is the intrinsic spectrum corresponding to two-
zone model. We have also shown in the blue shaded region the errors in $F_{2}$
and $\delta_{2}$ at the $1\sigma$ confidence interval.

The EHBL-like flaring event of 10th August 2016 was first observed in VHE by
FACT and follow-up observations were undertaken in MW by various telescopes.
Particularly in VHE it was observed by the MAGIC telescopes on MJD 57611 and
MJD 57612 (11th and 12th of August 2016). For MJD 57611 the source was
observed for 0.39 hours in the energy range $0.09\,TeV\lesssim
E_{\gamma}\lesssim 2.5\,TeV$ and high flux was recorded. On MJD 57612, the
source was observed for 0.48 hours and the observed integral VHE flux above
300 GeV was $(2.1\pm 0.4)\times 10^{-11}\mathrm{cm^{-2}\,s^{-1}}$ which is
about three times lower than the flux observed during the previous night. The
synchrotron peak frequency observed during this period was above $10^{18}$ Hz
and according to the classification scheme this belongs to the EHBL category.
The flaring events of MJD 57611 and MJD 57612 have similar physical conditions
to that of EHBL-like VHE flaring events of Mrk 421, Mrk 501 and 1ES 1959+650
and the latter events are explained with the two-zone photohadronic scenario
(Sahu et al., 2020a, 2021a, 2021b). Thus it is imperative to use again the
two-zone photohadronic model to interpret the VHE spectra of MJD 57611 and MJD
57612 where the zone-1 has spectral index $\delta_{1}$ and the zone-2 has
spectral index $\delta_{2}$. Both these zones have the common transition
energy $E^{intd}_{\gamma}$. We first tried to fit the spectrum of MJD 57611
with the two-zone photohadronic model as discussed above by using Eq. (6).
However, it is observed that the whole spectrum can be fitted perfectly well
only with the parameters corresponding to zone-2, i.e., $\delta_{2}=3.1$ and
its corresponding normalization constant $F_{2}=4.31\times
10^{-11}\mathrm{erg\,cm^{-2}\,s^{-1}}$ which is shown in Fig. 3. Thus, in this
case zone-1 appears to be depleted and therefore there is no need for
$E^{intd}_{\gamma}$.

Recently, from the EHBL-like VHE flaring of Mrk 421 (Sahu et al., 2021b) which
was observed for 13 consecutive days, we have shown that the VHE flaring event
decreases from a high emission state to a low emission state, and
correspondingly, the transition energy $E^{intd}_{\gamma}$ starts with $\sim
1$ TeV and slowly decreases to $\sim 0.25$ TeV. Consequently, during this
transition process, zone-2 became wider and spread into zone-1, making the
latter region narrower at the end of the flaring process.

It was expected that after the first observation of the EHBL-like behavior on
10th of August 2016 from 1ES 2344+514, the VHE spectrum observed in the
subsequent two days (11th and 12th) should be fitted with two different zones
with a transition energy $E^{intd}_{\gamma}$ and in the subsequent days
$E^{intd}_{\gamma}$ should have moved towards lower energy. However,
unexpectedly in the case of MJD 57611, $E^{intd}_{\gamma}$ does not exist and
zone-1 is absent in the spectrum. Thus, for the spectrum of MJD 57611, it is
possible that the zone-1 became too narrow or was completely taken over by
zone-2 in one day, making the spectrum fit well with the spectral index of
zone-2 only.

We have also fitted the spectrum using the one-zone photohadronic model and
good fit is obtained for $\delta=3.0$ and $F_{0}=4.76\times
10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$ which is consistent with a low
emission state. Even though the spectrum seems to be compatible with the low
emission state interpretation (with $\delta=3.0$) of a HBL, the physical
conditions of the flaring event are those of the EHBL type. Thus, the standard
HBL interpretation is discarded. For comparison, we have also shown this fit
in Fig. 3. The blue shaded region in this figure corresponds to the errors in
parameters $\delta_{2}=3.1^{+0.19}_{-0.13}$ and
$F_{2}=4.31^{+0.82}_{-1.02}\times 10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$ at
the $1\sigma$ confidence level. The integrated flux on MJD57611 and its
corresponding luminosity are respectively, $F_{\gamma}=1.25\times
10^{-10}\mathrm{erg\,cm^{-2}\,s^{-1}}$ and $L_{\gamma}\sim 6.0\times
10^{44}\,\mathrm{erg\,s^{-1}}$.

Figure 4: The VHE spectrum of MJD 57612 (12th August 2016) is fitted with the one-zone ($\delta=3.0$ and $F_{0}=1.93\times 10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$) and the two-zone ($\delta_{2}=3.1$ and $F_{2}=1.82\times 10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$) photohadronic model. The dashed curve is the intrinsic spectrum corresponding to the two-zone photohadronic model. The blue shaded region in this figure corresponds to the errors in parameters $\delta_{2}=3.1^{+0.42}_{-0.32}$ and $F_{2}=1.82^{+0.56}_{-0.65}\times 10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$ at the $1\sigma$ confidence level. We compare the photohadronic fit with the one-zone SSC model (SSC1a and SSC1b), the two component SSC model (SSC2) and the proton-synchrotron model (PS) of MAGIC Collaboration et al. (2020a) and MAGIC Collaboration et al. (2020c). The proton synchrotron model obtained by varying the maximum proton-synchrotron peak frequency, radius of the emitting region and the proton normalization is represented by the gray shaded region.  Table 1: In Fig. 4 the VHE spectrum of MJD 57612 is fitted with various models: the two different one-zone SSC model (SSC1a and SSC1b), the two-component SSC model (SSC2) (blob+core), the proton synchrotron model (PS) (MAGIC Collaboration et al., 2020a, c) and the photohadronic model (PH). The various parameters, like the bulk Lorentz factor ($\Gamma$), the spectral index ($\delta$), the blob Radius ($R^{\prime}_{b}$ in units of $10^{16}$ cm ) and the magnetic field ($B^{\prime}$ in G) used in these models are summarized below. The value of the spectral index $\delta$ given here is for the photohadronic model for MJD 57611 and MJD 57612 when the object was behaving as EHBL.  Parameter | SSC1a | SSC1b | SSC2 | PS | PH
---|---|---|---|---|---
|  |  | blob, core |  |
$\Gamma$ | 30 | 20 | 6, 4 | 30 | $\gtrsim 16$
$\delta$ | 2 | 2.07 | –,– | 1.5-2.5 | 3.1
$R^{\prime}_{b}$ | 1 | 1.22 | 1.07, 16.0 | 0.09-0.19 | 1
$B^{\prime}$ | 0.02 | 0.02 | 0.1, 0.1 | 48-63 | $\sim 10^{-4}$

The observed VHE spectrum of MJD 57612 was in the energy range
$0.14\,TeV\lesssim E_{\gamma}\lesssim 2.1\,TeV$ and it has exactly the same
behavior as MJD 57611. The observed spectrum on this night is fitted very well
by zone-2 parameters $\delta_{2}=3.1$ and $F_{2}=1.82\times
10^{-11}\mathrm{erg\,cm^{-2}\,s^{-1}}$. The errors in these parameters at
$1\sigma$ confidence level are given by $\delta_{2}=3.1^{+0.42}_{-0.32}$ and
$F_{2}=4.31^{+0.56}_{-0.65}\times 10^{-11}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$.
The best fit to the spectrum and the blue shaded region corresponding to the
errors in $F_{2}$ and $\delta_{2}$ are shown in Fig. 4 along with the
intrinsic spectrum. For comparison, we have shown the one-zone photohadronic
fit with $\delta=3.0$ and $F_{0}=1.93\times
10^{-11}\mathrm{erg\,cm^{-2}\,s^{-1}}$ which fits well to the observed
spectrum as shown in Fig. 4. However, in this case also, we have discarded
this fit for the same reason given for the spectrum of previous night (MJD
57611).

In previous studies, the VHE spectrum of MJD 57612 was fitted using the one-
zone SSC model, the two-component SSC model (SSC2) and the proton-synchrotron
model (PS) (MAGIC Collaboration et al., 2020a, c). In the one-zone SSC model,
two different sets of parameters were used to fit the spectrum. For
convenience, we make use of the following labels, SSC1a (MAGIC Collaboration
et al., 2020a) and SSC1b (MAGIC Collaboration et al., 2020c) as reference
labels. Although, the overall fit to the spectrum is good in all these models
as shown Fig. 4 , they behave differently in the high energy regime and also
the parameters used are different. In Table 1, we have summarized the
parameters used in these models. For $E_{\gamma}\gtrsim 1$ TeV all models
start behaving differently as shown in Fig. 4. The SSC2 fit falls faster than
the rest of the fits (MAGIC Collaboration et al., 2020c). The proton-
synchrotron model of Cerruti et al. (2015) is employed to fit the VHE
spectrum. By varying the maximum proton-synchrotron peak frequency, radius of
the emitting region and the proton normalization, the shaded region is
constructed in Fig. 4. Although the leptonic and the proton-synchrotron models
fit the data well, the values of the magnetic field used are different and the
difference commences for $E_{\gamma}\gtrsim 2$ TeV. The photohadronic fit lies
between the SSC2 fit and the SSC1a fit (MAGIC Collaboration et al., 2020a).
Also, for this night the integrated flux and the luminosity are respectively,
$F_{\gamma}=4.21\times 10^{-11}\mathrm{erg\,cm^{-2}\,s^{-1}}$ and
$L_{\gamma}\sim 2.02\times 10^{44}\,\mathrm{erg\,s^{-1}}$.

In the context of the photohadronic scenario, the characteristic feature of
the two flaring events (MJD 57611 and MJD 57612) is that, during the temporary
EHBL-like behavior of the source, their intrinsic spectra are proportional to
$E^{-0.1}_{\gamma}$ and the seed photon flux in the SSC region follows
$\Phi_{SSC}\propto\epsilon^{1.1}_{\gamma}$. The spectra of both the nights are
fitted very well with the parameters from zone-2 ($\delta_{2}$ and $F_{2}$)
only. It is important to mention that, had the source not been observed in the
X-ray region to see the synchrotron peak above $10^{18}$ Hz, it would have
easily mislead one to interpret the VHE flaring events as standard HBL flaring
events as both these spectra can be explained well with the photohadronic
model for $\delta=3.0$, corresponding to low emission state (as shown in Figs.
3 and 4). Previously, in the context of Mrk 421, Mrk 501 and 1ES 1959+650, we
have shown that the VHE spectrum from temporary EHBL-like behavior can be
explained by the two-zone photohadronic model where the spectral indices
$\delta_{1}$ and $\delta_{2}$ contribute to the fit of the spectrum (Sahu et
al., 2020a, 2021a, 2021b). However, to our knowledge, fitting the VHE spectrum
of a temporary EHBL with the spectral index $\delta_{2}$ only, so far has not
been observed in any of the EHBL-like events identified in other nearby
objects. Probably this is an indication of a new type of temporary EHBL whose
VHE spectrum is explained by a single zone-2 spectral index with
$\delta\gtrsim 3.1$.

In the photohadronic model, neutrino energy expectations are
$E_{\nu}=0.5\,E_{\gamma}$. On December 20, 1995, the HBL 1ES 2344+514 attained
the highest flux level and during this flaring period, the maximum
$\gamma$-ray energy observed was $E_{\gamma}\simeq 9.1\,TeV$. This corresponds
to $E_{\nu}\simeq 4.5\,TeV$ which is very low for neutrinos to be observed in
a neutrino detector like IceCube. Furthermore, the luminosity distance of the
HBL 1ES 2344+514 is $d_{L}\simeq 200\,Mpc$ which will further deplete the
neutrino flux. During the flaring events of December 2007 and August 2016, the
maximum observed $E_{\gamma}$ was much lower than $9.1\,TeV$ which correspond
to $E_{\nu}\ll 4.5\,TeV$. This implies that it is difficult to observe high
energy neutrinos from the above VHE flaring events.

##  5 Conclusion

In conclusion, the HBL 1ES 2344+514 was first detected in 1995 and since then
several episodes of VHE flaring have been reported. In 1996 a flare in the
X-ray was observed and the synchrotron peak of the SED was found to be in the
EHBL category. Similarly, an outburst in the VHE was observed by FACT during
August 2016 and once again the emission was EHBL-like.

Using the photohadronic model we have attempted to explain the VHE emissions
of the 20th of December 1995 observed by Whipple and the 7th of December 2007
observed by VERITAS respectively. As the VHE spectra of MJD 57611 and MJD
57612 are EHBL-like, we have used the two-zone photohadronic model to explain
these spectra. However, we observe that the low energy regime i.e. zone-1 of
these spectra are absent in these fittings and only zone-2 with
$\delta_{2}=3.1$ is relevant to fit the observed spectra. The photohadronic
fit to the VHE spectrum of MJD 57612 is also compared with the leptonic and
the proton-synchrotron models. To our knowledge, the VHE flaring events of MJD
57611 and MJD 57612 seem like new type of events originating from a transient
EHBL-like source and their observed VHE spectra can be deduced from Eq. (4) by
taking $\delta\gtrsim 3.1$. So far this type of behavior has not been observed
in any of the EHBL–like events from nearby objects. We believe that these
events are of the new type. However, in order to strengthen our claim, more
data on simultaneous observations on other potential sources are needed in the
future.

## Acknowledgements

We thank the referee for suggetions and constructive remarks to considerably
improve the manuscript. The work of S.S. is partially supported by DGAPA-UNAM
(Mexico) Project No. IN103522. Partial support from CSU-Long Beach is
gratefully acknowledged.

## Availability of data

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

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010, ApJ, 716, 30, doi: 10.1088/0004-637X/716/1/30
  * Acciari et al. (2011) Acciari, V. A., Aliu, E., Arlen, T., et al. 2011, ApJ, 738, 169, doi: 10.1088/0004-637X/738/2/169
  * Acciari et al. (2020) Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2020, ApJS, 247, 16, doi: 10.3847/1538-4365/ab5b98
  * Aharonian (2007) Aharonian, F. 2007, Astron. Astrophys., 470, 475.  https://arxiv.org/abs/0705.2946
  * Aharonian et al. (2007a) Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2007a, A&A, 473, L25.  https://arxiv.org/abs/0708.3021
  * Aharonian et al. (2007b) —. 2007b, A&A, 475, L9.  https://arxiv.org/abs/0709.4584
  * Ahnen et al. (2017) Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2017, A&A, 603, A31, doi: 10.1051/0004-6361/201629540
  * Ahnen et al. (2018) —. 2018, A&A, 620, A181, doi: 10.1051/0004-6361/201833704
  * Albert et al. (2007a) Albert, J., Aliu, E., Anderhub, H., et al. 2007a, ApJ, 669, 862, doi: 10.1086/521382
  * Albert et al. (2007b) —. 2007b, ApJ, 662, 892.  https://arxiv.org/abs/astro-ph/0612383
  * Aleksić et al. (2013) Aleksić, J., Antonelli, L. A., Antoranz, P., et al. 2013, A&A, 556, A67.  https://arxiv.org/abs/1211.2608
  * Aleksić et al. (2015a) Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2015a, A&A, 578, A22, doi: 10.1051/0004-6361/201424811
  * Aleksić et al. (2015b) —. 2015b, A&A, 573, A50, doi: 10.1051/0004-6361/201322906
  * Allen et al. (2017) Allen, C., Archambault, S., Archer, A., et al. 2017, MNRAS, 471, 2117, doi: 10.1093/mnras/stx1756
  * Barth et al. (2003) Barth, A. J., Ho, L. C., & Sargent, W. L. W. 2003, ApJ, 583, 134, doi: 10.1086/345083
  * Błażejowski et al. (2000) Błażejowski, M., Sikora, M., Moderski, R., & Madejski, G. M. 2000, ApJ, 545, 107, doi: 10.1086/317791
  * Catanese et al. (1998) Catanese, M., Akerlof, C. W., Badran, H. M., et al. 1998, ApJ, 501, 616, doi: 10.1086/305857
  * Cerruti (2020) Cerruti, M. 2020, Galaxies, 8, 72, doi: 10.3390/galaxies8040072
  * Cerruti et al. (2015) Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, MNRAS, 448, 910, doi: 10.1093/mnras/stu2691
  * Chandra et al. (2021) Chandra, S., Boettcher, M., Goswami, P., et al. 2021, ApJ, 918, 67, doi: 10.3847/1538-4357/ac01d1
  * Costamante et al. (2018) Costamante, L., Bonnoli, G., Tavecchio, F., et al. 2018, MNRAS, 477, 4257, doi: 10.1093/mnras/sty857
  * Costamante et al. (2001) Costamante, L., Ghisellini, G., Giommi, P., et al. 2001, A&A, 371, 512, doi: 10.1051/0004-6361:20010412
  * Dermer & Schlickeiser (1993) Dermer, C. D., & Schlickeiser, R. 1993, ApJ, 416, 458, doi: 10.1086/173251
  * Domínguez et al. (2013) Domínguez, A., Finke, J. D., Prada, F., et al. 2013, ApJ, 770, 77, doi: 10.1088/0004-637X/770/1/77
  * Domínguez et al. (2011) Domínguez, A., Primack, J. R., Rosario, D. J., et al. 2011, MNRAS, 410, 2556, doi: 10.1111/j.1365-2966.2010.17631.x
  * Elvis et al. (1992) Elvis, M., Plummer, D., Schachter, J., & Fabbiano, G. 1992, ApJS, 80, 257, doi: 10.1086/191665
  * Foffano et al. (2019) Foffano, L., Prandini, E., Franceschini, A., & Paiano, S. 2019, MNRAS, 486, 1741, doi: 10.1093/mnras/stz812
  * Franceschini et al. (2008) Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837, doi: 10.1051/0004-6361:200809691
  * Gao et al. (2013) Gao, H., Lei, W.-H., Wu, X.-F., & Zhang, B. 2013, MNRAS, 435, 2520, doi: 10.1093/mnras/stt1461
  * Gilmore et al. (2012) Gilmore, R. C., Somerville, R. S., Primack, J. R., & Domínguez, A. 2012, MNRAS, 422, 3189.  https://arxiv.org/abs/1104.0671
  * Giommi et al. (2000) Giommi, P., Padovani, P., & Perlman, E. 2000, MNRAS, 317, 743.  https://arxiv.org/abs/astro-ph/9907377
  * Godambe et al. (2007) Godambe, S. V., Rannot, R. C., Baliyan, K. S., et al. 2007, Journal of Physics G Nuclear Physics, 34, 1683, doi: 10.1088/0954-3899/34/7/009
  * MAGIC Collaboration et al. (2020a) MAGIC Collaboration, Acciari, V. A., Ansoldi, S., et al. 2020a, MNRAS, 496, 3912.  https://arxiv.org/abs/2006.06796
  * MAGIC Collaboration et al. (2020b) —. 2020b, A&A, 638, A14, doi: 10.1051/0004-6361/201935450
  * MAGIC Collaboration et al. (2020c) —. 2020c, A&A, 640, A132.  https://arxiv.org/abs/2006.04493
  * Maraschi et al. (1992) Maraschi, L., Ghisellini, G., & Celotti, A. 1992, ApJ, 397, L5, doi: 10.1086/186531
  * Perlman et al. (1996) Perlman, E. S., Stocke, J. T., Schachter, J. F., et al. 1996, ApJS, 104, 251, doi: 10.1086/192300
  * Sahu (2019) Sahu, S. 2019, Rev. Mex. Fis., 65, 307, doi: 10.31349/RevMexFis.65.307
  * Sahu et al. (2020a) Sahu, S., López Fortín, C. E., Castañeda Hernández, L. H., Nagataki, S., & Rajpoot, S. 2020a, ApJ, 901, 132, doi: 10.3847/1538-4357/abb089
  * Sahu et al. (2021a) Sahu, S., López Fortín, C. E., Castañeda Hernández, L. H., & Rajpoot, S. 2021a, ApJ, 906, 91, doi: 10.3847/1538-4357/abc9c6
  * Sahu et al. (2020b) Sahu, S., López Fortín, C. E., Iglesias Martínez, M. E., Nagataki, S., & Fernández de Córdoba, P. 2020b, MNRAS, 492, 2261, doi: 10.1093/mnras/staa023
  * Sahu et al. (2019) Sahu, S., López Fortín, C. E., & Nagataki, S. 2019, ApJ, 884, L17, doi: 10.3847/2041-8213/ab43c7
  * Sahu et al. (2021b) Sahu, S., López Fortín, C. E., Valadez Polanco, I. A., & Rajpoot, S. 2021b, ApJ, 914, 120, doi: 10.3847/1538-4357/abfd9a
  * Salamon & Stecker (1998) Salamon, M. H., & Stecker, F. W. 1998, ApJ, 493, 547.  https://arxiv.org/abs/astro-ph/9704166
  * Schroedter et al. (2005) Schroedter, M., Badran, H. M., Buckley, J. H., et al. 2005, ApJ, 634, 947, doi: 10.1086/496968
  * Shah et al. (2021) Shah, Z., Ezhikode, S. H., Misra, R., & Rajalakshmi, T. R. 2021, MNRAS, 504, 5485, doi: 10.1093/mnras/stab1244
  * Sikora et al. (1994) Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153, doi: 10.1086/173633
  * Zech & Lemoine (2021) Zech, A., & Lemoine, M. 2021, A&A, 654, A96, doi: 10.1051/0004-6361/202141062