# Sub-GeV Gamma Rays from Nearby Seyfert Galaxies and Implications for Coronal
Neutrino Emission

Kohta Murase11affiliation: Department of Physics; Department of Astronomy &
Astrophysics; Center for Multimessenger Astrophysics, Institute for
Gravitation and the Cosmos, The Pennsylvania State University, University
Park, PA 16802, USA 22affiliation: School of Natural Sciences, Institute for
Advanced Study, Princeton, NJ 08540, USA 33affiliation: Center for
Gravitational Physics and Quantum Information, Yukawa Institute for
Theoretical Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan
Christopher M. Karwin44affiliation: NASA Postdoctoral Program Fellow, NASA
Goddard Space Flight Center, Greenbelt, MD, 20771, USA Shigeo S.
Kimura55affiliation: Frontier Research Institute for Interdisciplinary
Sciences; Astronomical Institute, Graduate School of Science, Tohoku
University, Sendai 980-8578, Japan Marco Ajello66affiliation: Department of
Physics and Astronomy, Clemson University, Clemson, SC 29634, USA Sara
Buson77affiliation: Fakultät für Physik und Astronomie, Julius-Maximilians-
Universität Würzburg, Würzburg D-97074, Germany

###### Abstract

Recent observations of high-energy neutrinos by IceCube and $\gamma$ rays by
the Fermi Large Area Telescope (LAT) and the MAGIC telescope have suggested
that neutrinos are produced in $\gamma$-ray opaque environments in the
vicinity of supermassive black holes. In this work, we present 20 MeV – 1 TeV
spectra of three Seyfert galaxies whose nuclei are predicted to be active in
neutrinos, NGC 4151, NGC 4945 and the Circinus Galaxy, using 14.4 years of the
Fermi LAT data. In particular, we find evidence of sub-GeV excess emission
that can be attributed to $\gamma$ rays from NGC 4945, as was also seen in NGC
1068. These spectral features are consistent with predictions of the
magnetically-powered corona model, and we argue that NGC 4945 is among the
brightest neutrino active galaxies detectable for KM3Net and Baikal-GVD. On
the other hand, in contrast to other reported results, we do not detect
$\gamma$ rays from NGC 4151, which constrains neutrino emission from the
accretion shock model. Future neutrino detectors such as IceCube-Gen2 and MeV
$\gamma$-ray telescopes such as AMEGO-X will be crucial for discriminating
among the theoretical models.

##  1 Introduction

Since the IceCube Collaboration announced the first evidence for high-energy
cosmic neutrinos in 2013 (Aartsen et al., 2013a, b), the detection of many
astrophysical TeV–PeV neutrinos have been reported. They are produced by
hadronuclear ($pp$) and photohadronic ($p\gamma$) interactions, inevitably
producing a similar amount of $\gamma$ rays. These $\gamma$ rays can be
absorbed by the Breit-Wheeler process ($\gamma+\gamma\rightarrow e^{+}+e^{-}$)
and reprocessed to lower-energy $\gamma$ rays, which can contribute to the
cosmic $\gamma$-ray background measured by Fermi LAT (Ackermann et al., 2015).
If all the neutrino sources were $\gamma$-ray transparent, $\gamma$ rays
accompanied by the large all-sky neutrino flux at neutrino energies of
$E_{\nu}\sim 1-10$ TeV (Aartsen et al., 2015, 2020a) would overshoot the
observed diffuse isotropic $\gamma$-ray background. This necessitates the
existence of hidden cosmic-ray (CR) accelerators that are opaque to GeV-TeV
$\gamma$ rays, and the vicinity of supermassive black holes (SMBHs) is among
the most promising sites (Murase et al., 2016).

The existence of hidden neutrino sources has further been supported by the
recent IceCube observation of neutrino signals from the nearby Seyfert galaxy,
NGC 1068 (Aartsen et al., 2020b; Abbasi et al., 2022). The GeV-TeV
$\gamma$-ray fluxes seen by LAT (Abdollahi et al., 2020; Ajello et al., 2023)
and MAGIC upper limits (Acciari et al., 2019) are much lower than the TeV
neutrino flux measured with IceCube. The region has to be compact for the
system to be opaque to GeV $\gamma$ rays (Murase, 2022), as expected in disk-
corona models of active galactic nuclei (AGNs) (Murase et al., 2020;
Kheirandish et al., 2021). This is the main consequence of the GeV-TeV
$\gamma$-ray and TeV-PeV neutrino observations. It has been predicted that
X-ray bright AGNs are promising targets for neutrino observations (see Murase
& Stecker, 2022, as a review), while TeV $\gamma$ rays are cascaded down to
sub-GeV energies, which will be ideal targets for MeV–GeV $\gamma$-ray
telescopes.

In this work, we investigate multimessenger implications of three nearby
Seyfert galaxies, NGC 4151, NGC 4945, and the Circinus Galaxy, using Fermi LAT
data. We analyze 14.4 years of Fermi LAT data to obtain the $\gamma$-ray
fluxes from these galaxies. Then, we compare the sub-GeV $\gamma$-ray data
with CR-induced electromagnetic cascade emission and discuss corresponding
neutrino fluxes. We use $Q_{x}=Q/10^{x}$ in cgs units and assume cosmological
parameters with $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and $h=0.7$.

##  2 Neutrino and $\gamma$-ray Emission from the Hearts of AGNs

Multimessenger data from NGC 1068 suggest that high-energy neutrinos from
Seyfert galaxies are likely to originate from compact regions with a size of
$R\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}30-100~{}R_{S}$ (Murase, 2022), and X-ray
emitting coronae or possible shock regions, located at the vicinity of SMBHs,
are the most plausible site of the neutrino production (Murase et al., 2020;
Inoue et al., 2020; Eichmann et al., 2022; Inoue et al., 2022). In this work,
we consider neutrinos and $\gamma$ rays from such coronal regions with
$R=10-30~{}R_{S}$.

The spectral energy distribution (SED) of AGNs is modeled empirically (Murase
et al., 2020), based on the parameterization as a function of the SMBH mass
$M_{\rm BH}$ and the intrinsic X-ray luminosity $L_{X}$. We use $M_{\rm
BH}=1.0\times 10^{7}~{}M_{\odot}$ and $L_{X}=2.6\times 10^{42}~{}{\rm
erg}~{}{\rm s}^{-1}$ for NGC 4151 (Bentz et al., 2022; Koss et al., 2022),
$M_{\rm BH}=1.4\times 10^{6}~{}M_{\odot}$ and $L_{X}=3.3\times 10^{42}~{}{\rm
erg}~{}{\rm s}^{-1}$ (Greenhill et al., 1997; Puccetti et al., 2014) for NGC
4945, and $M_{\rm BH}=1.7\times 10^{6}~{}M_{\odot}$ and
$L_{X}=10^{42.5}~{}{\rm erg}~{}{\rm s}^{-1}$ for the Circinus Galaxy
(Greenhill et al., 2003; Tanimoto et al., 2019), respectively. We checked that
the adopted SED model (see Figure 2 of Murase et al., 2020) agrees with the
observed SEDs of the AGNs considered in this work.

