# Photohadronic Modelling of the 2010 Gamma-ray Flare from Mrk 421

Alberto Rosales de León, Anthony M. Brown and Paula M. Chadwick
1Centre for Advanced Instrumentation (CfAI), Department of Physics, University
of Durham, South Road, Durham DH1 3LE, UK Contact e-mail: alberto.rosales-de-
<EMAIL_ADDRESS>
(Published 2020 December 12. Accepted 2020 December 9. Received 2020 December
1; in original form 2020 August 3)

###### Abstract

Blazars are a subclass of active galactic nuclei (AGN) that have a
relativistic jet with a small viewing angle towards the observer. Recent
results based on hadronic scenarios have motivated an ongoing discussion of
how a blazar can produce high energy neutrinos during a flaring state and
which scenario can successfully describe the observed gamma-ray behaviour.
Markarian 421 is one of the closest and brightest objects in the extragalactic
gamma-ray sky and showed flaring activity over a 14-day period in March 2010.
In this work, we describe the performed analysis of Fermi-LAT data from the
source focused on the MeV range (100 MeV - 1 GeV), and study the possibility
of a contribution coming from the $p\gamma$ interactions between protons and
MeV SSC target photons to fit the very high energy (VHE) gamma-ray emission.
The fit results were compared with two leptonic models (one-zone and two-zone)
using the Akaike Information Criteria (AIC) test, which evaluates goodness-of-
fit alongside the simplicity of the model. In all cases the photohadronic
model was favoured as a better fit description in comparison to the one-zone
leptonic model, and with respect to the two-zone model in the majority of
cases. Our results show the potential of a photohadronic contribution to a
lepto-hadronic origin of gamma-ray flux of blazars. Future gamma-ray
observations above tens of TeV and below 100 MeV in energy will be crucial to
test and discriminate between models.

###### keywords:

gamma-rays: galaxies – galaxies: active – BL Lacerate objects: individual: Mrk
421 – galaxies: jets.

††pubyear: 2020††pagerange: Photohadronic Modelling of the 2010 Gamma-ray
Flare from Mrk 421–References

##  1 Introduction

Blazars are a sub-class of radio-loud Active Galactic Nuclei (AGN) with a
relativistic jet pointing close to our line of sight. These objects have a
highly variable spectrum and flaring states, periods of enhanced activity on
time scales that can go from hours to months. This variability makes it
extremely difficult to model the broad-band spectral energy distribution (SED)
of these sources.

Leptonic models have been used for many years to fit the SEDs from blazars
(see, for example, Bloom & Marscher (1996); Tavecchio et al. (1998);
Sahayanathan & Godambe (2012); Aleksić et al. (2015b)). These models consider
a population of relativistic electrons to be responsible for the
characteristic two peaked SED of a blazar. In this approach, the first peak
(covering radio to X-rays) can be explained by synchrotron emission;
meanwhile, the second peak (X-rays and gamma rays) may be composed of
different contributions coming from inverse Compton interactions between the
electrons and a photon field. One option is synchrotron-self Compton
contributions (SSC) coming from the emission region inside the jet (Maraschi
et al., 1992; Dermer & Schlickeiser, 1993; Sikora et al., 1994; Bloom &
Marscher, 1996; Tavecchio et al., 1998); another option is to consider
external Compton-scattering contributions with the target photons coming from
the broad-line region (BLR), the accretion disk, the dusty torus or any other
external radiation field (Böttcher et al., 2013; Finke, 2018).

The different versions of the leptonic scenarios mentioned above have been
applied to the spectra of several objects, but the challenging observations
raise the question of whether a hadronic component is necessary to explain the
full spectrum. This idea has once again caught the interest of the scientific
community after a possible correlation between the blazar TXS 0506+056 and a
neutrino alert (IC-170922A). On September 22nd 2017, an extensive multi-
wavelength campaign was triggered by the high energy neutrino-induced muon
track event. The blazar TXS 0506+056 was reported to be $0.1^{\circ}$ from the
best-fitting neutrino direction by the Fermi Large Area Telescope (Fermi-LAT)
Collaboration (Tanaka et al., 2017). This source was in a flaring state at the
time and had a considerably brightened in the GeV band since April 2017. After
the follow-up observations, models associating neutrino and gamma-ray
production during the flaring state of the source found the gamma-ray emission
was correlated with the neutrino alert at a statistical significance of
$3\sigma$ (IceCube Collaboration et al., 2018b). In addition, the IceCube
Collaboration performed an independent analysis using prior data. The result
was an excess of high energy neutrino events coming from the direction of the
source, with respect to atmospheric backgrounds, at a significance level of
$\sim 3.5\sigma$ between September 2014 and March 2015 (IceCube Collaboration
et al., 2018a). These findings motivated an ongoing discussion of how a blazar
can reproduce the experimental data and which model (or models) can
successfully describe the observed behaviour (e.g. Righi et al. (2019); Murase
et al. (2018); Palladino et al. (2019); Padovani et al. (2019); Rodrigues et
al. (2019); Cerruti et al. (2019); Petropoulou et al. (2020); Winter & Gao
(2019); Liu et al. (2019); Halzen et al. (2019)).

The possible neutrino/gamma-ray connection exhibited between the IC-170922A
alert and TXS 0506+056 is one motivation to explore hadronic contributions,
but there are other hints and plausible evidence for this scenario. The blazar
3HSP J095507.9+355101 was recently reported to be in a flaring state and
$0.62^{\circ}$ away from the best-fit position of the alert IceCube200107A, a
muon track event (Giommi et al., 2020). There is also PKS B1424-418, a Flat
Spectrum Radio Quasar (FSRQ) which was in temporal and positional coincidence
with a high-energy starting event (HESE) on December 2012, a cascade-like
event with a reconstructed energy of 2 PeV but a median positional uncertainty
of $\sim 16^{\circ}$ that gives an estimated $\sim 5\%$ chance of coincidence
(Kadler et al., 2016). In addition, there are several studies searching for
gamma-ray counterparts and predicting consistent limits within the IC neutrino
flux so far (Krauß et al., 2014, 2015; Brown et al., 2015; Padovani et al.,
2016; Glüsenkamp, 2016); the successful SED modelling of blazars during
flaring episodes using hadronic models (Mücke & Protheroe, 2001; Diltz et al.,
2015; Diltz & Böttcher, 2016; Sahu et al., 2017, 2018a; Sahu et al., 2018b);
and hadronic emission has been proposed as an explanation for the spectral
hardening in TeV energy gamma-ray spectra, behaviour which has been observed
in some blazars (e.g. W Comae, 3C 66A (Böttcher et al., 2013), 1ES 0229+200
(Tavecchio et al., 2009), 1ES 1101-232 and H 2356-309 (Aharonian et al.,
2006)).

Hadronic scenarios propose that protons are accelerated to relativistic
energies in blazar jets. A group of models invoke photohadronic ($p\gamma$)
interactions, which involve collisions between the high energy protons and a
target photon field. Another option is hadronuclear interactions ($pp$), where
a matter target such as a gas cloud is required. Both scenarios lead to photo-
meson production, from which gamma rays and neutrinos are generated in the
decay process. The decay products will also emit radiation, including proton-
synchrotron emission, photo-pion production, electron-positron synchrotron
triggered pair cascades, or even synchrotron from the charged decay products
(muons and pions). Depending on the hadronic model, the physical conditions
and chosen parameters, there might be a dominant component, for instance
proton-synchrotron radiation (Mücke et al., 2003; Dimitrakoudis et al., 2012)
or photo-pion production (Mannheim & Biermann, 1992; Mannheim, 1993).

In the hadronuclear models, $pp$ interactions can occur, if the high energy
protons accelerated along the jet reach a matter target, for example the gas
clouds in the broad line region (BLR) around the AGN (Dar & Laor, 1997; Araudo
et al., 2010; Liu et al., 2019). This process has a lower interaction cross
section compared to the prominent $\Delta-$resonance of the photohadronic
process, so a high density target is required to improve the efficiency of the
hadronic interactions.

If a model considers mixed contributions from hadronic and leptonic origin,
then it can also be referred as lepto-hadronic (including some of the
previously given examples: Araudo et al. (2010); Diltz & Böttcher (2016);
Rodrigues et al. (2019); Cerruti et al. (2019); Petropoulou et al. (2020); Liu
et al. (2019); among others).