High-energy neutrinos are produced through meson production by $p\gamma$ and
$pp$ interactions. We calculate neutrino and $\gamma$-ray emission,
considering CR-induced cascades, with the AGN module of the Astrophysical
Multimessenger Emission Simulator (AMES) (e.g., Murase, 2018; Zhang & Murase,
2023). Details of our steady-state corona model are described in Murase et al.
(2020), but we use the detailed cross sections of the photomeson production
and Bethe-Heitler pair production processes as in Murase (2022). The model
predicts that TeV neutrinos must be accompanied by sub-GeV $\gamma$ rays,
predicting $\gamma$-ray and neutrino energy fluxes of
$E_{\gamma}F_{E_{\gamma}}\sim(0.1-1)E_{\nu}F_{E_{\nu}}$ within 1 order of
magnitude. This is a consequence of the multimessenger connection, which is
insensitive to model parameters.

###  2.1 Magnetically-powered coronae

Magnetohydrodynamic (MHD) simulations for accretion flows revealed that the
magnetorotational instability (MRI) and subsequent dynamos amplify the
magnetic field and develop strong turbulence in accretion flows (e.g., Balbus
& Hawley, 1998; Kimura et al., 2019). It is believed that magnetic dissipation
in the accretion flows heat the plasma and form a hot region called a corona.
The observed X-ray emission with a power-law spectrum originates from the
Comptonization of disk photons by thermal or bulk electrons.

The MHD turbulence scatters CRs, and the scattered particles randomly change
their energies, and effectively gain energy from the turbulence. In the
coronal region of AGNs, Alfvén and turbulent velocities are mildly
relativistic, enabling the turbulence to accelerate CRs up to TeV-PeV
energies. In this process, the resulting CR spectrum can be hard with a
spectral index of $s\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}2$ (for $dN_{\rm
CR}/d\varepsilon_{p}\propto\varepsilon_{p}^{-s}$, where $\varepsilon_{p}$ is
the CR proton energy in the source frame), and the maximum CR energy typically
lies in the range of $\sim 0.1-1$ PeV with a reasonable acceleration
efficiency ($\eta_{\rm tur}\sim 10-100$). The stochastic acceleration by
large-scale turbulence has been studied by test-particle simulations with MHD
(Kimura et al., 2016, 2019; Sun & Bai, 2021) and large-scale particle-in-cell
(PIC) simulations (Zhdankin et al., 2019; Comisso & Sironi, 2019, 2022).
Alternatively, magnetic reconnections accelerate non-thermal particles
efficiently, especially in the relativistic regime ($\sigma_{\rm
mag}=B^{2}/(8\pi mc^{2})>1$) (Hoshino & Lyubarsky, 2012; Guo et al., 2020).
PIC simulations have demonstrated that the development of MRI turbulence
triggers magnetic reconnections, leading to non-thermal particle production
(e.g., Hoshino, 2015; Kunz et al., 2016; Bacchini et al., 2022). Kheirandish
et al. (2021) suggested such a magnetic reconnection model to explain the
neutrino data of NGC 1068 (see also Mbarek et al., 2023; Fiorillo et al.,
2023).

In this work, we consider the magnetically-powered corona model with
stochastic acceleration (Murase et al., 2020), where we solve the Fokker-
Planck equation to obtain the CR spectrum and calculate the resulting neutrino
and $\gamma$-ray spectra taking into account electromagnetic cascades through
$\gamma\gamma\rightarrow e^{+}e^{-}$, synchrotron and inverse-Compton emission
processes. The coronal plasma is assumed to be geometrically thick and
advected to a SMBH with the infall velocity $V_{\rm
infall}=\alpha\sqrt{GM_{\rm BH}/R}$, where $\alpha=0.1$ is the viscosity
parameter. The magnetic field is obtained through the plasma beta set to
$\beta=1$, where synchrotron cascade emission is important for magnetized
coronae with $\beta\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}1$ and the cascaded $\sim 0.1–10$ MeV
$\gamma$-ray flux is a robust signature (Murase et al., 2020). Then, the
remaining principal parameters of the model are for CR inputs regarding the CR
normalization and maximum energy. We consider two specific cases. For Model A
we assume the CR injection efficiency111The corresponding CR pressure is 50%
of the thermal pressure with the virial temperature for
$L_{X}=3\times{10}^{43}~{}{\rm erg}~{}{\rm s}^{-1}$ and $d=10$ Mpc (Ajello et
al., 2023). and $\eta_{\rm tur}=70$ that are used to explain the IceCube data
of NGC 1068 (see Figure 4 of Abbasi et al., 2022) with $R=30~{}R_{S}$. For
Model B, the CR pressure $P_{\rm CR}$ is set to 10% of the thermal pressure
with the virial temperature ($P_{\rm vir}$). Because the coronal pressure
would be dominated by magnetic fields, this assumption is reasonable. The
emission radius is assumed to be $R=10~{}R_{S}$, in which more neutrinos and
$\gamma$ rays may be produced due to $p\gamma$ interactions. For the three
AGNs studied in this work, Model B gives more optimistic predictions for the
neutrino and $\gamma$-ray fluxes.

###  2.2 Accretion shocks

Accretion shocks have been discussed as a CR acceleration site for a long time
(e.g., Stecker et al., 1991; Inoue et al., 2020; Anchordoqui et al., 2021). In
this model, a spherical accretion flow, aside from the accretion disk flow, is
assumed to form a shock with a velocity of $V_{s}\approx V_{\rm
ff}=\sqrt{GM_{\rm BH}/R}$ in the vicinity of the SMBH, and
$s\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}2$ is expected based on diffusive shock
acceleration theory (e.g., Drury, 1983). However, the neutrino data of NGC
1068 suggest that the particle acceleration efficiency, $\eta_{\rm sho}$, has
to be much larger than the canonical value in the Bohm limit ($\eta_{\rm
sho}\sim 1-10$) not to violate the constraint on the neutrino break/cutoff
(see Figure 5 of Kheirandish et al., 2021). Thus, in this work, we ad hoc fix
the proton maximum energy to 100 TeV (implying $\eta_{\rm sho}\gg 10$) and
consider the largest CR luminosity $L_{\rm CR}$ allowed by the LAT data. Then,
we calculate neutrino and $\gamma$-ray spectra assuming the steady-state
injection of CRs for given $R$ and SEDs.

In the accretion shock model, magnetic fields are expected to be weaker than
the magnetically-powered corona model, and we set $B=30$ G throughout this
work (Inoue & Khangulyan, 2023). For such low magnetizations, the results on
$\gamma$-ray spectra are largely insensitive to CR spectra as long as
electromagnetic cascades proceed mainly through $\gamma\gamma\rightarrow
e^{+}e^{-}$ and inverse-Compton emission. Note that we only consider CR-
induced cascade $\gamma$ rays without any primary non-thermal electron
component, unlike Inoue & Khangulyan (2023) that focused on $\gamma$ rays from
primary electrons only with attenuation (without cascades). The primary
electron component would add an additional sub-GeV gamma-ray component, which
would result in a lower neutrino flux in order to be consistent with the given
sub-GeV gamma-ray data.

##  3 Data Analysis

Table 1: Summed Likelihood Components

Energy [GeV] | Z Max [deg] | PSF Types
---|---|---
$0.02-0.1$ | 80 | 3
$0.1-0.3$ | 90 | 2, 3
$0.3-1$ | 100 | 1, 2, 3
$1-1000$ | 105 | 0, 1, 2, 3

Note. — These energy bins distinguish different components of the joint
likelihood analysis, however, the underlying fits still use 8 energy bins per
decade, as decribed in the text.

Figure 1: Left: TS map of the NGC 4151 field. Contours show the 95%
localization uncertainty. Red contours are from the 4FGL-DR3 catalog and the
pink contour is calculated based on the peak in the TS map. Blue circles show
the locations of the BL Lacs which are associated with the 4FGL-DR3 sources
(4FGL 1211.6+3901 towards the lower left ($0.43^{\circ}$ separation) and 4FGL
1210.3+3928 closest to NGC 4151 ($0.08^{\circ}$ separation)). Right: Same as
the left plot but zoomed in on NGC 4151. Overlaid are multiwavelength images
of NGC 4151 in both optical (black and white image) and radio (H I, shown with
the rainbow colormap).

We study three nearby Seyfert galaxies, NGC 4151, NGC 4945, and the Circinus
Galaxy. These objects were all reported as LAT sources and are also very
bright in the hard X-ray band. In the magnetically-powered corona model, the
neutrino luminosity is roughly proportional to the X-ray luminosity (Murase et
al., 2020; Kheirandish et al., 2021), and thus NGC 1068 and NGC 4151 (NGC 4945
and Circinus) are predicted to be among the brightest sources in the northern
(southern) neutrino sky among Seyfert galaxies in the Swift BASS catalog with
LAT detection222With the Swift BASS catalog considering the 0.2-10 keV and
14-195 keV bands for the intrinsic X-ray flux, the brightest Seyfert galaxies
in the northern sky including the near-horizon are NGC 1068 and NGC 4151, and
NGC 1068 is expected to be the most promising neutrino source especially with
the NuSTAR flux. In the southern sky, NGC 4945, ESO 138-1 and Circinus are the
brightest. We focus on NGC 4945 and Circinus, which were reported as LAT
sources, and the high luminosity of ESO 138-1 would suppress coronal
$\gamma$-ray emission more severely.. Prior to this work, the results of NGC
1068 (being the brightest in the IceCube sky) are presented in Ajello et al.
(2023).

We analyze data collected by the LAT (Atwood et al., 2009, 2013; Bruel et al.,
2018) from 2008 August 04 to 2023 January 5 (14.4 years). We select an energy
range of 20 MeV to 1 TeV, and bin the data using 8 energy bins per decade. The
pixel size is $0.08^{\circ}$. We use a $10^{\circ}\times 10^{\circ}$ region of
interest (ROI) centered on each respective galaxy. The standard data filters
are used: DATA_QUAL$>$0 and LAT_CONFIG==1. The analysis is performed using
Fermipy (v1.2)333Available at https://fermipy.readthedocs.io/en/latest/, which
utilizes the underlying Fermitools (v2.2.0).

We select photons corresponding to the P8R3_SOURCE_V3 instrument response. In
order to optimize the sensitivity of our analysis we implement a joint
likelihood fit with the four point spread function (PSF) event types available
in the Pass 8 data set444For more information on the different PSF types see
https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data/LAT_DP.html..
The data are divided into quartiles corresponding to the quality of the
reconstructed direction, from the lowest quality quartile (PSF0) to the best
quality quartile (PSF3). Each sub-selection has its own binned likelihood
instance that is combined in a summed likelihood function for the ROI. This is
easily implemented in Fermipy by specifying the components section in the
configuration file. We include different event types depending on the energy
interval. This is motivated by the energy-dependence of the LAT instrument
response. Additionally, in order to reduce the contamination from the Earth’s
limb, we apply an energy-dependent cut on the maximum zenith angle. These
selections are similar to those used for the 4FGL-DR3 catalog (Abdollahi et
al., 2022), and they are summarized in Table 1. Note, however, that the 4FGL-
DR3 did not use photon energies below 50 MeV. Each PSF type also has its own
corresponding isotropic spectrum, namely, iso_P8R3_SOURCE_V3_PSFi_v1, for $i$
ranging from 0 to 3. The isotropic spectrum file is defined from 34 MeV to 878
GeV, and so we use a power-law extrapolation for energies outside this range.
The Galactic diffuse emission is modeled using the standard component
(gll_iem_v07). It is only defined between 50 MeV to 814 GeV, and so likewise
we use a power law extrapolation. The point source emission is modeled using
the 4FGL-DR3 catalog (gll_psc_v28) (Abdollahi et al., 2022). In order to
account for photon leakage from sources outside of the ROI due to the PSF of
the detector, the model includes all 4FGL sources within a $15^{\circ}\times
15^{\circ}$ region. The energy dispersion correction (edisp_bins=–1) is
enabled for all sources except the isotropic component.

Before calculating SEDs, we first perform an initial fit for each of the three
ROIs, in which we optimize the model, freeing the normalization and index of
the Galactic diffuse, the normalization of the isotropic, and all point
sources (normalization and index) within 3∘. We then find new sources using
the Fermipy function find_sources, which generates TS maps and identifies new
sources based on peaks in the TS. The TS maps are generated using a power-law
spectral model with a photon index of $2.0$. The minimum separation between
two point sources is set to $0.5^{\circ}$, and the minimum TS for including a
source in the model is set to 16. Finally, we perform a second fit in a
similar way as the first, with the exception that we free all point sources
within 5∘, having a TS$>25$.

Among the three Seyfert galaxies, the LAT observations of NGC 4151 are
complicated by the presence of two nearby BL Lacs, 4FGL J1211.6+3901 and 4FGL
1210.3+3928. In order to determine the most likely origin of the corresponding
$\gamma$-ray emission, we re-localize the point source closest to NGC 4151.
This is done following a similar procedure as the ROI optimization described
above. We use the same data selection, with the exception that we use a pixel
size of 0.04∘. In the first step of the optimization, we only include point
sources from the 4FGL-DR3, i.e. we do not include a point source at the
location of NGC 4151. Then, before finding new sources, we remove 4FGL
1210.3+3928 from the model, thus allowing the fitting algorithm to find the
best-fit location of the associated source. The resulting localization is
shown in Figure 1. The pink contour shows the 95% localization uncertainty,
based on the peak in the TS map. Red contours show the 95% localization
uncertainty from the 4FGL-DR3, and the blue circles show the BL Lacs which are
associated with the 4FGL sources. As can be clearly seen, the $\gamma$-ray
source near NGC 4151 is most likely associated with 4FGL 1210.3+3928 (and
hence, the BL Lac), consistent with the 4FGL catalog. We also verified that
the results remain consistent when using energies $E_{\gamma}>1$ GeV.
Additionally, we checked with the preliminary version of the 4FGL-DR4, also
finding consistent results. We therefore keep 4FGL 1210.3+3928 in the model
when calculating the SED for NGC 4151.

After optimizing the ROIs, SEDs are calculated for each source using the
Fermipy SED analysis. This method computes an SED by performing independent
fits for the flux normalization in each energy bin. For the calculation, we
combine the original energy bins into larger bins, using 2 bins per decade
between 100 MeV to 100 GeV, as well as a single bin between $20-100$ MeV, and
a single bin between $100-1000$ GeV. The normalization in each bin is fit
using a power-law spectral model with a fixed index of 2.0. The parameters of
the background components are held fixed at their best-fit values from the
baseline fit. Further details of the SED analyses for the three Seyfert
galaxies are provided in Appendix A. The resulting SEDs of NGC 4151, NGC 4945
and Circinus are shown in Figures 2 and 3.