In this work, we test the potential of a dominant photohadronic contribution
within a lepto-hadronic scenario to fit the very high energy (VHE) gamma-ray
observations (E > 100 GeV) during a flaring blazar state. We account for the
$p\gamma$ interactions through the $\Delta-$resonance approximation; the
subsequent decay products include gamma rays and neutrinos in the following
way:

$p+\gamma\rightarrow\Delta^{+}\rightarrow\bigg{\\{}\begin{array}[]{c}p\pi^{0},\;\pi^{0}\rightarrow\gamma\gamma\\\
n\pi^{+},\;\pi^{+}\rightarrow\mu^{+}\nu_{\mu},\;\mu^{+}\rightarrow
e^{+}+\nu_{e}+\bar{\nu_{\mu}}\,.\\\ \end{array}$ (1)

The $\Delta^{+}$ particle decays into $(p+\pi^{0})$ in 2/3 of all cases while
goes to $(n+\pi^{+})$ in 1/3 of all cases (Hümmer et al., 2010). The model
considered is described in Section 5 and a broader discussion of the
$\Delta-$resonance can be found in Mücke et al. (1999); Gaisser et al. (1995).

The blazar Markarian 421 (Mrk 421; RA=66.114∘, Dec=38.209∘, z=0.031) is one of
the closest and brightest objects in the extragalactic VHE sky. It was the
first extragalactic source detected using Imaging Atmospheric Cherenkov
Telescopes (IACTs) (Punch et al., 1992) and has been regularly monitored since
then. Mrk 421 has been measured during flaring states on several occasions
(e.g. Błażejowski et al. (2005); Abdo et al. (2011); Aleksić et al. (2015a,
b)), and the recorded multi-wavelength (MWL) data from radio to high energy
gamma rays makes it an ideal candidate to test different production mechanisms
and their evolution during a flare. The 2010 flaring activity from Mrk 421
(Aleksić et al., 2015b) provides a rich dataset. We performed an analysis in
the MeV energy range (100 MeV-1 GeV) with the updated instrument response
functions (IRFs) from the Fermi-LAT to obtain an input seed photon spectrum to
a photohadronic model that could provide a good fit to the VHE gamma-ray data.

This paper is structured as follows: in Section 2 we review the 2010 flare, in
Section 3 we describe the Fermi-LAT analysis undertaken, and in Section 4 we
describe the results of that analysis. Section 5 describes the photohadronic
model and the method used to fit this to the data, and Section 6 describes the
results of this fit. Finally, we discuss our results in Section 7 and conclude
in Section 8.

##  2 Flaring Activity in 2010

Mrk 421 exhibited flaring activity over a 14-day period in 2010 from March 10
to March 22 (MJD 55264-55277). At the time, a multi-instrument campaign was
performed which included the gamma-ray space telescope Fermi-LAT and three
IACTs: the Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) telescope
system, the Very Energetic Radiation Imaging Telescope Array System (VERITAS)
and the Whipple gamma-ray telescope.

MAGIC took 11 observations in stereoscopic mode with exposure times ranging
from 10 to 80 min each, which led to 4.7 h of good-quality data with a zenith
angle range of $5^{\circ}-30^{\circ}$. The data collected were taken in dark
conditions and were not affected by moonlight, but the data recorded on MJD
55272 and 55275 suffered from bad weather and were therefore removed from the
MWL observations (Aleksić et al., 2015b). For more details on the MAGIC
telescope system see Aleksić et al. (2012).

VERITAS monitored the source on MJD 55260, 55265, and 55267-55274 with a 10
min run per day. The observations were performed at zenith angles
$18^{\circ}-23^{\circ}$ to benefit from the lowest possible energy threshold.
Further information about the VERITAS instrument can be found in Perkins et
al. (2009).

The Whipple telescope performed 10 observations in ON/OFF and TRK (tracking)
modes (Pichel, 2009), lasting from one to six hours each on MJD 55267-55271
and MJD 55273-55277. The dataset collected for this flaring period amounts to
36 h. More information and details about the Whipple telescope can be found in
Kildea et al. (2007).

The VHE gamma-ray data from ground based IACTs published in Aleksić et al.
(2015b) were used to test a possible dominant contribution from photohadronic
interactions. In this paper, we focus on the modelling of the VHE gamma-ray
observations. In this case, the Fermi analysis and the IACT data were the two
key elements needed, and data from other wavelengths were not critical for our
calculations (see Aleksić et al. (2015b) for a full description of the MWL
observations). The light curves from MAGIC, VERITAS and Whipple above 200 GeV
are shown in Figure 1. In terms of simultaneity of the data, the shorter
observation times of the IACTs are embedded in the 2-day bins period defined
for the Fermi analysis, noting that there is a 7-h time difference between the
VERITAS/Whipple and MAGIC observations due to their different longitudes. The
variability reported in the gamma-ray data corresponds to daily changes in the
VHE emission; no intra-night variability was reported on the days studied.

The photohadronic model does not aim to describe the whole SED with purely
hadronic components. This approach relies on a standard leptonic scenario to
explain the low energy peak of the SED and provide seed photons for the
$p\gamma$ interactions. To describe the blob from which the gamma-ray photons
are produced, we have adopted the values of the physical parameters in the
one-zone model of Aleksić et al. (2015b) relating to magnetic field, Doppler
factor and the radius of the emission region. These parameters are fixed
during the flaring events, and the evolution of the spectral parameters
provides snapshots on the different days considered for modelling.

Mrk 421 was highly active during other months in 2010, and the VERITAS
Collaboration reported another flare in February 2010 (MJD 55234-55240), the
brightest ever observed from this object in VHE gamma rays (Abeysekara et al.,
2020). They concluded that the time variability of the source is difficult to
explain using a single-zone SSC model. This result provides another motivation
to try to extend the current models and look for hadronic contributions.

##  3 Fermi Analysis

Launched in June 2008, the Large Area Telescope (LAT) on board the Fermi
satellite is a pair conversion telescope covering the energy range from 20 MeV
to more than 300 GeV. The LAT’s combination of a wide field of view (FoV
$\simeq 2.4$ sr) and large effective area allows it to scan the entire gamma-
ray sky approximately every 3 hours (Atwood et al., 2009). The data analysis
described here was performed with the Pass8v6 version of the IRF and the
v11r5p3 Science Tools software with Fermipy (Wood et al., 2017). This IRF
provides a full reprocessing of the entire mission dataset, including improved
event reconstruction, a wider energy range, better energy measurements, and
significantly increased effective area in comparison to previous versions. As
a result, we were able to extend our analysis into a lower energy range than
previously possible, resulting in a more comprehensive spectrum of the source
during the flaring period.

The Fermipy package provides a set of tools and an interface (GTAnalysis) to
perform the data preparation, modelling, statistics and analysis tasks. The
4FGL-DR2 catalog (Abdollahi et al., 2020), containing the positions and
spectral information of the known gamma-ray sources, was used during the
analysis.

Our data reduction steps considered all ‘SOURCE’111See
https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data/LAT_DP.html
class events photons in an energy range of 100 MeV to 1 GeV between 2010 March
9th to 22nd (MJD 55264-55277). The flaring period studied was divided into
shorter 2-day intervals; this time period was the minimum to obtain enough
photon events to calculate the SED points and upper limits. The Fermi-LAT data
are publicly-available and were downloaded from the LAT data server
system222See https://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi.

The analysis was performed inside a 15∘ region of interest (RoI) around Mrk
421’s coordinates, a 90∘ zenith cut angle was applied to avoid gamma rays
coming from Earth’s atmosphere. The gamma-ray data were then binned using 5
bins per decade in energy and 0.1∘ spatial bin size. To remove sub-optimal
data, only the events within good time intervals (GTI) were analysed, these
were selected by the ‘gtmktime’ tool filters ‘(DATA_QUAL>0) &&
(LAT_CONFIG==1)’ and a cut above 52∘ in rocking angle.

A model consisting of gamma-ray point sources and a background with a Galactic
diffuse and an extragalactic component was employed. All the sources listed in
the 4FGL catalog inside a 20∘ neighbourhood from the centre of the ROI were
included; this was to account for the possible contributions of sources near
the edge of our ROI. The spatial model, position and spectral parameters of
the sources were adopted from the 4FGL catalog. The Galactic diffuse emission
component used in the model was ‘gll_iem_v06.fits’333 See
https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html. The
extragalactic gamma-ray contribution coming from unresolved extragalactic
sources, and residual (misclassified) cosmic-ray emission was also included as
an isotropic spectral template parametrised in the file
‘iso_P8R2_SOURCE_V6_v06.txt’3. In the 4FGL catalog, Mrk 421 (4FGL
J1104.4+3812) is listed as a source with a log-parabola (LP) spectrum type.
However, over the short time intervals which we are considering, there were
not sufficient photon statistics to allow a log-parabola model to be
distinguished from a power-law (PL) model, so the latter was adopted. The
spectral models should not differ significantly in the range of interest for
our analysis and the SED extrapolation below 100 MeV.

A maximum likelihood method444The maximum-likelihood test statistics (TS) is
defined as $TS=2[logL-logL_{0}]$ where $L$ and $L_{0}$ are the likelihood when
the source is included or not, respectively Mattox et al. (1996). was used for
fitting the ROI; the spectral shape parameters of the sources were left free
to vary within a 5∘ radius around the ROI’s centre. The two background
components were also left free to vary during the maximum likelihood fitting.
The sources in the model that were considered insignificant (TS < 1) were
discarded. A second optimisation and fit steps were applied to the ROI using
the new model with the same free components and criteria used before.

To check for any point sources inside the ROI which are not listed in the 4FGL
catalog, the ‘find_sources’ routine built within Fermipy was implemented. No
significant extra source candidates (TS>16) were found during the analysis.

The SED points for each 2-day bin were calculated using the ‘sed()’ tool
included in Fermipy, which performs an independent maximum likelihood fit per
energy bin for the flux normalisation. In our analysis, the energy range
studied (100 MeV to 1 GeV) was split into 5 evenly spaced log-energy bins. The
same analysis steps described above were applied to each 2-day interval in the
flaring period studied (MJD 55264-55277) to obtain the corresponding SED and
PL spectral parameters.

Figure 1: Light curve of Markarian 421 during the 14 day flaring period in
March 2010. The upper plot shows the gamma-ray flux in the energy range 100
MeV < $E{\gamma}$ < 1 GeV and the points are calculated in couple of days
bins. In the bins with a TS<25 upper limits for the flux are shown. The lower
plot presents the light curves for MAGIC, VERITAS and Whipple above 200 GeV.

##  4 Fermi Results

Time | TS | Flux | $N_{\rm PL}$ | $\kappa$
---|---|---|---|---
MJD |  | [$10^{-7}$ ph cm-2 s-1] | [$10^{-11}$ MeV cm-2 s-1] |
55266-67 | 29 | 2.78 $\pm$ 0.98 | 2.11 $\pm$ 1.32 | 2.21 $\pm$ 0.44
55274-75 | 42 | 0.91 $\pm$ 0.33 | 4.51 $\pm$ 2.60 | 1.17 $\pm$ 0.61
55276-77 | 47 | 1.76 $\pm$ 0.50 | 5.94 $\pm$ 3.12 | 1.11 $\pm$ 0.61
55264-65 | 10 | < 2.34 | - | -
55268-69 | 6 | < 2.07 | - | -
55270-71 | 18 | < 2.27 | - | -
55272-73 | 19 | < 3.89 | - | -
Table 1: Summary table of the spectral parameters of Mrk 421. The 2-day bins
with significant TS values are listed on the top of the table. The 4th and 5th
columns correspond to the optimised parameters for a PL fit coming from the
Fermi analysis performed. The days with a low TS value are listed in the
bottom of the table together with the upper limits for the gamma-ray flux.

The Fermi analysis was done with the purpose of characterising the seed photon
spectrum of the source using an extrapolation below the 100 MeV energy range.
In the assumed scenario, the seed photons for the $p\gamma$ interactions are
expected to be between 2 and 168 MeV (see Section 5.1). For the Fermi
analysis, we selected a photon energy range from 100 MeV (the recommended
starting energy) and up to 1 GeV, then extrapolate to the lower energy range
of interest. Using the updated version of the IRFs, the shortest time bins
which allowed the spectrum to be obtained were 2-days in length.

As discussed in the previous section, the spectrum in the MeV range was
characterised using a PL model:

$\frac{\rm{d}\it{N}}{\rm{d}\epsilon_{\gamma}}=N_{\rm
PL}\epsilon_{\gamma}^{-\kappa},$ (2)

where the normalisation constant $N_{\rm PL}$ and the spectral index $\kappa$
act as free parameters that were optimised to get the best-fit values. The PL
description was then used to extend the spectrum below 100 MeV. The spectra of
the source alongside the PL extrapolation are shown in Figure 2 for the
selected bins.

The spectral parameters of the days on which the source was detected
significantly (TS>25) are shown in Table 1 and the corresponding light curves
from our Fermi analysis and in the VHE energy range are shown in Figure 1.
This analysis extends the previous results by Aleksić et al. (2015b), which
started at 300 MeV rather than 100 MeV. There are no significant flux changes
in the VHE band in the combined 2-day bins which we used for the Fermi
analysis and subsequent modelling. The remaining days presented low photon
statistics and were not considered for further VHE fitting with the
photohadronic model. In order to get our final result, the gamma-ray spectrum
of the source was studied in the MeV energy range, with the fitted spectrum
serving as an input for the photohadronic modelling. The time bins in our
analysis and that of Aleksić et al. (2015b) coincide; however, the extended
spectrum analysed with the updated IRF allowed us to calculate a PL
extrapolation in our range of interest.

Once the MeV region is characterised, the input seed photon spectrum for the
photohadronic model can be expressed as

$\Phi_{\rm
input}=\epsilon_{\gamma}^{2}\frac{\rm{d}\it{N}}{\rm{d}\epsilon_{\gamma}}=N_{\rm
PL}\epsilon_{\gamma}^{-\kappa+2}.$ (3)

The uncertainty of the seed photon spectrum will impact on the optimisation
process of the other free parameters within the model ($\alpha$ and
$A_{\gamma}$, see section 5.1) and therefore the final fitting result.

(a) MJD 55266-67

(b) MJD 55274-75

(c) MJD 55276-77

Figure 2: Fermi spectra (blue points) and power-law extrapolation (magenta
line) for the MeV range in 2-days bins: (a) MJD 55266-67, (b) MJD 55274-75,
(c) MJD 55276-77. The black dotted vertical lines are positioned at 2 and 168
MeV, which is the expected energy range for the seed photons. The spectral
parameters for the selected days are summarised in Table 1.

Alongside the spectral parameters, the light curve of Mrk 421 during the
flaring period (MJD 55264-55277) was calculated using the photons in the MeV
energy range (100 MeV to 1 GeV), also using 2-day temporal bins (see top of
Figure 1). The data reduction steps described in the previous section were
followed to prepare, optimise and fit the ROI.

The light curve was generated using the ‘gta.lightcurve’ method within
Fermipy, which performs a likelihood fit for each time bin. An optimised
region of 15 degrees was considered in the energy range of 100 MeV to 1 GeV
with 5 bins per decade in energy. The average photon flux for Mrk 421 was
$1.67\pm 0.28\times 10^{-7}$ ph cm-2 s-1. For each significant bin (TS>25),
the gamma-ray flux is reported in the third column of Table 1, in bins where
the TS value was below 25 we present upper limits.

##  5 Photohadronic Contributions to the Flares

In the lepto-hadronic scenario we consider, a one-zone leptonic model is
assumed to contribute to the SED via electron synchrotron and SSC radiation;
this is a standard leptonic interpretation. The low-energy peak from the SED
comes from the synchrotron component of leptonic origin, while the SSC
component is assumed to provide the target photon field in the MeV range. For
the photohadronic contribution to arise, it is hypothesised that protons are
accelerated into the single spherical emission region of radius
$R^{\prime}_{\rm f}$ (flaring blob), with a tangled magnetic field $B$,
propagating along the jet with a velocity $\beta_{\Gamma}c$ and an associated
bulk Lorentz factor $\Gamma$. The jet forms a small angle $\theta$ with
respect to the line of sight, which results in a Doppler boosting
characterised by the Doppler factor ${\cal D}$. It is proposed that during the
flaring episode, the blazar possesses a dense, compact inner jet structure
(Ghisellini et al., 2005; Marscher et al., 2008, 2010; MacDonald et al., 2015;
Homan et al., 2015; Walker et al., 2016). Geometrically this represents a
double conical shape, with a compact and smaller region enclosed by the jet
along its axis (for a schematic view see Figure 1 in Sahu et al. (2016)). The
inner compact region has a photon density $n^{\prime}_{\gamma,\rm{f}}$, which
is much higher than the outer region $n^{\prime}_{\gamma}$, this helps to
increase the efficiency of the photohadronic interactions. The prime notation
is adopted to refer the jet comoving reference frame.Inside the emission
region an electron population will produce synchrotron and SSC radiation
following the usual one-zone leptonic scenario. $p\gamma$ interactions can
emerge from the collisions between high energy protons and the internal photon
field, the SSC photons in the MeV range will serve as targets for the
interaction to get a TeV energy gamma rays from the photo-pion production. A
more comprehensive review of the photohadronic flaring model can be found in
Sahu et al. (2012, 2013).