Figure 2:  Multimessenger SEDs of NGC 4151, NGC 4945 and Circinus galaxy in
the magnetically-powered corona model. Black points and upper limits are
Fermi-LAT data in $\gamma$ rays (this work). The error bars on the LAT data
points are 1$\sigma$, and the ULs are plotted for bins with TS$<$4, and they
are shown at the $95\%$ C.L. The sky-blue shaded region represents IceCube
data in neutrinos (Glauch et al., 2023; Liu et al., 2023). Sensitivities for
e-ASTROGAM (De Angelis et al., 2017) and GRAMS (Aramaki et al., 2020) with an
effective exposure time of 1 yr are overlaid, together with the AMEGO-X
sensitivity for the 3 yr mission (Caputo et al., 2022). For Model A, we use
parameters to explain NGC 1068 as in Murase et al. (2020), where the emission
radius is set to $R=30R_{S}$ with $\eta_{\rm tur}=70$. For Model B, we assume
$R=10R_{S}$ and the CR pressure is set to $P_{\rm CR}/P_{\rm vir}=10$%, where
$\eta_{\rm tur}=100$ for NGC 4151 and $\eta_{\rm tur}=10$ for NGC 4945 and
Circinus are used.

Figure 3:  Same as Figure 2 but in the accretion shock model. The maximum
injection luminosities of CRs for $R=30R_{S}$ are $L_{\rm CR}=3\times
10^{42}~{}{\rm erg}~{}{\rm s}^{-1}$ (NGC 4151), $L_{\rm CR}=6\times
10^{42}~{}{\rm erg}~{}{\rm s}^{-1}$ (NGC 4945), and $L_{\rm CR}=10^{43}~{}{\rm
erg}~{}{\rm s}^{-1}$ (Circinus). For $R=10R_{S}$, they are $L_{\rm
CR}=10^{43}~{}{\rm erg}~{}{\rm s}^{-1}$, $L_{\rm CR}=2\times 10^{43}~{}{\rm
erg}~{}{\rm s}^{-1}$, and $L_{\rm CR}=3\times 10^{43}~{}{\rm erg}~{}{\rm
s}^{-1}$, respectively.

##  4 Implications and Discussions

GeV $\gamma$ rays mainly interact with X rays from the corona. The two-photon
annihilation optical depth for a photon index of $\Gamma_{X}\approx 2$ is
(Murase, 2022)

$\displaystyle\tau_{\gamma\gamma}(\varepsilon_{\gamma})$
$\displaystyle\approx$
$\displaystyle\eta_{\gamma\gamma}\sigma_{T}Rn_{X}{(\varepsilon_{\gamma}/\tilde{\varepsilon}_{\gamma\gamma-X})}^{\Gamma_{X}-1}$
(1) $\displaystyle\simeq$ $\displaystyle 24~{}{[R/(30R_{S})]}^{-1}M_{\rm
BH,7}^{-1}\tilde{L}_{X,42.4}(\varepsilon_{\gamma}/1~{}\rm
GeV),\,\,\,\,\,\,\,\,\,$

where $\eta_{\gamma\gamma}\sim 0.1$ is a numerical coefficient depending on
$\Gamma_{X}$, $\sigma_{T}\approx 6.65\times{10}^{-25}~{}{\rm cm}^{2}$,
$\tilde{\varepsilon}_{\gamma\gamma-X}=m_{e}^{2}c^{4}/\varepsilon_{X}\simeq
0.26~{}{\rm GeV}~{}{(\varepsilon_{X}/1~{}{\rm keV})}^{-1}$, $\varepsilon_{X}$
is the reference X-ray energy, and $n_{X}\approx\tilde{L}_{X}/(2\pi
R^{2}c\varepsilon_{X})$ is used. Here $\tilde{L}_{X}$ is the differential
X-ray luminosity. The numerical results shown in Figures 2 and 3 are
consistent with equation (1), and we may expect $\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}0.1$ GeV $\gamma$
rays. For NGC 4945 and Circinus, while GeV emission is detected, its origin
should be different, e.g., star-forming activities, as seen in NGC 1068
(Ajello et al., 2023). Note that hadronic emission from the starburst region
predicts a low-energy break in the $\gamma$-ray spectrum due to the decay of
neutral pions, although there could be leptonic contributions from inverse-
Compton and bremsstrahlung emission. However, detailed studies of the GeV
emission are beyond the scope of this work.

NGC 4151 is known to be a Seyfert 1.5 galaxy, showing the characteristic
features of both types 1 (dominated by broad-line components) and 2 (dominated
by narrow-line components). NGC 4151 is Compton thin unlike the other two
Seyferts and NGC 1068. The escape of GeV $\gamma$ rays is easier than in NGC
1068 due to its lower intrinsic X-ray luminosity, which makes the strong
limits in the GeV range important. As shown in Figures 2 and 3 left, we find
that the magnetically-powered corona and accretion shock models are consistent
with the upper limits obtained in this work. Nevertheless, the data give
interesting constraints on neutrino emission. In the magnetically-powered
corona model, Model B predicts a neutrino flux that can explain the possible
neutrino excess emission found in the IceCube data (Glauch et al., 2023; Liu
et al., 2023; Neronov et al., 2023). In the case of NGC 4151, the effective
maximum CR energy555In the diffusive-escape-limited case, the maximum energy
found in the CR distribution is effectively $\sim 10-30$ higher than the
energy inferred from equating the acceleration time with diffusive escape time
(Becker et al., 2006; Kimura et al., 2015). is limited by the diffusive escape
rather than the Bethe-Heitler cooling process that is relevant for luminous
AGNs including NGC 1068. On the other hand, due to the LAT upper limits in the
GeV range, the neutrino flux is unlikely to be explained in the accretion
shock model even if magnetic fields change. This demonstrates that
observations of NGC 4151-like galaxies are useful for discriminating between
the magnetically-powered corona and accretion shock models. When we include
not only 2-10 keV but also 14-195 keV data, the neutrino brightness of NGC
4151 is next to NGC 1068, where the ranking is higher than that in Kheirandish
et al. (2021), and NGC 4151 will be an important target for IceCube-Gen2 which
can reach a sensitivity of $E_{\nu}F_{E_{\nu}}\sim{10}^{-9}~{}{\rm GeV}~{}{\rm
cm}^{-2}~{}{\rm s}^{-1}$ (Aartsen et al., 2021), especially in the
magnetically-powered corona model (see Figure 2 left).