In this framework, the $\Delta-$resonance approximation is used. The
$\Delta^{+}$ particle has a mass of $m_{\Delta}=1.232$ GeV; this is the
threshold for interaction and corresponds to the production of the particle at
rest. Above this energy threshold, the cross-section of the process is
enhanced and this decay channel becomes dominant over other components. At its
peak, the cross-section of the $\Delta-$resonance reaches a value of
$\sigma_{\rm{peak}}\sim 500$ $\mu\text{barn}$, which is $\sim 5\times
10^{-28}\text{cm}^{2}$, being this bigger by a factor of $\sim 5$ than the
direct channel cross section production (Hümmer et al., 2010).

The threshold of the interaction dictates an energy relation between the
proton energy $E^{\prime}_{\rm p}$ and seed photon
$\epsilon^{\prime}_{\gamma}$ in the emission region reference frame:

$E^{\prime}_{\rm p}=\frac{m^{2}_{\Delta}-m^{2}_{\rm
p}}{2\epsilon^{\prime}_{\gamma}(1-\rm{cos}\phi)}\,,$ (4)

where $m_{\rm p}$ is the proton mass and $\phi$ is the angle formed between
the interacting particles. Since the proton will collide with the target
photons from all directions there is not a preferred angle of interaction and
$1-\rm{Cos}(\phi)\sim 1$. In the observer’s frame, due to the Doppler boosting
effect from the jet, the proton energy will be enhanced as:

$E_{\rm p}=\frac{\Gamma}{1+z}E^{\prime}_{\rm p}\,,$ (5)

where $E_{\rm p}$ is the energy which would be measured by the observer if the
proton could be able to escape the source and reach Earth without energy loss.
In a similar way, the target photon energy in the observer’s frame can be
expressed as:

$\epsilon_{\gamma}=\frac{D}{1+z}\epsilon^{\prime}_{\gamma}\,.$ (6)

Considering that each pion carries $\sim 20\%$ of the proton energy (Hümmer et
al., 2010) and in the photo-pion production 2 gamma rays are produced from the
$\pi^{0}$-decay, we have the following relation between the gamma-ray photon
energy $E_{\gamma}$ produced with a proton energy $E_{\rm p}$ in the
observer’s frame:

$E_{\gamma}=\frac{1}{10}\frac{D}{1+z}E^{\prime}_{\rm
p}=\frac{D}{10\Gamma}E_{\rm p}\,.$ (7)

###  5.1 Photohadronic Fit

In this scenario, high energy protons can be injected into a confined region
(a spherical blob) of radius $R^{\prime}_{\rm f}$ inside the blazar’s jet. For
simplicity, it is assumed that the internal jet region and the external jet
are moving with almost the same bulk Lorentz factor $\Gamma$. For blazars, we
consider that the Doppler factor and the Lorentz factor are approximately of
the same magnitude ${\cal D}\sim\Gamma$ (Oikonomou et al., 2019).

From the energy threshold condition to produce the $\Delta-$resonance, an
energy relation between the target photons $\epsilon_{\gamma}$ and the gamma-
ray photons $E_{\gamma}$ in the observer’s frame can be expressed as:

$E_{\gamma}\epsilon_{\gamma}\simeq\text{0.032}\frac{{\cal
D}^{2}}{(1+z)^{2}}\;\text{GeV}^{2}\,.$ (8)

The central region of an AGN possesses shocks that are able to accelerate
electrons and ions trough the Fermi mechanism; when one of these relativistic
particles crosses the shock from downstream to upstream or vice versa, it
gains energy (Baring, 1997). A PL injected spectrum for the protons is
considered: ${\rm d}N(E_{\rm p})/{\rm d}E_{\rm p}\propto E_{\rm p}^{-\alpha}$,
where the spectral index $\alpha$ is treated as a free parameter in the model.
The high energy protons will interact in the inner jet region where the seed
photon density is $n^{\prime}_{\gamma,\rm f}$. The gamma-ray spectrum obtained
at VHE will depend proportionally on the photon background and the injected
proton spectrum (Sahu et al., 2012, 2013):

$F_{\rm int}(E_{\gamma})\propto n^{\prime}_{\gamma,\rm f}E^{2}_{\rm
p}\frac{\rm{d}\it{N}_{\rm p}}{\rm{d}\it{E}_{\rm p}}\,.$ (9)

The seed photon density will impact on the efficiency of the $p\gamma$
process; a low value reduces the chances of interaction and therefore the
gamma-ray photon emission obtained by this method. This photon density in the
inner region of the jet is unknown, but we can set a very rough upper limit by
assuming that the Eddington luminosity ($L_{\rm Edd}$) of the source should
not be exceeded and that it is equally shared by the jet and the counter jet.
The upper limit on the seed photon density can be placed using:

$n^{\prime}_{\gamma,{\rm f}}\ll\frac{L_{\rm Edd}}{8\pi R^{\prime 2}_{\rm
f}\epsilon^{\prime}_{\gamma}}\,,$ (10)

where $L_{\rm Edd}$ for Mrk 421 is $\sim 2.5\times 10^{46}$ erg s-1 for a
black hole mass of $2\times 10^{8}M_{\odot}$ (Sahu et al., 2016). This gives
us a limit of $\sim 3\times 10^{16}$ ph cm-3. The photon density
$n^{\prime}_{\gamma,\rm f}$ is proportional to the luminosity
$L_{\gamma}(\epsilon_{\gamma})$, and inversely proportional to the seed photon
energy $\epsilon_{\gamma}$. The luminosity at a certain energy is proportional
to the observed flux $\Phi_{\rm input}(\epsilon_{\gamma})$, which is known
from the PL extrapolation obtained in Section 4, so we have that:

$n^{\prime}_{\gamma,\rm f}\propto\Phi_{\rm
input}(\epsilon_{\gamma})\epsilon_{\gamma}^{-1}\,.$ (11)

This means the intrinsic gamma-ray flux will follow:

$F_{\rm int}(E_{\gamma})\propto\Phi_{\rm
input}(\epsilon_{\gamma})\epsilon_{\gamma}^{-1}E_{\rm
p}^{2}\frac{\rm{d}\it{N}_{\rm p}}{\rm{d}\it{E}_{\rm p}}\,,$ (12)

considering the PL injected spectrum of the protons and using the energy
relations between the proton energy ($E_{\rm p}$), the seed photon energy
($\epsilon_{\gamma}$) and the energy of the gamma-ray photon ($E_{\gamma}$),
the intrinsic gamma-ray flux $F_{\rm int}$ coming from the $\pi^{0}$-decay can
be expressed as:

$F_{\rm int}(E_{\gamma})=A_{\gamma}\Phi_{\rm
input}(\epsilon_{\gamma})\left(\frac{E_{\gamma}}{\text{TeV}}\right)^{-\alpha+3}\,,$
(13)

where $A_{\gamma}$ is a dimensionless normalisation constant that absorbs the
information from the various proportional relations given above, and $\alpha$
is the power index from the assumed proton spectrum. In this methodology,
$A_{\gamma}$ and $\alpha$ are optimised to fit the VHE gamma-ray data day by
day.

When calculating the gamma-ray spectra, we must account for the attenuation of
the high energy gamma rays due to the pair production effect with the
Extragalactic Background Light (EBL). The EBL provides an attenuation factor
of the form $e^{-\tau_{\gamma\gamma}}$, where $\tau_{\gamma\gamma}$ is known
as the optical depth, which increases at higher energies. In this case, we
apply the model of Domínguez et al. (2011). Including this attenuation factor
in the expression for the gamma-ray flux (Equation 13), we get:

$F_{\gamma}(E_{\gamma})=A_{\gamma}\Phi_{\rm
input}(\epsilon_{\gamma})\left(\frac{E_{\gamma}}{\text{TeV}}\right)^{-\alpha+3}{\rm
e}^{-\tau_{\gamma\gamma}(E_{\gamma},{\rm z})}\,.$ (14)