For NGC 4945, as shown in the middle panels of Figures 2 and 3, we find that
GeV and higher-energy emission has a break at $E_{\gamma}\sim 1$ GeV, which is
consistent with a pionic $\gamma$-ray component from
$\pi^{0}\rightarrow\gamma\gamma$, as expected in models relying on wind-torus
interactions (Inoue et al., 2022) and/or star-forming activities (Aguilar-Ruiz
et al., 2021). In addition, the LAT data around $\sim 30$ MeV shows a clear
excess compared to the pionic $\gamma$-ray component. However, it is important
to note that the LAT performance below 50 MeV quickly deteriorates, and the
instrument response at these low energies has not been well characterized. We
have performed a number of tests to validate the robustness of the observed
low-energy spectral feature in NGC 4945 (summarized in Appendix A), and none
of them have shown any indication of the signal being spurious. But
nevertheless, this feature should still be interpreted with caution due to the
limitations of the LAT performance below 50 MeV. That said, this signal is
very interesting because a similar excess flux over the pionic $\gamma$-ray
component was reported for NGC 1068 (Ajello et al., 2023). More intriguingly,
these sub-GeV $\gamma$-ray signatures of the two most promising neutrino-
bright AGNs are consistent with predictions of the magnetically-powered corona
model, and the corresponding LAT data can be explained if the CR pressure is
10% of the thermal pressure with the virial temperature (see Figure 2 middle).
For Circinus, as shown in Figures 2 and 3 right, we only have upper limits at
sub-GeV energies, which are also consistent with both of the corona and
accretion shock models.

The corona model predicts that NGC 4945 and the Circinus galaxy are Seyfert
galaxies that are among the most promising for neutrino detection, and the
all-flavor neutrino flux reaches
$E_{\nu}F_{E_{\nu}}\sim{10}^{-8}-{10}^{-7}~{}{\rm GeV}~{}{\rm cm}^{-2}~{}{\rm
s}^{-1}$ as seen in Figure 2. They are promising targets for northern-
hemisphere neutrino detectors such as KM3Net (Adrian-Martinez et al., 2016),
Baikal-GVD (Avrorin et al., 2020), P-ONE (Agostini et al., 2020) and Trident
(Ye et al., 2022). For example, the sensitivity of KM3Net is
$E_{\nu}F_{E_{\nu}}\sim{10}^{-9}~{}{\rm GeV}~{}{\rm cm}^{-2}~{}{\rm s}^{-1}$.
We also encourage IceCube searches for starting tracks and showers in the
southern sky666Aartsen et al. (2020a) gave
$E_{\nu}F_{E_{\nu}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}2\times{10}^{-8}~{}{\rm GeV}~{}{\rm
cm}^{-2}~{}{\rm s}^{-1}$ for an $E_{\nu}^{-2}$ spectrum but it is not
sensitive to energies below $10-30$ TeV due to the severe muon background..

##  5 Summary

We studied $\gamma$-ray emission from three nearby Seyfert galaxies, NGC 4151,
NGC 4945 and the Circinus Galaxy. The theoretical models predict that the
neutrino flux is approximately proportional to the intrinsic X-ray luminosity.
NGC 1068 and NGC 4151 are the most promising neutrino sources in the IceCube
sky, while NGC 4945 and Circinus are the best in the southern sky. We analyzed
the 20 MeV – 1 TeV spectra of the three Seyfert galaxies, NGC 4151, NGC 4945
and Circinus, using 14.4 years of the Fermi LAT data. We derived upper limits
on $\gamma$-ray emission from NGC 4151, which gives strong constraints on the
models of high-energy neutrino emission from this source. The recent hint of a
neutrino excess toward NGC 4151, reported by the IceCube Collaboration (Glauch
et al., 2023; Liu et al., 2023), appears to be consistent with the
magnetically-powered corona model. Future observations of NGC 4151 and other
Compton thin AGNs with next-generation neutrino detectors such as IceCube-Gen2
are crucial for discriminating the different models. We also found evidence of
sub-GeV $\gamma$-ray emission in the direction of NGC 4945. If it originates
from NGC 4945, a component other than pionic $\gamma$-ray emission is needed.
Intriguingly, this excess is consistent with CR-induced electromagnetic
cascades expected in the magnetically-powered corona model, although other
interpretations such as high-energy emission from jets would also be possible.
Cascade signatures expected in the MeV and sub-GeV $\gamma$-ray bands are
important for future MeV $\gamma$-ray detectors such as AMEGO-X (Caputo et
al., 2022) e-ASTROGAM (De Angelis et al., 2017), and GRAMS (Aramaki et al.,
2020).

We thank the Topical Workshop: NGC 1068 as cosmic laboratory sponsored by
SFB1258 and Cluster of Excellence ORIGINS, at which we started this project.
The work of K.M. is supported by the NSF grants No. AST-1908689, No.
AST-2108466 and No. AST-2108467, and KAKENHI No. 20H01901 and No. 20H05852.
CMK’s research was supported by an appointment to the NASA Postdoctoral
Program at NASA Goddard Space Flight Center, administered by Oak Ridge
Associated Universities under contract with NASA. S.S.K. acknowledges the
support by KAKENHI No. 22K14028, No. 21H04487, No. 23H04899, and the Tohoku
Initiative for Fostering Global Researchers for Interdisciplinary Sciences
(TI-FRIS) of MEXT’s Strategic Professional Development Program for Young
Researchers. The Fermi LAT Collaboration acknowledges generous ongoing support
from a number of agencies and institutes that have supported both the
development and the operation of the LAT as well as scientific data analysis.
These include the National Aeronautics and Space Administration and the
Department of Energy in the United States, the Commissariat à l’Energie
Atomique and the Centre National de la Recherche Scientifique / Institut
National de Physique Nucléaire et de Physique des Particules in France, the
Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in
Italy, the Ministry of Education, Culture, Sports, Science and Technology
(MEXT), High Energy Accelerator Research Organization (KEK) and Japan
Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg
Foundation, the Swedish Research Council and the Swedish National Space Board
in Sweden. Additional support for science analysis during the operations phase
is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy
and the Centre National d’Études Spatiales in France. This work performed in
part under DOE Contract DE-AC02-76SF00515. Work at NRL is supported by NASA.

##  Appendix A Details of $\gamma$-ray Data Analyses

###  A.1 NGC 4151

NGC 4151 is a nearby Seyfert 1.5 galaxy, located at a distance of $d=15.8$ Mpc
(Yuan et al., 2020). X-ray observations of NGC 4151 indicate that it contains
an ultra-fast outflow (UFO) (Tombesi et al., 2010, 2012). Evidence for high-
energy $\gamma$-ray emission from UFOs, including NGC 4151, was first reported
in Ajello et al. (2021), based on Fermi-LAT observations. In that work, a
stacking analysis was performed using a small sub-sample of the UFO
population, resulting in a $5.1\sigma$ detection. Although NGC 4151 itself was
below the LAT detection threshold, it had an individual significance of
4.2$\sigma$, which was the highest of the sample. It should be noted that when
excluding NGC 4151 from the UFO sample, the population was still detected at
the level of 3.5$\sigma$, with similar spectral properties. This showed that
the signal was not due to NGC 4151 alone. It should also be noted that in the
stacking analysis the source locations are fixed at their optical positions.

Following the detection of the sub-threshold UFO sample, the individual
detection of NGC 4151 was recently reported in Peretti et al. (2023).
Importantly, NGC 4151 is located at a distance of only 0.18∘ from the fourth
most significant hot spot in the search for northern neutrino point sources
(Abbasi et al., 2022), and the $\gamma$-ray detection could potentially have
significant implications for the nature of the IceCube neutrinos.