###  5.2 Contribution to Mrk 421’s Flare

During the 2010 flaring period of Mrk 421, the VHE data recorded by the IACTs
lies in an energy range of 80 GeV to 5 TeV (Aleksić et al., 2015b). The energy
relation from the Equation 8 indicates that the seed photon energy in the
$p\gamma$ interaction is between 2 and 168 MeV (in the observer’s reference
frame). The reduction in the source flux combined with the low sensitivity of
the LAT below 100 MeV are an impediment to obtaining precise measurements at
these energies, and therefore we use the results obtained in Section 4 to
estimate the flux coming from the seed photons ($\Phi_{\rm
input}(\epsilon_{\gamma})$). The PL input is shown on Figure 2.

The energy range of the seed photons and the threshold condition for the
$\Delta-$resonance can be used to estimate proton energy. From Equation 7, we
have that $E_{\rm p}\sim 10E_{\gamma}$ in the observer’s reference frame; if
measured from Earth, these high energy protons, boosted by the blazar’s jet,
will be detected in an energy range of $\sim$ 800 GeV to 50 TeV, which
corresponds in the emission region reference frame to $\sim$ 40 GeV <
$E^{\prime}_{\rm p}$ < 2.45 TeV. This is the energy range of the protons to
reach the threshold condition for the $\Delta-$resonance.

Time | $A_{\gamma}$ | $\alpha$ | Preferred | $\Delta{\rm AIC}_{\rm{SSC},\it{p}\gamma}$ | $\Delta{\rm AIC}_{\rm{two-zone,SSC},\it{p}\gamma}$
---|---|---|---|---|---
MJD |  |  | Model |  |
55266 | 5.02 $\pm$ 2.74 | 3.12 $\pm$ 0.07 | two-zone SSC | 25.45 | -48.78
55267 | 27.24$\pm$ 12.79 | 3.41 $\pm$ 0.09 | $p\gamma$ | 6.11 | 9.10
55274 | 0.19 $\pm$ 0.01 | 2.31 $\pm$ 0.03 | inconclusive | 2.54 | 0.73
55276 | 0.10 $\pm$ 0.02 | 2.17 $\pm$ 0.03 | $p\gamma$ | 26.40 | 2.04
55277 | 0.18 $\pm$ 0.02 | 2.32 $\pm$ 0.03 | $p\gamma$ | 5.92 | 2.90
Table 2: Summary table of the photohadronic fit for each day in Figure 3. The
optimised values for the normalisation constant $A_{\gamma}$ and the power
index $\alpha$ are shown in the second and third column. The AIC difference
between the one-zone and two-zone SSC model with respect to the photohadronic
($p\gamma$) model is shown on the fifth and sixth column respectively. An
inconclusive result is obtained, if the AIC difference between the two models
with the lowest values is less than 2.

The emission region has some physical parameters (magnetic field, Doppler
factor, radius of the spherical blob) which are the same as those used in the
calculation of the photohadronic component. These are fixed parameters taken
from the one-zone emission region leptonic model of Aleksić et al. (2015b) and
are: magnetic field of $B=38$ mG, a Doppler factor ${\cal D}=21$ and a radius
of the emission region log${}_{10}(R^{\prime}_{\rm f}[\text{cm}])=16.72$. The
variability reported in Aleksić et al. (2015b) corresponds to daily changes in
the VHE emission, a time scale which in principle is related to the proton
injection. The power index $\alpha$ and the normalisation constant
$A_{\gamma}$ were estimated daily for each VHE dataset considered (days at the
top of Table 1). These two free parameters were optimised using a $\chi^{2}$
minimisation method within Scipy Python package (Virtanen et al., 2020).

##  6 Model Fit Results

The best-fit value parameters for the photohadronic fit ($A_{\gamma}$,
$\alpha$) are given in Table 2. The photohadronic component ($p\gamma$ model)
for the best-fit values is shown as a magenta continuous curve in Figure 3;
the one-zone leptonic model from Aleksić et al. (2015b) (SSC model) is shown
as a dash-dot black line, and the two-zone model (two-zone SSC) used in the
same paper is represented by a dashed red line. The two-zone SSC model assumes
one quiescent blob producing the steady emission, and a smaller independent
blob responsible for the temporal evolution of the SED (flaring blob). The VHE
energy points in Figure 3 are a combination from the IACT observations during
the flare (MAGIC and VERITAS, see section 2).

To compare quantitatively these models, we perform an Akaike Information
Criterion (AIC) test (Akaike, 1974), which can be used to determine if a model
fit is significantly better than another (Harris et al., 2014; Bozdogan, 1987;
Forster, 2000). The AIC takes into account both the goodness of the fit, and
the simplicity of the model. This is done by assessing the likelihood and the
number of free parameters adopted, AIC is defined as:

${\rm AIC}_{\rm s}=-2ln(L_{\rm s})+2k_{\rm fs}\,,$ (15)

where $L_{\rm s}$ is the likelihood of the model used and $k_{\rm fs}$ is the
number of free parameters in the model. In our optimisation process, the
$\chi^{2}$ statistic was used to obtain our best-fit results and is adopted as
our likelihood function. Most of the parameters used in the photohadronic
model are fixed. There are 4 free parameters considered within the model,
which vary to fit the daily VHE observations: 2 from the PL input of the seed
photons and 2 from the calculation of the gamma-ray flux $F_{\gamma}$ from the
photohadronic component.

According to Aleksić et al. (2015b), 11 parameters were used in the one-zone
leptonic model, 9 of which come from a broken power-law function (with two
breaks) in the electron energy distribution required to satisfactorily
describe the SED during the flaring period. From these, 5 parameters were left
free to vary during the flaring days and the rest were fixed. For the two-zone
leptonic model, 20 parameters were used to describe the SED. After fixing the
parameters of the quiescent blob, only 4 parameters were left free to vary in
the flaring blob. We use the number of free parameters in each model to
calculate the corresponding AIC value.

The AIC is based on a theoretical framework within information theory and
provides a way to evaluate both the simplicity and accuracy of the model fits.
The difference between the AIC of two models p,q is expressed as:

$\Delta{\rm AIC}_{\rm p,q}={\rm AIC}_{\rm p}-{\rm AIC}_{\rm q}\,.$ (16)

The AIC difference enables the models considered to be compared and ranked.
The model with the lowest AIC represents the best description of the empirical
data available. Any model comparison with a $\Delta{\rm AIC}_{\rm p,q}>2$
above the minimal AIC value is considered significantly worse (Burnham &
Anderson, 2001; Lewis et al., 2011).

A direct comparison through the AIC difference test between the $p\gamma$
model and the two leptonic models considered is shown on Table 2. Among the
IACTs observations from the 2010 flaring period from which we were able to fit
a spectrum (Figure 3), we found that the photohadronic fit was the preferred
model in 3 out of 5 cases according to the AIC comparison test. On one day
(MJD 55266), the two-zone SSC model has the minimum AIC by a large difference
with respect to the other models. Meanwhile on MJD 55274, the AIC test was
inconclusive, as no difference larger than 2 was found between the $p\gamma$
and the two-zone SSC model.

In all cases, the comparison of the $p\gamma$ model with the one SSC model
results in a $\Delta{\rm AIC}_{\rm SSC,\it{p}\gamma}>2$, which means the
$p\gamma$ model is favoured over the one-zone SSC model and represents a
significantly better fit (Table 2). Nevertheless, this does not necessarily
means that the $p\gamma$ model fits all the datasets well.

During the first days of the observations (MJD 55264-67) the source was in its
maximum emission state. The gamma-ray flux detected, as shown in Table 1 for
the seed photon spectrum reaches a peak on this day, and the spectral index
associated is also the highest at $\kappa\sim 2.2$. This behaviour also
relates to the optimised free parameters ($\alpha$ and $A_{\gamma}$) which
presents their maximum values during these couple of days.

The fitting for MJD 55266 presents a large difference between the calculated
AIC values, the two-zone SSC resulting in the lowest AIC and therefore
favoured as the best fit for data from this day. The SSC model over-predicts
the VHE flux just above 100 GeV, whereas the $p\gamma$ model over-predicts the
flux above the TeV energy limit and is flatter below that point, thus
underestimating the flux peak. Neither of these two models represent a good
fit to the data, which is reflected in the AIC difference between them and the
two-zone SSC model.