The LAT observations of NGC 4151 are complicated by the presence of two nearby
BL Lacs. In fact, one of these sources (4FGL J1211.6+3901, located 0.43∘ away)
was carefully analyzed in Peretti et al. (2023), and it was determined that
the detection of NGC 4151 was indeed robust against the presence of this
nearby $\gamma$-ray-bright blazar. The point source model in Peretti et al.
(2023) was originally based on the 4FGL-DR2 catalog; however, in the most
recent catalog (4FGL-DR3) there is now another source (4FGL 1210.3+3928) which
is located even closer, at a separation of 0.08∘, associated with the BL Lac,
1E 1207.9+3945. It therefore seems that the $\gamma$-ray emission thought to
be coming from NGC 4151 is in fact dominated by a nearby blazar.

Figure 4: NGC 4945 results for 3 different ROI sizes, as specified in the
legends. Results from the 4FGL-DR3 are also overlaid. The left panel shows the
SED, the middle panel shows the corresponding significance ($\sqrt{TS}$), and
the right panel shows the count residuals. In addition to the nominal binning,
the first energy bin is divided into two smaller bins, shown with open markers
corresponding to the ROI size. This is done to make a more direct comparison
to the 4FGL-DR3, as discussed in the text. In the middle panel we indicate the
predicted significance based on the additional exposure time, for both the
noise-limited case ($\sqrt{t}$), and the signal-limited case ($t$). The error
bars on the data points are 1$\sigma$, and the ULs are plotted for bins with
TS$<$4, and they are shown at the $95\%$ C.L.

Figure 5: Spatial residuals for the full energy range (top row) and first
energy bin (bottom row). From left to right, the columns are $10^{\circ}$,
$20^{\circ}$, and $26^{\circ}$ ROIs, respectively. Red, cyan, and green
contours are shown at the $\pm 3\sigma,\pm 4\sigma,\mathrm{and}\pm 5\sigma$
levels.

Based on the localization results (presented in Section 3), we keep the nearby
source (4FGL 1210.3+3928) in the model. We add another source at the optical
center of NGC 4151 and refit the data. Unsurprisingly, NGC 4151 is no longer
detected (TS $\approx$ 0), and so we calculate upper limits (ULs), shown in
the left panel of Figure 2, with corresponding data provided in Table 2. The
bins have TS=0, and so we calculate the ULs using a Bayesian approach (as
opposed to a frequentist approach). Specifically, we use the calc_int method
from the IntegralUpperLimit class, available in the Fermi Science Tools. This
method calculates an integral upper limit by integrating the likelihood
function up to a point which contains a given fraction of the total
probability. We verified that the field is well-modeled, with good data-model
agreement seen in the fractional count residuals and spatial residuals. We
note that in a revised version of their work, Peretti et al. (2023) exclude on
the basis of SED modeling that 1E 1207.9+3945 contributes to the emission of
the source they associate to NGC 4151. In our work (and in the upcoming 4FGL-
DR4) we clearly find that the position of the $\gamma$-ray source 4FGL
1210.3+3928 is compatible with 1E 1207.9+394.

Table 2: NGC 4151 SED Data  Energy | Energy Low | Energy High | 95% Flux UL ($\times 10^{-7}$) | TS
---|---|---|---|---
[GeV] | [GeV] | [GeV] | [$\mathrm{MeV\ cm^{-2}\ s^{-1}}$] |
$0.0447$ | $0.02$ | $0.1$ | 10.3 | 0.0
$0.173$ | $0.1$ | $0.3$ | 2.6 | 0.0
$0.548$ | $0.3$ | $1.0$ | 1.1 | 0.0
$1.78$ | $1.0$ | $3.16$ | 0.7 | 0.0
$5.62$ | $3.16$ | $10$ | 1.2 | 0.0
$17.8$ | $10$ | $31.6$ | 0.8 | 0.0
$56.2$ | $31.6$ | $100$ | 3.0 | 0.0
$316.0$ | $100$ | $1000$ | 5.3 | 0.0

###  A.2 NGC 4945

We calculate the SED for NGC 4945, located at the distance $d=3.6$ Mpc
(Monachesi et al., 2016), using the same fitting approach as was used for NGC
4151 (in this case there is no need to re-localize). Results for this are
shown in the middle panel of Figure 2, and the corresponding data are provided
in Table 3. The source is detected with a TS=830. Overall, the SED is
consistent with the values from the 4FGL-DR3. However, compared to the 4FGL,
the two lowest energy bins now have significant signals (TS$>$9), as opposed
to ULs. For the second energy bin, this is quite reasonable given the
additional exposure. Interestingly, extending the lowest energy bin down to 20
MeV also results in a significant signal.

The low-energy detection of NGC 4945 has important implications for the corona
model, and so we perform additional tests to check the robustness of the
emission. The 68% containment angle at 20 MeV for PSF3 events is $\sim
13^{\circ}$, and the 95% containment angle is $\sim 30^{\circ}$. However, our
ROI size is only $10^{\circ}\times 10^{\circ}$, and the model includes all
sources within $15^{\circ}\times 15^{\circ}$. Although the ROI is smaller than
the PSF at 20 MeV, it should be sufficient for the analysis, given the
following consideration. The lowest energy bin spans the energy range $20-100$
MeV, and below 100 MeV the effective area quickly falls off. Specifically, the
(total) effective area at 100 MeV is $\sim 4.4\times$ higher compared to 20
MeV. More of the events will therefore be coming towards the upper part of the
bin, where the 68% containment angle is $\sim 3.3^{\circ}$. Moreover, for a
Gaussian PSF, the distribution peaks towards the center. On the other hand,
the behaviour of the LAT PSF below 50 MeV has not been well characterized. We
therefore test the robustness of the emission with respect to the ROI size.

We repeat the analysis for NGC 4945 using two different regions, motivated by
the approximate $68\%$ and $95\%$ containment radii. For one test we use a
$20^{\circ}\times 20^{\circ}$ ROI, and include all sources within a
$40^{\circ}\times 40^{\circ}$ region. For another test we use a
$26^{\circ}\times 26^{\circ}$ ROI, and include all sources within a
$60^{\circ}\times 60^{\circ}$ region. Note that for the latter test, the data
cover the entire $68\%$ containment region, and the model sources cover the
entire $95\%$ containment region. Figure 4 shows the resulting SEDs,
significance ($\sqrt{\mathrm{TS}}$), and fractional count residuals for all
three fits. Additionally, we overlay the SED and corresponding TS from the
4FGL-DR3. We find consistent results for all three fit variations. The flux in
the first bin is compatible with the $95\%$ ULs from the 4FGL-DR3. The flux
for the second energy bin is slightly lower for the $10^{\circ}$ ROI, but it
is still consistent with the other fits within uncertainties. The fractional
count residuals show that the data are well modeled, with agreement better
than $5\%$ for the full energy range. The data are specifically well modeled
in the first two energy bins. The spatial residuals are shown in Figure 5
(calculated with the fermipy tool, residmap), for both the entire energy range
(top row) and the first energy bin (bottom row). For the full energy range,
the $10^{\circ}$ ROI is very well modeled, whereas the larger ROIs show some
excesses/deficits towards the edges of the field. For the first energy bin,
all fit variations show good data-model agreement.

In order to better understand the reason for the discrepancy in the first
energy bin with respect to the 4FGL-DR3, we divide it into two smaller bins to
match the 4FGL binning, and re-calculated the SED. The lower part of the bin
spans the energy range $20-45$ MeV, and the upper part of the bin spans the
energy range $45-100$ MeV. We make the split at 45 MeV (compared to 50 MeV for
the lower edge of the 4FGL bin), as required by our initial energy binning.
This test is performed for all three ROI variations, and the results are shown
in Figure 4 with open markers, corresponding to the respective ROI size.