On the following day (MJD 55267), using the corresponding PL input to
calculate the photohadronic contribution, the overall behaviour of the VERITAS
observations can be fitted with the $p\gamma$ model, but the PL behaviour from
the input takes over if extended below 100 GeV and predicts an increase in
flux; this is a non-physical effect that is outside our validity range and
represents a caveat of the model. Nonetheless, the $p\gamma$ model scored the
lowest AIC value among the 3 models and the difference in this day is enough
to consider it the best fit.

For the final days of the flaring period (MJD 55274 to 55277) three days of
VHE data were available, there being no VHE observations on MJD 55275. The PL
spectral parameters used as an input for the $p\gamma$ model are very similar
for these 3 days. The spectral index during this period is $\kappa\sim 1.1$
and the power index of the proton spectrum is around $\alpha\sim 2.2-2.3$.

During this later part of the flaring period, the lowest AIC values found were
from the $p\gamma$ model, indicating these are the best fits to the datasets
and ranking the model as the preferred option on MJD 55276 and 55277.

On MJD 55274 the three models predict a similar flux in the 100 GeV to 1 TeV
range, but the photohadronic contribution differs at higher energies where it
is expected a larger contribution. This is the only day on which the models
were tested with data from both MAGIC and VERITAS, although the $p\gamma$
model reached the lowest value from the three, the AIC difference with the
two-zone SSC model was not significant enough ($\Delta{\rm AIC}_{\rm
SSC,\it{p}\gamma}>2$) for the $p\gamma$ model to be selected as the preferred
model, therefore this day is listed as inconclusive in Table 2.

On MJD 55276, the photohadronic model is able to reproduce accurately the VHE
gamma-ray data and has the lowest AIC in the period studied. In contrast, the
SSC model over-predicts the MAGIC observations, and although the two-zone SSC
is a better approximation to the data, the AIC difference in favour of the
$p\gamma$ model is enough to select it as the preferred model.

On the final day studied (MJD 55277) both leptonic models are a good
approximation to the data; nevertheless, the $p\gamma$ model results in a much
lower AIC value due to the accuracy of the fit and its simplicity (in terms of
free parameters), remaining as the preferred model. The behaviour of the
$p\gamma$ model and the influence of the PL input approximation will be
discussed in the next section.

(a) MJD 55266 (b) MJD 55267

(c) MJD 55274 (d) MJD 55276

(e) MJD 55277

Figure 3: Photohadronic fit for the VHE gamma-ray data on flaring days with
significant TS values: (a) MJD 55266, (b) MJD 55267, (c) MJD 55274, (d) MJD
55276, (e) MJD 55277. The photohadronic component calculated from the PL input
is shown in magenta for the valid energy range of the model, which extends
roughly down to 80 GeV. The one-zone (two-zone) SSC model from Aleksić et al.
(2015b) is shown as a dash-dot black (dashed red) line. The calculated AIC
values for the three models are included for comparison. VHE data points are
from MAGIC and VERITAS observations.

##  7 Discussion and outlook

The fits obtained using the $p\gamma$ model over the flaring period are a
significantly better description of the VHE data than a purely one-zone
leptonic model, due to a combination of the fit quality obtained and the
simplicity of the model in terms of the approximations used and the number of
free parameters adopted. The two-zone SSC model is more competitive, being
selected as the preferred model on MJD 55266 and scoring a similar AIC value
to $p\gamma$ model on MJD 55274. It is worth noting that for the AIC
calculation we are using an optimistic estimation of the number of free
parameters in the leptonic models, following Aleksić et al. (2015b) in which
only 4 (5) from a total of 20 (11) parameters are considered as free for the
two-zone (one-zone) SSC leptonic model during the flaring.

An interesting feature is that the best fits for a photohadronic contribution
are obtained in the last 4 days of the flare. In a more complex hadronic
scenario, besides the photo-pion production, there is gamma-ray radiation
emitted by proton-synchrotron cooling and the synchrotron radiation of the
secondary charged particles. These hadronically induced gamma rays are usually
in competition with synchrotron and inverse Compton photons radiated by
primary electrons considered in the usual leptonic model (Rachen & Mészáros,
1998). On some days the hadronic component could be dominant at VHE and then
followed by a dominant SSC leptonic component. If the proton injection occurs
randomly, there is no preferred time for this to happened during the flare.
The interplay between these mechanisms could lead to a time-dependent model
with a dominant component at VHE (Diltz & Böttcher, 2016; Dimitrakoudis et
al., 2012). The current stage of the model relies on the leptonic
contributions to provide the target photon field at MeV energies (presumably
the SSC component) and ignores other hadronic components based on the dominant
process of photo-pion production through the $\Delta-$resonance. This work
aims to be a first step into a more comprehensive lepto-hadronic modelling of
flaring blazars.

A common problem with hadronic modelling is the high proton energy required to
produce the observed gamma-ray emission. However, in this case the high
frequency of the seed photos considered in the $p\gamma$ interactions lowers
the energy threshold for the accelerated protons to 40 GeV < $E^{\prime}_{\rm
p}$ < 2.45 TeV in the comoving frame (the emission region), which is below the
extreme energies considered in other hadronic models (Mannheim & Biermann,
1992; Mücke et al., 2003). This feature from the $p\gamma$ model can be
considered an advantage that would facilitate the conditions for objects like
blazars to produce VHE gamma rays and neutrinos from hadronic interactions.
Using the Hillas criterion (Hillas, 1984; Ptitsyna & Troitsky, 2010; Meli et
al., 2008) with the considered parameters of the emission region ($B$, ${\cal
D}$, $R^{\prime}_{\rm f}$), a proton could be accelerated up to
$E^{\prime}_{\rm p,max}\sim 650$ TeV at the source, more than sufficient for
the $p\gamma$ model to produce the VHE gamma-ray emission.

An interesting feature of the $p\gamma$ model (Figure 3) is the increase of
the energy flux above TeV energies. This is where the one-zone SSC model and
the $p\gamma$ model differ and needs to be tested in future observations. The
forthcoming Cherenkov Telescope Array (CTA) (CTA Consortium et al., 2019) will
be critical to differentiating between these scenarios. The expected
sensitivity and capabilities of the CTA observatories will enable detection of
the gamma-ray flux from Mrk 421 and other near blazars with unprecedented
accuracy and instantaneous sensitivity above tens of TeV in energy.

CTA will also take part in the Neutrino Target of Opportunity (NToO) program
to look for gamma-ray counterparts to neutrino alerts in follow-up
observations. The large field-of-view (FoV) and the rapidly re-position
capabilities of CTA’s telescopes working alongside the real time alert program
from IceCube (Blaufuss et al., 2019) look very promising for the upcoming
years. A more detailed discussion of the current work and development of the
NToO for CTA can be found in Satalecka et al. (2020).

From the kinematics of Equation 1, each $\pi^{+}$ will produce 3 neutrinos and
one $e^{+}$, which will carry 1/4 of the $\pi^{+}$ energy each. In the other
channel, the $\pi^{0}$ will produce a pair of photons, so that the observed
gamma-ray photon energy and the neutrino energy satisfy $E_{\nu}\sim
E_{\gamma}/2$ (Mücke et al., 1999). Also from the kinematics of the decay
chain, we have that $F_{\nu}\sim\frac{3}{4}F_{\pi^{+}}=\frac{3}{8}F_{\gamma}$
(Sahu et al., 2012, 2013). After neutrino oscillations, the expected flavor
ratio at earth will be $\nu_{e}:\nu_{\mu}:\nu_{\tau}=1:1:1$, and so the
estimated neutrino flux for only muon neutrinos will be a third of the all
flavour flux, then $F_{\nu_{\mu}}\sim\frac{1}{8}F_{\gamma}$.

For the VHE gamma-ray photons in consideration from these observations, the
energy threshold condition for the $\Delta-$resonance leads to a starting
point for the neutrino energy range (in the observer’s reference frame) of
around 0.04 TeV < $E_{\nu,\rm min}$ < 2.5 TeV. From the maximum proton energy
condition given by the Hillas criterion, the maximum neutrino energy in the
observer’s frame is expected to be $E_{\nu,\rm max}\sim 680$ TeV.

If we set an upper limit of $F_{\gamma}<10^{-4}$ MeV cm-2 s-1 (which is around
the maximum value reached during the flare) and if we assume a neutrino
spectrum of the type
$\frac{\rm{d}\it{N}}{\rm{d}\it{E}_{\nu}}=A_{\nu}E_{\nu}^{-2}$, then the
constant $A_{\nu}$ is estimated as a fraction of the gamma-ray flux. From the
highest flux state of the source, we have an approximation of $A_{\nu}\sim
3.1\times 10^{-12}$ TeV cm-2 s-1. Then, by integrating the neutrino spectrum
using the effective areas $A_{\rm eff}(E_{\nu})$ of the 59-strings
configuration (IC-59) operating in 2010 (Aartsen et al., 2013; Halzen &
Hooper, 2005), we can set an upper limit of $N_{\rm events}<0.14$ for the
expected number neutrino events that would have been detected during the Mrk
421 flaring event. While in this paper, we are focusing on the explanation of
the VHE gamma-ray emission rather than neutrino production, the possible
detection of neutrino events from flaring blazars is an interesting topic that
will be investigated in a future work on the photohadronic contributions.

On MJD 55267, 55274 and 55277, the photohadronic model would behave like a PL
that increases, if we extend the fit below $10^{5}$ MeV. This is related to
the type of seed photon input that we are using (a PL approximation) and does
not represent a realistic physical description of the SED. This is a caveat of
the current model which is focused on the contribution of the photohadronic
component at VHE and it is not designed to extend continuously to lower
energies. A more complex input model would improve the overall final fit but
at the expense of a larger number of free parameters. For the PL description
of the seed photon spectrum, 2 free parameters were optimised over bins
lasting a couple of days due to the lack of photon statistics; this was the
shortest time period over which it was possible to obtain Fermi spectra in the
> 100 MeV energy range. The parameters $\alpha$ and $A_{\gamma}$ coming from
the photohadronic contribution at VHE were calculated on daily intervals: this
discrepancy between the time bins was dictated by the differences in
instantaneous sensitivity between the space and ground-based telescopes and
might affect the accuracy of the final fit.

The uncertainty of the seed photon spectrum impacts the optimisation process
for the free model parameters $\alpha$ and $A_{\gamma}$ which are chosen to
get the best-fit values to the VHE data and therefore the minimum AIC value
for the model. Better data in the MeV energy range (1 MeV to 100 MeV) would
enable a better description of the input seed photons and hence an improvement
in the photohadronic fit and more reliable predictions. On this regard, the
All-sky Medium Energy Gamma-ray Observatory (AMEGO) mission (McEnery et al.,
2019), which is planned to operate from 200 keV to >10 GeV with $\sim 5\times$
better angular resolution than Fermi-LAT, would be very helpful. The
combination of MeV gamma-ray photon data with precise measurements from gamma
rays in the 10s of TeV regime would be ideal to further test hadronic emission
models.

##  8 Conclusions

The hadronic modelling of flaring episodes from blazars is a complex challenge
that has gained relevance in recent years. In this work, the scenario of a
dominant hadronic contribution in the VHE region of the SED coming from
$p\gamma$ interactions during the flaring period of Mrk 421 in 2010 was
studied. A photohadronic model with 4 free parameters was used and the gamma-
ray flux calculated using the $\Delta-$resonance approximation. For the target
photon spectrum, we used a PL description estimated from an analysis performed
of Fermi-LAT data. The injected proton spectrum assumed for the model was also
characterised by a PL. We were able to fit the VHE gamma-ray data with the
$p\gamma$ model on the days with sufficient photon statistics and according to
the AIC test, in all cases the $p\gamma$ model was favoured as a better fit
description than a one-zone leptonic model, and in comparison with the two-
zone SSC model, the $p\gamma$ was favoured by the AIC test on 3 out of the 5
days fitted (MJD 55267, 55276, 55277). The AIC test was inconclusive on MJD
55274 because the difference between the $p\gamma$ and the two-zone model was
not meaningful. On MJD 55266 the two-zone model was favoured as a better
description of the observations.

Our results therefore show that a dominant contribution from the photohadronic
component can be used to successfully fit the observations of a blazar flaring
episode, which shows the potential of including $p\gamma$ interactions in
blazar modelling. However, other contributions coming from leptonic processes,
synchrotron emission from the charged particles in the hadronic decay chain
and cascading effects can also play an important role at VHE. These will be
investigated in future works in the search to complement our model towards a
lepto-hadronic description.

To explore the neutrino/gamma-ray connection in the upcoming years, the next
generation of gamma-ray and neutrino observatories, such as CTA, AMEGO and
IceCube-Gen2, will play a crucial role; the improvements in observations at
VHE and follow-up programs will make possible to test hadronic components and
discriminate between pure leptonic and hadronic scenarios.

## Acknowledgements

We thank the anonymous referee for their comments and suggestions that helped
us to improve the quality and clarity of this article. ARdL acknowledges the
support of the National Council for Science and Technology from Mexico
(CONACYT). AMB and PMC acknowledge the financial support of the UK Science and
Technology Facilities Council consolidated grant ST/P000541/1.

## Data availability