For our nominal ROI size of $10^{\circ}$ we find a significance
($\sqrt{\mathrm{TS}}$) of 5.8 and 2.6 for the two bins, respectively. The
significance of the upper part of the bin is now reasonable given that of the
4FGL, and the flux is also compatible with the 4FGL upper limits.
Additionally, we find a significant signal in the lower portion of the bin.
However, the corresponding count residuals in the lower portion of the bin are
slightly elevated at the level of $\sim 5\%$. This can plausibly be attributed
to some mis-modeling of either the isotropic emission or the Galactic diffuse
emission (which at these energies is dominated by inverse-Compton radiation).
Note that the isotropic model uses a power law extrapolation below 34 MeV, and
the Galactic diffuse model uses a power law extrapolation below 50 MeV. For
the $20^{\circ}$ ROI, the flux of the two smaller bins is found to be
consistent with that of the combined bin, and likewise the count residuals are
close to zero. For the $26^{\circ}$ ROI, the upper portion of the bin is also
compatible with the combined bin, however, the lower portion of the bin is
over-modeled by $\sim 3\%$, and the SED calculation results in only an upper
limit, albeit at the same level as the observed signal.

As another test, we have made the SED calculation using a spectral index of
3.0 (instead of the commonly used value of 2.0), consistent with our model
predictions, and also freeing all point sources within $5^{\circ}$, having a
TS$>$16\. This test was ran for the $20^{\circ}$ ROI. We have found that the
low-energy bin ($E_{\gamma}<100$ MeV) is still detected with a similar flux,
having a TS of 19.0. Finally, we have tested the SED calculation when freeing
all sources within $7^{\circ}$, having a spectral index $>2.5$. In this case,
we have found consistent results for the low-energy bin, having a TS of 27.5.

Overall, these tests show that the emission between $20-100$ MeV appears to be
robust, although there is clearly a systematic uncertainty in the measured
flux of a factor of $\sim 2-3$. However, this result should still be
interpreted with caution, due to the generally poor performance of the LAT
below 50 MeV.

Table 3: NGC 4945 SED Data  Energy | Energy Low | Energy High | Flux ($\times 10^{-7}$) | 1$\sigma$ Flux Error ($\times 10^{-7}$) | 95% Flux U.L ($\times 10^{-7}$) | TS
---|---|---|---|---|---|---
[GeV] | [GeV] | [GeV] | [$\mathrm{MeV\ cm^{-2}\ s^{-1}}$] | [$\mathrm{MeV\ cm^{-2}\ s^{-1}}$] | [$\mathrm{MeV\ cm^{-2}\ s^{-1}}$] |
$0.0447$ | $0.02$ | $0.1$ | 31.3 | 6.5 | 42.1 | 23.0
$0.173$ | $0.1$ | $0.3$ | 7.8 | 2.1 | 11.3 | 14.0
$0.548$ | $0.3$ | $1.0$ | 16.7 | 1.2 | 18.7 | 251
$1.78$ | $1.0$ | $3.16$ | 11.1 | 0.9 | 12.6 | 270
$5.62$ | $3.16$ | $10$ | 7.6 | 1.0 | 9.3 | 171
$17.8$ | $10$ | $31.6$ | 4.7 | 1.1 | 6.9 | 73.3
$56.2$ | $31.6$ | $100$ | 4.1 | 1.9 | 7.9 | 21.1
$316.0$ | $100$ | $1000$ | 1.5 | 1.6 | 5.8 | 7.6

###  A.3 Circinus

We also calculate the SED for Circinus, located at $d=4.2$ Mpc (Freeman et
al., 1977), and results are shown in the right panel of Figure 2, and the
corresponding data are provided in Table 4. The source is detected with a
TS=130. Overall, the SED is consistent with the values from the 4FGL-DR3. The
fractional count residuals are very good, being close to zero for all
energies.

Table 4: Circinus SED Data  Energy | Energy Low | Energy High | Flux ($\times 10^{-7}$) | 1$\sigma$ Flux Error ($\times 10^{-7}$) | 95% Flux U.L ($\times 10^{-7}$) | TS
---|---|---|---|---|---|---
[GeV] | [GeV] | [GeV] | [$\mathrm{MeV\ cm^{-2}\ s^{-1}}$] | [$\mathrm{MeV\ cm^{-2}\ s^{-1}}$] | [$\mathrm{MeV\ cm^{-2}\ s^{-1}}$] |
$0.0447$ | $0.02$ | $0.1$ | 22.8 | 13.9 | 45.8 | 2.7
$0.173$ | $0.1$ | $0.3$ | 6.5 | 4.0 | 13.1 | 2.7
$0.548$ | $0.3$ | $1.0$ | 6.8 | 1.6 | 9.5 | 18.6
$1.78$ | $1.0$ | $3.16$ | 5.2 | 1.0 | 6.9 | 32.5
$5.62$ | $3.16$ | $10$ | 4.3 | 0.9 | 5.8 | 40.6
$17.8$ | $10$ | $31.6$ | 2.6 | 1.0 | 4.5 | 15.6
$56.2$ | $31.6$ | $100$ | 2.7 | 1.5 | 5.8 | 14.7
$316.0$ | $100$ | $1000$ | 2.3 | 2.0 | 7.0 | 3.5