This work has made use of public Fermi data obtained from the High Energy
Astrophysics Science Archive Research Center (HEASARC), provided by NASA
Goddard Space Flight Center. The day-by-day broadband SEDs data used from
Aleksić et al. (2015b) is publicly available at:
http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/578/A22. This research has made
use of Fermipy (Wood et al., 2017) to analyse the data from the LAT and Scipy
(https://www.scipy.org/) a Python-based ecosystem of open-source software for
mathematics, science, and engineering. Any extra data generated as part of
this project may be shared on a reasonable request to the corresponding
author.

## References

  * Aartsen et al. (2013) Aartsen M. G., et al., 2013, ApJ, 779, 132
  * Abdo et al. (2011) Abdo A. A., et al., 2011, ApJ, 736, 131
  * Abdollahi et al. (2020) Abdollahi S., et al., 2020, ApJS, 247, 33
  * Abeysekara et al. (2020) Abeysekara A. U., et al., 2020, ApJ, 890, 97
  * Aharonian et al. (2006) Aharonian F., et al., 2006, Nature, 440, 1018
  * Akaike (1974) Akaike H., 1974, IEEE Transactions on Automatic Control, 19, 716
  * Aleksić et al. (2012) Aleksić J., et al., 2012, Astroparticle Physics, 35, 435
  * Aleksić et al. (2015a) Aleksić J., et al., 2015a, A&A, 576, A126
  * Aleksić et al. (2015b) Aleksić J., et al., 2015b, A&A, 578, A22
  * Araudo et al. (2010) Araudo A. T., Bosch-Ramon V., Romero G. E., 2010, A&A, 522, A97
  * Atwood et al. (2009) Atwood W. B., et al., 2009, ApJ, 697, 1071
  * Baring (1997) Baring M. G., 1997, in Giraud-Heraud Y., Tran Thanh van J., eds, Very High Energy Phenomena in the Universe; Moriond Workshop. p. 97 (arXiv:astro-ph/9711177)
  * Blaufuss et al. (2019) Blaufuss E., Kintscher T., Lu L., Tung C. F., 2019, in 36th International Cosmic Ray Conference (ICRC2019). p. 1021 (arXiv:1908.04884)
  * Błażejowski et al. (2005) Błażejowski M., et al., 2005, ApJ, 630, 130
  * Bloom & Marscher (1996) Bloom S. D., Marscher A. P., 1996, ApJ, 461, 657
  * Böttcher et al. (2013) Böttcher M., Reimer A., Sweeney K., Prakash A., 2013, ApJ, 768, 54
  * Bozdogan (1987) Bozdogan H., 1987, Psychometrika, 52, 345
  * Brown et al. (2015) Brown A. M., Adams J., Chadwick P. M., 2015, MNRAS, 451, 323
  * Burnham & Anderson (2001) Burnham K. P., Anderson D. R., 2001, Wildlife Research, 28, 111
  * CTA Consortium et al. (2019) CTA Consortium et al., 2019, Science with the Cherenkov Telescope Array. WSP, doi:10.1142/10986
  * Cerruti et al. (2019) Cerruti M., Zech A., Boisson C., Emery G., Inoue S., Lenain J. P., 2019, MNRAS, 483, L12
  * Dar & Laor (1997) Dar A., Laor A., 1997, ApJ, 478, L5
  * Dermer & Schlickeiser (1993) Dermer C. D., Schlickeiser R., 1993, ApJ, 416, 458
  * Diltz & Böttcher (2016) Diltz C., Böttcher M., 2016, ApJ, 826, 54
  * Diltz et al. (2015) Diltz C., Böttcher M., Fossati G., 2015, ApJ, 802, 133
  * Dimitrakoudis et al. (2012) Dimitrakoudis S., Mastichiadis A., Protheroe R. J., Reimer A., 2012, A&A, 546, A120
  * Domínguez et al. (2011) Domínguez A., et al., 2011, MNRAS, 410, 2556
  * Finke (2018) Finke J. D., 2018, ApJ, 860, 178
  * Forster (2000) Forster M. R., 2000, Journal of Mathematical Psychology, 44, 205
  * Gaisser et al. (1995) Gaisser T. K., Halzen F., Stanev T., 1995, Physics Reports, 258, 173
  * Ghisellini et al. (2005) Ghisellini G., Tavecchio F., Chiaberge M., 2005, A&A, 432, 401
  * Giommi et al. (2020) Giommi P., Padovani P., Oikonomou F., Glauch T., Paiano S., Resconi E., 2020, arXiv e-prints, p. arXiv:2003.06405
  * Glüsenkamp (2016) Glüsenkamp T., 2016, in European Physical Journal Web of Conferences. p. 05006 (arXiv:1502.03104), doi:10.1051/epjconf/201612105006
  * Halzen & Hooper (2005) Halzen F., Hooper D., 2005, Astroparticle Physics, 23, 537
  * Halzen et al. (2019) Halzen F., Kheirandish A., Weisgarber T., Wakely S. P., 2019, ApJ, 874, L9
  * Harris et al. (2014) Harris J., Chadwick P. M., Daniel M. K., 2014, MNRAS, 441, 3591
  * Hillas (1984) Hillas A. M., 1984, ARA&A, 22, 425
  * Homan et al. (2015) Homan D. C., Lister M. L., Kovalev Y. Y., Pushkarev A. B., Savolainen T., Kellermann K. I., Richards J. L., Ros E., 2015, ApJ, 798, 134
  * Hümmer et al. (2010) Hümmer S., Rüger M., Spanier F., Winter W., 2010, ApJ, 721, 630
  * IceCube Collaboration et al. (2018a) IceCube Collaboration et al., 2018a, Science, 361, 147
  * IceCube Collaboration et al. (2018b) IceCube Collaboration et al., 2018b, Science, 361, eaat1378
  * Kadler et al. (2016) Kadler M., et al., 2016, Nature Physics, 12, 807
  * Kildea et al. (2007) Kildea J., et al., 2007, Astroparticle Physics, 28, 182
  * Krauß et al. (2014) Krauß F., et al., 2014, A&A, 566, L7
  * Krauß et al. (2015) Krauß F., et al., 2015, in 5th International Fermi Symposium. (arXiv:1502.02147)
  * Lewis et al. (2011) Lewis F., Butler A., Gilbert L., 2011, Methods in Ecology and Evolution, 2, 155
  * Liu et al. (2019) Liu R.-Y., Wang K., Xue R., Taylor A. M., Wang X.-Y., Li Z., Yan H., 2019, Phys. Rev. D, 99, 063008
  * MacDonald et al. (2015) MacDonald N. R., Marscher A. P., Jorstad S. G., Joshi M., 2015, ApJ, 804, 111
  * Mannheim (1993) Mannheim K., 1993, A&A, 269, 67
  * Mannheim & Biermann (1992) Mannheim K., Biermann P. L., 1992, A&A, 253, L21
  * Maraschi et al. (1992) Maraschi L., Ghisellini G., Celotti A., 1992, ApJ, 397, L5
  * Marscher et al. (2008) Marscher A. P., et al., 2008, Nature, 452, 966
  * Marscher et al. (2010) Marscher A. P., et al., 2010, ApJ, 710, L126
  * Mattox et al. (1996) Mattox J. R., et al., 1996, ApJ, 461, 396
  * McEnery et al. (2019) McEnery J., et al., 2019, in BAAS. p. 245 (arXiv:1907.07558)
  * Meli et al. (2008) Meli A., Becker J. K., Quenby J. J., 2008, A&A, 492, 323
  * Mücke & Protheroe (2001) Mücke A., Protheroe R. J., 2001, Astroparticle Physics, 15, 121
  * Mücke et al. (1999) Mücke A., Rachen J. P., Engel R., Protheroe R. J., Stanev T., 1999, Publ. Astron. Soc. Australia, 16, 160
  * Mücke et al. (2003) Mücke A., Protheroe R. J., Engel R., Rachen J. P., Stanev T., 2003, Astroparticle Physics, 18, 593
  * Murase et al. (2018) Murase K., Oikonomou F., Petropoulou M., 2018, ApJ, 865, 124
  * Oikonomou et al. (2019) Oikonomou F., Murase K., Padovani P., Resconi E., Mészáros P., 2019, MNRAS, 489, 4347
  * Padovani et al. (2016) Padovani P., Resconi E., Giommi P., Arsioli B., Chang Y. L., 2016, MNRAS, 457, 3582
  * Padovani et al. (2019) Padovani P., Oikonomou F., Petropoulou M., Giommi P., Resconi E., 2019, MNRAS, 484, L104
  * Palladino et al. (2019) Palladino A., Rodrigues X., Gao S., Winter W., 2019, ApJ, 871, 41
  * Perkins et al. (2009) Perkins J. S., Maier G., The VERITAS Collaboration 2009, arXiv e-prints, p. arXiv:0912.3841
  * Petropoulou et al. (2020) Petropoulou M., et al., 2020, ApJ, 891, 115
  * Pichel (2009) Pichel A., 2009, arXiv e-prints, p. arXiv:0908.0010
  * Ptitsyna & Troitsky (2010) Ptitsyna K. V., Troitsky S. V., 2010, Physics Uspekhi, 53, 691
  * Punch et al. (1992) Punch M., et al., 1992, Nature, 358, 477
  * Rachen & Mészáros (1998) Rachen J. P., Mészáros P., 1998, Phys. Rev. D, 58, 123005
  * Righi et al. (2019) Righi C., Tavecchio F., Inoue S., 2019, MNRAS, 483, L127
  * Rodrigues et al. (2019) Rodrigues X., Gao S., Fedynitch A., Palladino A., Winter W., 2019, ApJ, 874, L29
  * Sahayanathan & Godambe (2012) Sahayanathan S., Godambe S., 2012, MNRAS, 419, 1660
  * Sahu et al. (2012) Sahu S., Zhang B., Fraija N., 2012, Phys. Rev. D, 85, 043012
  * Sahu et al. (2013) Sahu S., Oliveros A. F. O., Sanabria J. C., 2013, Phys. Rev. D, 87, 103015
  * Sahu et al. (2016) Sahu S., Miranda L. S., Rajpoot S., 2016, Eur. Phys. J. C, 76, 127
  * Sahu et al. (2017) Sahu S., de León A. R., Miranda L. S., 2017, Eur. Phys. J. C, 77, 741
  * Sahu et al. (2018a) Sahu S., de León A. R., Nagataki S., 2018a, Eur. Phys. J. C, 78, 484
  * Sahu et al. (2018b) Sahu S., de León A. R., Nagataki S., Gupta V., 2018b, Eur. Phys. J. C, 78, 557
  * Satalecka et al. (2020) Satalecka K., et al., 2020, PoS, ICRC2019, 784
  * Sikora et al. (1994) Sikora M., Begelman M. C., Rees M. J., 1994, ApJ, 421, 153
  * Tanaka et al. (2017) Tanaka Y. T., Buson S., Kocevski D., 2017, The Astronomer’s Telegram, 10791, 1
  * Tavecchio et al. (1998) Tavecchio F., Maraschi L., Ghisellini G., 1998, ApJ, 509, 608
  * Tavecchio et al. (2009) Tavecchio F., Ghisellini G., Ghirlanda G., Costamante L., Franceschini A., 2009, MNRAS, 399, L59
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Walker et al. (2016) Walker R. C., Hardee P. E., Davies F., Ly C., Junor W., Mertens F., Lobanov A., 2016, Galaxies, 4, 46
  * Winter & Gao (2019) Winter W., Gao S., 2019, in 36th International Cosmic Ray Conference (ICRC2019). p. 1032 (arXiv:1909.06289)
  * Wood et al. (2017) Wood M., Caputo R., Charles E., Di Mauro M., Magill J., Perkins J. S., Fermi-LAT Collaboration 2017, in 35th International Cosmic Ray Conference (ICRC2017). p. 824 (arXiv:1707.09551)