## References

  * Aartsen et al. (2013a) Aartsen, M., et al. 2013a, Phys.Rev.Lett., 111, 021103
  * Aartsen et al. (2013b) —. 2013b, Science, 342, 1242856
  * Aartsen et al. (2015) Aartsen, M. G., et al. 2015, Astrophys. J., 809, 98
  * Aartsen et al. (2020a) —. 2020a, Phys. Rev. Lett., 125, 121104
  * Aartsen et al. (2020b) —. 2020b, Phys. Rev. Lett., 124, 051103
  * Aartsen et al. (2021) —. 2021, J. Phys. G, 48, 060501
  * Abbasi et al. (2022) Abbasi, R., et al. 2022, Science, 378, 538.  https://www.science.org/doi/10.1126/science.abg3395
  * Abdollahi et al. (2020) Abdollahi, S., et al. 2020, Astrophys. J. Suppl., 247, 33
  * Abdollahi et al. (2022) —. 2022, Astrophys. J. Supp., 260, 53
  * Acciari et al. (2019) Acciari, V. A., et al. 2019, Astrophys. J., 883, 135
  * Ackermann et al. (2015) Ackermann, M., et al. 2015, Astrophys.J., 799, 86
  * Adrian-Martinez et al. (2016) Adrian-Martinez, S., et al. 2016, J. Phys., G43, 084001
  * Agostini et al. (2020) Agostini, M., et al. 2020, Nature Astron., 4, 913
  * Aguilar-Ruiz et al. (2021) Aguilar-Ruiz, E., Fraija, N., Joshi, J. C., Galvan-Gamez, A., & de Diego, J. A. 2021, Phys. Rev. D, 104, 083013
  * Ajello et al. (2023) Ajello, M., Murase, K., & McDaniel, A. 2023, Astrophys. J. Lett., 954, L49
  * Ajello et al. (2021) Ajello, M., Baldini, L., Ballet, J., et al. 2021, Astrophys. J., 921, 144
  * Anchordoqui et al. (2021) Anchordoqui, L. A., Krizmanic, J. F., & Stecker, F. W. 2021, PoS, ICRC2021, 993
  * Aramaki et al. (2020) Aramaki, T., Hansson Adrian, P., Karagiorgi, G., & Odaka, H. 2020, Astropart. Phys., 114, 107
  * Atwood et al. (2013) Atwood, W., Albert, A., Baldini, L., et al. 2013, arXiv e-prints, arXiv:1303.3514
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, Astrophys. J., 697, 1071
  * Avrorin et al. (2020) Avrorin, A. D., et al. 2020, PoS, ICRC2019, 1011
  * Bacchini et al. (2022) Bacchini, F., Arzamasskiy, L., Zhdankin, V., et al. 2022, Astrophys. J., 938, 86
  * Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1
  * Becker et al. (2006) Becker, P. A., Le, T., & Dermer, C. D. 2006, Astrophys. J., 647, 539
  * Bentz et al. (2022) Bentz, M. C., Williams, P. R., & Treu, T. 2022, Astrophys. J., 934, 168
  * Bruel et al. (2018) Bruel, P., Burnett, T. H., Digel, S. W., et al. 2018, arXiv e-prints, arXiv:1810.11394
  * Caputo et al. (2022) Caputo, R., et al. 2022, arXiv:2208.04990
  * Comisso & Sironi (2019) Comisso, L., & Sironi, L. 2019, Astrophys. J., 886, 122
  * Comisso & Sironi (2022) —. 2022, Astrophys. J. Lett., 936, L27
  * De Angelis et al. (2017) De Angelis, A., et al. 2017, Exper. Astron., 44, 25
  * Drury (1983) Drury, L. O. 1983, Rept. Prog. Phys., 46, 973
  * Eichmann et al. (2022) Eichmann, B., Oikonomou, F., Salvatore, S., Dettmar, R.-J., & Becker Tjus, J. 2022, Astrophys. J., 939, 43
  * Fiorillo et al. (2023) Fiorillo, D. F. G., Petropoulou, M., Comisso, L., Peretti, E., & Sironi, L. 2023, arXiv:2310.18254
  * Freeman et al. (1977) Freeman, K. C., Karlsson, B., Lynga, G., et al. 1977, A&A , 55, 445
  * Glauch et al. (2023) Glauch, T., et al. 2023, PoS, ICRC2023, 1052
  * Greenhill et al. (1997) Greenhill, L. J., Moran, J. M., & Herrnstein, J. R. 1997, Astrophys. J. Lett., 481, L23
  * Greenhill et al. (2003) Greenhill, L. J., Booth, R. S., Ellingsen, S. P., et al. 2003, Astrophys. J., 590, 162
  * Guo et al. (2020) Guo, F., Liu, Y.-H., Li, X., et al. 2020, Phys. Plasmas, 27, 080501
  * Hoshino (2015) Hoshino, M. 2015, Phys. Rev. Lett., 114, 061101
  * Hoshino & Lyubarsky (2012) Hoshino, M., & Lyubarsky, Y. 2012, Space Sci. Rev. , 173, 521
  * Inoue et al. (2022) Inoue, S., Cerruti, M., Murase, K., & Liu, R.-Y. 2022, arXiv:2207.02097
  * Inoue & Khangulyan (2023) Inoue, Y., & Khangulyan, D. 2023, arXiv:2304.04138
  * Inoue et al. (2020) Inoue, Y., Khangulyan, D., & Doi, A. 2020, Astrophys. J. Lett., 891, L33
  * Kheirandish et al. (2021) Kheirandish, A., Murase, K., & Kimura, S. S. 2021, Astrophys. J., 922, 45
  * Kimura et al. (2015) Kimura, S. S., Murase, K., & Toma, K. 2015, Astrophys.J., 806, 159
  * Kimura et al. (2016) Kimura, S. S., Toma, K., Suzuki, T. K., & Inutsuka, S.-i. 2016, Astrophys. J., 822, 88
  * Kimura et al. (2019) Kimura, S. S., Tomida, K., & Murase, K. 2019, Mon. Not. Roy. Astron. Soc., 485, 163
  * Koss et al. (2022) Koss, M. J., et al. 2022, Astrophys. J. Suppl., 261, 1
  * Kunz et al. (2016) Kunz, M. W., Stone, J. M., & Quataert, E. 2016, Phys. Rev. Lett., 117, 235101
  * Liu et al. (2023) Liu, Q., et al. 2023, 18th International Conference on Topics in Astroparticle & Underground Physics (TAUP 2023).  https://indico.cern.ch/event/1199289/contributions/5447243
  * Mbarek et al. (2023) Mbarek, R., Philippov, A., Chernoglazov, A., Levinson, A., & Mushotzky, R. 2023, arXiv:2310.15222
  * Monachesi et al. (2016) Monachesi, A., Bell, E. F., Radburn-Smith, D., et al. 2016, Mon. Not. Roy. Astron. Soc., 457, 1419
  * Murase (2018) Murase, K. 2018, Phys. Rev., D97, 081301
  * Murase (2022) —. 2022, Astrophys. J. Lett., 941, L17
  * Murase et al. (2016) Murase, K., Guetta, D., & Ahlers, M. 2016, Phys. Rev. Lett., 116, 071101
  * Murase et al. (2020) Murase, K., Kimura, S. S., & Mészáros, P. 2020, Phys. Rev. Lett., 125, 011101\.  https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.125.011101
  * Murase & Stecker (2022) Murase, K., & Stecker, F. W. 2022, arXiv:2202.03381
  * Neronov et al. (2023) Neronov, A., Savchenko, D., & Semikoz, D. V. 2023, arXiv:2306.09018
  * Peretti et al. (2023) Peretti, E., Peron, G., Tombesi, F., et al. 2023, arXiv:2303.03298
  * Puccetti et al. (2014) Puccetti, S., et al. 2014, Astrophys. J., 793, 26
  * Stecker et al. (1991) Stecker, F. W., Done, C., Salamon, M. H., & Sommers, P. 1991, Phys.Rev.Lett., 66, 2697
  * Sun & Bai (2021) Sun, X., & Bai, X.-N. 2021, Mon. Not. Roy. Astron. Soc., 506, 1128
  * Tanimoto et al. (2019) Tanimoto, A., Ueda, Y., Odaka, H., et al. 2019, arXiv:1904.08945
  * Tombesi et al. (2012) Tombesi, F., Cappi, M., Reeves, J., & Braito, V. 2012, MNRAS, 422, L1
  * Tombesi et al. (2010) Tombesi, F., Cappi, M., Reeves, J., et al. 2010, A&A, 521, A57
  * Ye et al. (2022) Ye, Z. P., et al. 2022, arXiv:2207.04519
  * Yuan et al. (2020) Yuan, W., Fausnaugh, M. M., Hoffmann, S. L., et al. 2020, Astrophys. J., 902, 26
  * Zhang & Murase (2023) Zhang, B. T., & Murase, K. 2023, Mon. Not. Roy. Astron. Soc., 524, 76
  * Zhdankin et al. (2019) Zhdankin, V., Uzdensky, D. A., Werner, G. R., & Begelman, M. C. 2019, Phys. Rev. Lett., 122, 055101