#  High-energy electromagnetic, neutrino, and cosmic-ray emission by stellar-
mass black holes in disks of active galactic nuclei

Hiromichi Tagawa11affiliationmark:  22affiliationmark:  33affiliationmark:  ,
Shigeo S Kimura22affiliationmark:  44affiliationmark:  , Zoltán
Haiman11affiliationmark:  11affiliationmark:  Department of Astronomy,
Columbia University, 550 W. 120th St., New York, NY, 10027, USA
22affiliationmark:  Astronomical Institute, Graduate School of Science, Tohoku
University, Aoba, Sendai 980-8578, Japan
33affiliationmark:  National Astronomical Observatory of Japan, National
Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
44affiliationmark:  Frontier Research Institute for Interdisciplinary
Sciences, Tohoku University, Sendai 980-8578, Japan
E-mail<EMAIL_ADDRESS>
###### Abstract

Some Seyfert galaxies are detected in high-energy gamma rays, but the
mechanism and site of gamma-ray emission are unknown. Also, the origins of the
cosmic high-energy neutrino and MeV gamma-ray backgrounds have been veiled in
mystery since their discoveries. We propose emission from stellar-mass BHs
(sBHs) embedded in disks of active galactic nuclei (AGN) as their possible
sources. These sBHs are predicted to launch jets due to the Blandford-Znajek
mechanism, which can produce intense electromagnetic, neutrino, and cosmic-ray
emissions. We investigate whether these emissions can be the sources of cosmic
high-energy particles. We find that emission from internal shocks in the jets
can explain gamma rays from nearby radio-quiet Seyfert galaxies including
NGC1068, if the Lorentz factor of the jets ($\Gamma_{\rm j}$) is high. On the
other hand, for moderate $\Gamma_{\rm j}$, the emission can significantly
contribute to the background gamma-ray and neutrino intensities in the
$\sim{\rm MeV}$ and $\lesssim{\rm PeV}$ bands, respectively. Furthermore, for
moderate $\Gamma_{\rm j}$ with efficient amplification of the magnetic field
and cosmic-ray acceleration, the neutrino emission from NGC1068 and the
ultrahigh-energy cosmic rays can be explained. These results suggest that the
neutrino flux from NGC1068 as well as the background intensities of ${\rm
MeV}$ gamma rays, neutrinos, and the ultrahigh-energy cosmic rays can be
explained by a unified model. Future MeV gamma-ray satellites will test our
scenario for neutrino emission.

###### Subject headings:

Stellar mass black holes (1611), Active galactic nuclei (16), Accretion (14),
Black hole physics (159), Jets (870), Galactic center (565)

##  1\. Introduction

Our Universe is filled with high-energy particles, including charged particles
(cosmic rays; CRs), neutrinos, and gamma rays, but the origins of these cosmic
high-energy particles are unknown. Recently, significant progress have been
made in high-energy neutrino astrophysics. IceCube reported the detection of
extraterrestrial neutrinos in 2013 (Aartsen et al., 2013), and has been
improving the measurement of the cosmic high-energy neutrino background for 10
TeV–100 PeV (Aartsen et al., 2015, 2020a). They also identified a nearby
Seyfert galaxy, NGC 1068, as a cosmic neutrino source (Aartsen et al., 2020b;
IceCube Collaboration et al., 2022). In addition, they reported a hint of
association between neutrino signals and radio-quiet active galactic nuclei
(AGN) (Abbasi et al., 2022). However, possible mechanisms for CR acceleration
and subsequent neutrino production sites in radio quiet AGN remain unclear.
This has motivated investigations of non-thermal phenomena in hot accretion
flows (Kimura et al., 2021b; Gutiérrez et al., 2021), hot coronae (Murase et
al., 2020; Eichmann et al., 2022), accretion shocks (Inoue et al., 2020), disk
winds (Inoue et al., 2022), and jets from accreting binaries (Sridhar et al.,
2022).

High-energy gamma rays ($>100$ MeV) have also been detected from radio-quiet
AGNs (Wojaczyński et al., 2015), but the origin of these gamma-rays are also
controversial. The majority of these gamma-ray detected AGNs exhibit
signatures of starburst activity, which causes gamma-ray production via
hadronuclear interactions (e.g., Ajello et al., 2020). The neutrino-emitting
AGN, NGC 1068, is also detected in high-energy gamma rays and shows starburst
activity. On the other hand, the gamma-ray signals from NGC 4945 show spectral
variation correlated with its X-ray flux, which implies that gamma rays may be
associated with AGN activity (Wojaczyński & Niedźwiecki, 2017). AGN wind
interacting with the dusty torus can be a possible site of high-energy gamma-
ray production (Inoue et al., 2022).

Finally, the origins of the unresolved cosmic MeV gamma-ray background and
ultrahigh-energy cosmic rays (UHECRs) have been unknown for a long time
(Inoue, 2014; Alves Batista et al., 2019). Radio-quiet AGNs have been proposed
as candidates for both the cosmic MeV gamma-ray background (e.g., Inoue et
al., 2013; Kimura et al., 2021b) and a source of UHECRs (Pe’Er et al., 2009).

Here, we propose a novel scenario for high-energy emission from radio-quiet
AGNs, where we consider relativistic jets launched from stellar-mass black
holes (sBHs) embedded in AGN disks. It has been predicted that stars and
compact objects (COs) including stellar-mass black holes (sBHs) are embedded
in AGN disks due to capture by dynamical interactions of nuclear star clusters
(Miralda-Escude & Gould 2000; Lu et al. 2013) with the AGN disk (Ostriker,
1983; Syer et al., 1991) and in-situ star formation (Levin & Beloborodov,
2003; Goodman & Tan, 2004; Thompson et al., 2005; Levin, 2007). There are
several observations supporting this picture (Artymowicz et al., 1993; Levin &
Beloborodov, 2003; Milosavljević & Loeb, 2004; Tagawa et al., 2020b).
Recently, the evolution of COs in AGN disks has attracted significant
attention as these are promising environments for some of the sBH-sBH (e.g.
Bartos et al., 2017; Stone et al., 2017; McKernan et al., 2018; Yang et al.,
2019; Tagawa et al., 2020b) and sBH-neutron star (NS) mergers (McKernan et
al., 2020; Tagawa et al., 2021) reported as gravitational wave (GW) events by
LIGO/Virgo/KAGRA (Abbott et al., 2021; The LIGO Scientific Collaboration et
al., 2021).

Many recent studies in the wake of the LIGO/Virgo/KAGRA discoveries have
investigated emission from transients emerging in AGN disks (McKernan et al.,
2019; Graham et al., 2020; Perna et al., 2021a, b; Zhu et al., 2021c, a, b;
Yang et al., 2022b; Moranchel-Basurto et al., 2021; Grishin et al., 2021;
Kimura et al., 2021a; Wang et al., 2021b; Yuan et al., 2022). Closest to the
present study, Wang et al. (2021a) considered emission from shocks emerging
from interactions of Blandford-Znajek (BZ) jets, (Blandford & Znajek, 1977)
launched from accreting sBHs, with gas in the AGN’s broad line region. In
Tagawa et al. (2022) (Paper I), we estimated the structure of the cavity
created by the BZ jet and the dynamical evolution of gas around the BHs. In
Tagawa et al. (2023a) (Paper II), we investigated the properties of emission
released when a jet, launched from merging sBHs embedded in an AGN disk,
breaks out from the disk. In this paper, we consider in more detail the high-
energy radiation from jets launched from solitary sBHs in AGN disks. Non-
thermal electrons accelerated at the internal shock emit broadband
electromagnetic radiation via synchrotron and inverse Compton scattering. Non-
thermal protons, accelerated together with electrons, will produce neutrinos
via hadronic interactions and might become cosmic rays after escaping from the
system. We evaluate the possibility that gamma rays and neutrinos from nearby
Seyfert galaxies, including NGC1068, are produced in such jets. We also
estimate their contributions to the diffuse cosmic gamma-ray, neutrino, and
UHECR background intensities.

Figure 1.—  A schematic picture of emission from internal shocks in a jet
launched from an sBH accreting gas in an AGN disk.

##  2\. Model

To assess the observability of electromagnetic (EM), neutrino, and CR emission
produced in shocks around the BZ jets, we first model the properties of the
jets launched from rapidly accreting and spinning sBHs in AGN disks (see Fig.
1 for a schematic illustration). The jet is predicted to be launched from an
accreting sBH embedded in an AGN disk (Appendix A.1). We focus on sBHs, which
are expected to be common in AGN disks and presumably dominate the total jet
luminosity, and we adopt the same values for our model parameters as in the
fiducial model of Paper I. For these parameters, through Eq. (A1) and Eq. (1)
of Paper I, the jet power of $L_{\rm j}\sim 2\times 10^{42}~{}{\rm erg/s}$ is
derived, which we adopt as a fiducial value.

In the fiducial model (model M1), we set the fraction of postshock energy
carried by the post-shock magnetic field and by electrons and protons to
$\epsilon_{\rm B}=0.01$, $\epsilon_{e}=0.02$, $\epsilon_{\rm CR}=0.05$ (e.g.,
Panaitescu & Kumar 2001; Santana et al. 2014; Spitkovsky 2008; Sironi et al.
2013; Tomita et al. 2019), respectively; the power-law slope for injected
electrons and protons accelerated by the first order Fermi process to $p=2.5$;
the Lorentz factor of the jet to $\Gamma_{\rm j}=30$ as derived in Appendix
A.2; the variability timescale of the jet to $T_{\rm vari}=10^{-3}~{}{\rm s}$;
and the opening angle of the injected jet to $\theta_{\rm j}=0.2$ (e.g.
Pushkarev et al., 2009; Hada et al., 2013, 2018; Berger, 2014). Here note that
the parameters $\epsilon_{\rm B}$, $\epsilon_{e}$, $\epsilon_{\rm CR}$, $p$,
and $\Gamma_{\rm j}$ are highly uncertain and expected to be distributed in
wide ranges of values as $\epsilon_{\rm B}\sim 10^{-5}$–$0.3$,
$\epsilon_{e}\sim 10^{-2}$–$0.5$, $\epsilon_{\rm CR}\sim 0.05$–$0.2$, $p\sim
2$–$3$, and $\Gamma_{\rm j}\sim 2$–$100$ depending on sources (e.g. Panaitescu
& Kumar, 2001; Sironi et al., 2013; Santana et al., 2014; Caprioli &
Spitkovsky, 2014; Troja et al., 2019; Matsumoto et al., 2020; Caprioli et al.,
2020). We also show the results with $\Gamma_{\rm j}=4$, $\epsilon_{e}=0.05$,
$\epsilon_{\rm B}=0.005$, $\epsilon_{\rm CR}=0.05$, and $p=2.2$ (model M2),
and those with $\Gamma_{\rm j}=4$, $\epsilon_{e}=0.05$, $\epsilon_{\rm
B}=0.3$, $\epsilon_{\rm CR}=0.15$, $p=2.2$, and $L_{\rm j}=10^{43}~{}{\rm
erg/s}$ (model M3), as results are sensitive to these parameters (see Appendix
A.5 for the parameter space where our non-thermal emission models are
applicable. ) Note that $L_{\rm j}=2\times 10^{42}~{}{\rm erg/s}$ and $L_{\rm
j}=10^{43}~{}{\rm erg/s}$ adopted in models M1–M2 and M3 are, respectively,
predicted for the jets from sBHs at the distance from the SMBH being $1$ and
$0.01~{}{\rm pc}$ (Paper I), where sBHs are typically accumulated due to a
long migration timescale and gap formation (Tagawa et al., 2020b; Perna et
al., 2021b).

Figure 2.—  The spectral energy distribution for non-thermal emission from
internal shocks of the jet for models M1 (solid), M2 (dashed), and M3
(dotted).  Figure 3.—  The cumulative number of gamma-ray emitting sources
predicted to be detectable by Fermi-LAT as a function of the luminosity
distance for the fiducial model (solid black line), a model with a smaller
number of accreting sBHs per AGN disk ($N_{\rm BH,acc}=5$, dashed orange
line), and that with a weaker gamma-ray luminosity ($L_{\gamma,\rm
iso}=10^{40}\,{\rm erg/s}$, dashed blue line). The number of gamma-ray
sources, inferred from observations, in nearby radio-quiet Seyfert galaxies
for moderate and optimistic cases are shown in dotted red and brown lines.

##  3\. Gamma rays

During the propagation of the jet, its kinetic energy is considered to be
dissipated. We assume that the fraction, $\epsilon_{\rm CR}$ and
$\epsilon_{e}$, of the postshock energy of the jet is used to accelerate
protons and electrons, respectively, via the first-order Fermi process (Bell,
1978; Blandford & Eichler, 1987). Then, the non-thermal electrons emit
broadband radiation via synchrotron and inverse Compton scattering. The
spectral shapes of the non-thermal emission produced at internal shocks of the
jet are presented in Fig. 2 (see Appendix A.3 for their derivation). In the
fiducial model (M1), synchrotron emission and synchrotron-self Compton
scattering, respectively, produce the bright emission in optical–$\rm MeV$ and
in X-ray–$\rm GeV$ bands. In addition, synchrotron self-absorption and
$\gamma\gamma$ annihilation create upper and lower cutoffs, respectively. The
emission by hadronic processes is also significantly absorbed by
$\gamma\gamma$ annihilation as discussed in Section § 5.2. Below we discuss
whether gamma-ray emission from the jets can be detected from nearby Seyfert
galaxies. In addition, we discuss the contribution of the jets to the diffuse
gamma-ray background.

###  3.1. Emission from radio-quiet Seyfert galaxies

We consider the possibility that the gamma rays from radio-quiet galaxies are
originated from the BZ jets launched from sBHs in AGN disks. According to Fig.
2, the isotropic-equivalent gamma-ray luminosity from the jet at $\sim
1$–$10~{}{\rm GeV}$ (Wojaczyński et al., 2015) for model M1 is $L_{\rm GeV,\rm
iso}\sim 3\times 10^{40}~{}{\rm erg~{}s^{-1}}$. Then, the intrinsic gamma-ray
luminosity is

$\displaystyle L_{\rm GeV}=p_{\theta}L_{\rm GeV,iso}=$ $\displaystyle\sim
5\times 10^{38}~{}{\rm erg~{}s^{-1}}\left(\frac{L_{\rm GeV,iso}}{3\times
10^{40}~{}{\rm erg~{}s^{-1}}}\right)\left(\frac{\theta_{\rm
j}}{0.2}\right)^{2},$ (1)

where $p_{\theta}=\theta_{\rm j}^{2}/2$ is the probability that the jet is
directed towards an observer. By using Fermi-LAT with the sensitivity of
$(E_{\gamma}F_{E_{\gamma}})_{\rm LAT}\sim 3\times 10^{-13}~{}{\rm
erg~{}cm^{-2}~{}s^{-1}}$ at $\sim{\rm GeV}$ (Atwood et al., 2009; Funk et al.,
2013), the detectable distance for the gamma rays is

$\displaystyle d_{L,\rm det}=\left(\frac{L_{\rm GeV,\rm
iso}}{4\pi(E_{\gamma}F_{E_{\gamma}})_{\rm sens}}\right)^{1/2}$
$\displaystyle\sim 27\,{\rm Mpc}\left(\frac{L_{\rm GeV,\rm iso}}{3\times
10^{40}\,{\rm erg/s}}\right)^{1/2}\left(\frac{(E_{\gamma}F_{E_{\gamma}})_{\rm
LAT}}{3\times 10^{-13}~{}{\rm erg~{}cm^{-2}s^{-1}}}\right)^{-1/2}.$ (2)

The detectable number of such gamma-ray sources within the distance $d_{L}$ is
roughly estimated as

$\displaystyle N_{\rm det}(d_{L})=p_{\theta}\frac{4\pi[{\rm
min}(d_{L},d_{L,\rm det})]^{3}}{3}n_{\rm rotBH,acc}~{}$ $\displaystyle\sim
8\left(\frac{\theta_{\rm j}}{0.2}\right)^{2}\left(\frac{{\rm
min}(d_{L},d_{L,\rm det})}{27\,{\rm Mpc}}\right)^{3}\left(\frac{n_{\rm
rotBH,acc}}{5\times 10^{-3}\,{\rm Mpc}^{-3}}\right)~{}~{}~{}$ (3)

where $n_{\rm rotBH,acc}=n_{\rm AGN}N_{\rm rotBH,acc}$ is the number density
of rapidly accreting and spinning sBHs in AGN disks, $N_{\rm rotBH,acc}$ is
their typical number in a single AGN disk and $n_{\rm AGN}$ is the AGN space
density. In Eq. (3.1), we assume that the jets (and the sBH spins) are
randomly oriented (Tagawa et al., 2020b, a), and we adopt $n_{\rm AGN}=5\times
10^{-4}\,{\rm Mpc^{-3}}$ considering AGNs with X-ray luminosity of
$L_{X}\gtrsim 10^{42}\,{\rm erg/s}$ (Ueda et al., 2014). Assuming that at any
given time, the active fraction of sBHs is $f_{\rm active}\sim 0.1$ as roughly
estimated in Paper I, we adopt

$\displaystyle N_{\rm rotBH,acc}=N_{\rm rotBH,AGN}f_{\rm active}$
$\displaystyle\sim 10\left(\frac{N_{\rm
rotBH,AGN}}{100}\right)\left(\frac{f_{\rm active}}{0.1}\right),$ (4)

where $N_{\rm rotBH,AGN}\sim 100$ is based on typical numbers found in models
for the AGN-embedded sBH population (Tagawa et al., 2020b, 2021). Although the
distribution of the positions of sBHs in an AGN disk is uncertain and is
influenced by phases and properties of AGNs, we assume that most sBHs are at
$1~{}{\rm pc}$ (models M1 and M2) and $\sim 10\%$ of them are at $0.01~{}{\rm
pc}$ (model M3). Here, the former component corresponds to the sBHs captured
by an AGN disk before considerable migration, while the latter corresponds to
the sBHs accumulated at gaps, which are predicted to form in low opacity
regions of an AGN disk (Thompson et al., 2005; Tagawa et al., 2020b). The low
number density of active sBHs at $\sim 0.01~{}{\rm pc}$ is required so that
the neutrino flux for model M3 does not exceed the background neutrino flux (§
4.2), while if it is lowered further, the probability to detect neutrinos from
NGC1068 by model M3 becomes low (§ 4.3).

Fig. 3 compares the detectable number of gamma-ray sources predicted by the
models to the number of detected gamma-ray sources from nearby radio-quiet
Seyfert galaxies observed by Fermi-LAT. As fiducial candidates (dotted red
line), we include NGC4151, NGC6814, and NGC4258, for which the likelihood
ratio of the non-detection to the detection of gamma rays are $\sim 2\times
10^{-4}$, $\sim 10^{-7}$, and $\sim 9\times 10^{-3}$, respectively. In the
”optimistic” distribution (dotted brown), we additionally include three
Seyfert galaxies with starbursts, Circinus, NGC1068, and NGC4945, as the gamma
rays from these galaxies are not fully explained by starbursts (Hayashida et
al., 2013; Eichmann & Becker Tjus, 2016; Wojaczyński & Niedźwiecki, 2017).
From the figure, we find that the observed distribution of gamma-ray sources
is consistent with our models with $n_{\rm rotBH,acc}\sim 2\times 10^{-3}$ –
$5\times 10^{-3}~{}{\rm Mpc}^{-3}$ and $L_{\rm GeV,\rm iso}\sim 10^{40}$ –
$3\times 10^{40}\,{\rm erg/s}$ (solid black, dashed blue, and dashed orange
lines). Conversely, if the gamma rays from the radio-quiet Seyfert galaxies
are not originated from sBHs (e.g., $N_{\rm det}\sim 0$), $n_{\rm rotBH,acc}$
can be roughly constrained to be

$\displaystyle n_{\rm rotBH,acc}\lesssim 5\times 10^{-4}\,{\rm
Mpc}^{-3}\left(\frac{N_{\rm det}}{1}\right)\left(\frac{\theta_{\rm
j}}{0.2}\right)^{-2}$ $\displaystyle\left(\frac{L_{\rm GeV,iso}}{3\times
10^{40}~{}{\rm
erg/s}}\right)^{-3/2}\left(\frac{(E_{\gamma}F_{E_{\gamma}})_{\rm LAT}}{3\times
10^{-13}~{}{\rm erg~{}cm^{-2}s^{-1}}}\right)^{3/2}.$ (5)

Note that the gamma-ray luminosity at $\sim$ GeV bands is obscured by
$\gamma\gamma$ annihilation for low Lorentz factors (e.g. dashed and dotted
lines with $\Gamma_{\rm j}=4$ in Fig. 2). In this case, the jets from sBHs in
AGN disks cannot explain the GeV gamma rays from Seyfert galaxies, and another
explanation would be required.

Figure 4.—  The contribution to the gamma-ray flux from NGC1068 (panel b) and
the background gamma-ray intensity (panel a) by the internal shocks for models
M1 (solid black), M2 (dashed black), and M3 (dotted black). In panel a, solid
points represent the intensity observed by Fermi-LAT (red, Ackermann et al.
2015), COMPTEL (teal, Weidenspointner et al. 2000), SMM (brown, Watanabe et
al. 1997), Swift BAT (cyan, Ajello et al. 2008), and RXTE (purple, Revnivtsev
et al. 2003). Thin solid, dashed, and dotted brown lines are the intensities
predicted by models for Seyfert galaxies (Gilli et al., 2007), blazars (Giommi
& Padovani, 2015), and star-forming galaxies (Lacki et al., 2014), presented
in De Angelis et al. (2021). In panel b, the red, blue, and purple points are
the observed gamma-ray intensities adopted from Abdollahi et al. (2020),
Ajello et al. (2017), and (Acciari et al., 2019), respectively. The
sensitivities for future telescopes in $\sim{\rm MeV}$ bands are presented by
colored dashed lines.

###  3.2. Cosmic gamma-ray backrground intensity

We next estimate the contribution of gamma-ray emission from the global
population of jets to the diffuse background intensity in MeV bands as

$\displaystyle E_{\gamma}^{2}\Phi_{\gamma}\sim f_{z}L_{\rm
MeV,iso}p_{\theta}n_{\rm rotBH,acc}\frac{c}{4\pi H_{0}}$ $\displaystyle\sim
5\times 10^{-6}{\rm GeV~{}cm^{-2}~{}s^{-1}~{}sr^{-1}}$
$\displaystyle\times\left(\frac{L_{\rm MeV,iso}p_{\theta}}{2\times
10^{40}\,{\rm erg/s}}\right)\left(\frac{n_{\rm rotBH,acc}}{5\times
10^{-3}\,{\rm Mpc^{-3}}}\right)\left(\frac{f_{z}}{2}\right)$ (6)

(Fig. 4a), where $L_{\rm MeV,iso}$ is the isotropic-equivalent gamma-ray
luminosity around MeV bands, $H_{0}=67.8\,{\rm km~{}s^{-1}~{}Mpc^{-1}}$ is the
Hubble constant (Planck Collaboration et al., 2016) , and $f_{z}=2$ is a
correction factor for redshift evolution (Appendix A.6). We find that the
gamma-ray flux in our fiducial model is generally an order of magnitude or
more below the observed background intensity. However, model M2 can explain
the gamma-ray background intensity in the narrow $\sim 1$–$10~{}{\rm MeV}$
bands (dashed black line and brown points in Fig. 4a, Ackermann et al. 2015).
The origin of the background in this energy range has not been understood, as
other, previously proposed astrophysical contributions significantly
underpredict the level of the background (see thin solid, dashed, and dotted
brown lines in Fig. 4a). It is notable that our model can also explain the
neutrino background intensities as shown in Section
$\S~{}\ref{sec:UHECRs_Neutrinos}$ below.

###  3.3. Emission from NGC1068

We next consider whether the emission from internal shocks of the jets can
explain the gamma-ray emission from NGC1068. Note that NGC1068 is a type II
AGN, and its intrinsic X-ray luminosity is $\sim 7\times 10^{43}~{}{\rm
erg/s}$ (Marinucci et al., 2016), which is significantly brighter than the
X-ray emission by the jets from sBHs in the AGN disk (Fig. 2).

Fig. 4b compares the observed gamma-ray flux from NGC1068 and the predicted
fluxes in models M1–M3. Model M1 (solid line) significantly contributes to the
gamma-ray emission from NGC1068 between $\sim$100 MeV to $\sim$100 GeV because
model M1 avoids $\gamma\gamma$ annihilation owing to the high Lorentz factor.
On the other hand, models M2 and M3 have not been constrained by the current
gamma-ray observations. Also, the emission in infrared to X-ray bands is
significantly absorbed by dust (although this effect is not incorporated in
the predictions in Fig. 4 b) and the fluxes in the $\sim{\rm keV}$ – $\rm MeV$
bands receive significant contributions from coronae in the AGN. Thus, to test
models M2 and M3, MeV gamma rays are useful. MeV gamma-rays can be detected
with future gamma-ray telescopes, such as the Compton Spectrometer and Imager
(COSI) (Tomsick et al., 2019), the All-sky Medium Energy Gamma-ray Observatory
eXplorer (AMEGO)(Caputo et al., 2022), Gamma-Ray and AntiMatter Survey (GRAMS)
(Aramaki et al., 2020), eASTROGAM(De Angelis et al., 2021) and the Lunar
Occultation eXplorer (LOX)(Miller et al., 2019).

Figure 5.—  The contribution to the the background intensity (panel a) and the
flux from NGC1068 (panel b) for a single neutrino flavor by internal shocks
for models M1 (solid black), M2 (dashed black), and M3 (dotted black). The
observed neutrino intensity (Aartsen et al., 2020a) is presented by red points
(panel a), and the blue shaded region represents the 1, 2, and 3$~{}\sigma$
uncertainty on the spectrum measured by Aartsen et al. (2020b) (panel b).

##  4\. UHECRs and Neutrinos

In this section, we estimate whether high-energy protons and neutrinos
produced from the jets can explain the observed background fluxes and the
neutrino flux from NGC1068.

###  4.1. UHECRs

Since the adiabatic expansion is the most efficient cooling process for
protons accelerated in the jets in the models, the maximum proton energy
accelerated in internal shocks of the jets is given by the comparison between
the acceleration timescale in the Bohm limit (in which the particle mean free
path is assumed to be equal to the Larmor radius) and the expansion timescale
as

$\displaystyle E_{p,\rm max}=E^{\prime}_{p,\rm max}\Gamma_{j}=eB_{\rm
j}^{\prime}Z_{\rm diss}$ $\displaystyle\sim 10^{19}{\rm eV}\left(\frac{B_{\rm
j}^{\prime}}{4\times 10^{7}~{}{\rm G}}\right)\left(\frac{Z_{\rm diss}}{1\times
10^{9}~{}{\rm cm}}\right)$ (7)

where $B_{\rm j}^{\prime}$ is the magnetic field strength of the shocked jet,
$Z_{\rm diss}=2T_{\rm vari}c\Gamma_{\rm j}^{2}$ is the typical dissipation
radius for internal shocks, $e$ is the electric charge, and the primes denote
quantities in the fluid comoving frame. The maximum proton energies for models
M1, M2, and M3 are, respectively, $1\times 10^{17}$, $7\times 10^{17}$ and
$1\times 10^{19}~{}{\rm eV}$. It is notable that the CR production at a few
$10^{17}~{}{\rm eV}$ predicted in model M1 may be related to the observations
of CRs around this energy band (Buitink et al., 2016). On the other hand, to
produce high $E_{p,\rm max}$, low $\Gamma_{\rm j}$ and high $\epsilon_{\rm B}$
are required, and model M3 can explain the energy of the UHECRs.

The background intensity of UHECRs produced in the jets from sBHs in AGN disks
is estimated as

$\displaystyle E_{p}^{2}\Phi_{p}\sim\frac{\epsilon_{\rm CR}L_{\rm
j}}{R_{p}}n_{\rm rotBH,acc}f_{z}\frac{c}{4\pi H_{0}}$ $\displaystyle\sim
1\times 10^{-7}{\rm GeV~{}cm^{-2}~{}s^{-1}~{}sr^{-1}}$
$\displaystyle\times\left(\frac{\epsilon_{\rm
CR}}{0.15}\right)\left(\frac{L_{\rm j}}{10^{43}\,{\rm erg/s}}\right)$
$\displaystyle\times\left(\frac{n_{\rm rotBH,acc}}{5\times 10^{-4}\,{\rm
Mpc^{-3}}}\right)\left(\frac{R_{p}}{300}\right)^{-1}\left(\frac{f_{z}}{2}\right),$
(8)

where $E_{p}$ is proton energy, $R_{p}^{-1}=(E_{p}/E_{p,\rm
min})^{2-p}(p-2)/[1-(E_{p,\rm min}/E_{p,\rm max})^{p-2}]$, and $E_{p,\rm
min}=\Gamma_{\rm j}m_{p}c^{2}$ is the minimum proton energy. This value is
consistent with the observed value for the cosmic-ray intensity of $\sim
10^{-7}~{}{\rm GeV~{}cm^{-2}~{}sr^{-1}~{}s^{-1}}$ at $E_{p}\sim 10^{18}~{}{\rm
eV}$ (Blasi, 2013; Zweibel, 2013), while it is not obvious whether
extragalactic CRs of $E_{p}<10^{18}$ eV can propagate to the Earth within the
Hubble time (Kimura et al., 2015). Since the observed cosmic-ray background
intensity of $\lesssim 10^{-8}~{}{\rm GeV~{}cm^{-2}~{}sr^{-1}~{}s^{-1}}$ at
$E_{p}\sim 10^{20}~{}{\rm eV}$ (e.g. Zweibel, 2013) is lower by about an order
of magnitude compared to the production rate in models M1–M3, and the energy
can be predicted by model M3, the UHECRs could be explained by model M3. Also,
the high density of sBHs in AGN disks ($\sim 5\times 10^{-4}$–$5\times
10^{-3}~{}{\rm Mpc}^{-3}$) is consistent with the source density expected for
producing the UHECRs (Takami & Sato, 2009). On the other hand, since the heavy
metal abundance ratio of the UHECR compositions (Aab et al., 2014) is higher
than that of AGNs (Shields, 1996; Xu et al., 2018), we need to consider some
heavy-element enhancement process, such as re-acceleration of galactic CRs
(Caprioli, 2015; Kimura et al., 2018; Zirakashvili et al., 2022). Since the
cosmic-ray flux and abundance ratio are expected to be changed during
propagating to the Earth, detailed consistency of the flux and the metal
abundance distribution expected in this model are worth investigating in the
future.

###  4.2. Cosmic neutrino background

Interactions of high-energy protons lead to the production of pions, which
generate neutrinos via decay processes. We discuss the contribution of such
neutrino emission to the background neutrino flux and the neutrino flux from
NGC1068.

As cooling processes for baryons, we consider adiabatic cooling, $pp$
reaction, and $p\gamma$ reaction, where pions could be produced by inelastic
nuclear $pp$ reactions and photohadronic ($p\gamma$) reactions, respectively,
and the pions subsequently decay into neutrinos. Using the fraction of protons
producing pions through $pp$ and $p\gamma$ reactions ($f_{pp}$ and
$f_{p\gamma}$, Appendix A.4), the jet power, and the suppression factor due to
pion cooling $f_{\pi,\rm sup}$ (Kimura, 2022), the diffuse neutrino flux
produced from internal shocks of jets can be roughly estimated as (e.g.,
Razzaque et al. 2004; Murase et al. 2016)

$\displaystyle\epsilon_{\nu}^{2}\Phi_{\nu}\sim\frac{3K}{4(1+K)}\frac{{\rm
min}[1,{\max}(f_{p\gamma},f_{pp})]f_{\pi,\rm sup}\epsilon_{\rm CR}L_{\rm
j}}{R_{p}}$ $\displaystyle n_{\rm rotBH,acc}f_{z}\frac{c}{4\pi H_{0}}$
$\displaystyle\sim 1\times 10^{-13}{\rm GeV~{}cm^{-2}~{}s^{-1}~{}sr^{-1}}$
$\displaystyle\left(\frac{\epsilon_{\rm
CR}}{0.05}\right)\left(\frac{f_{pp}}{3\times
10^{-6}}\right)\left(\frac{f_{\pi,\rm sup}}{1}\right)\left(\frac{L_{\rm
j}}{2\times 10^{42}\,{\rm erg/s}}\right)$ $\displaystyle\left(\frac{n_{\rm
rotBH,acc}}{5\times 10^{-3}\,{\rm
Mpc^{-3}}}\right)\left(\frac{R_{p}}{500}\right)^{-1}\left(\frac{f_{z}}{2}\right)$
(9)

(solid black lines in Fig. 5), where $K=1$ and $K=2$, respectively, denote the
average ratio of charged to neutral pion for photohadronic ($p\gamma$) and
inelastic hadronuclear ($pp$) reaction, and $f_{\pi,\rm sup}$ is the fraction
of the energy loss before pions decay. Since the neutrino background intensity
at $\sim 10^{5}~{}{\rm GeV}$ is observed to be $\sim 10^{-7}$–$10^{-8}~{}{\rm
GeV~{}cm^{-2}~{}s^{-1}~{}sr^{-1}}$ (red points in Fig. 5a, Aartsen et al.
2015, 2020a), neutrinos produced by the internal shocks make only a minor
contribution to the background neutrino intensity in the fiducial model. On
the other hand, for low $\Gamma_{\rm j}$, $f_{pp}$ and $f_{p\gamma}$ are high,
and then the internal shocks could significantly contribute to the background
intensity (dashed and dotted black lines in Fig. 5a). As AGNs are expected to
be major production sites for neutrinos (Bartos et al., 2021) especially
around $E_{\nu}\sim 100~{}{\rm TeV}$ (Abbasi et al., 2022), the model with the
jets from sBHs in AGN disks may be a promising scenario for producing the
background neutrinos in the Universe.

###  4.3. Neutrino emission from NGC1068

We here discuss the possibility that the neutrinos from NGC1068 are produced
at internal shocks of jets launched from sBHs embedded in AGN disks. Fig. 5b
shows the neutrino flux from one jet in our models (models M2 and M3) and the
observed neutrino flux from NGC1068. The neutrino flux predicted by model M2
is not high enough to explain the observed flux (dashed black line). In the
model with high $\epsilon_{\rm B}$ (model M3), pion cooling via synchrotron
radiation is efficient, and therefore the neutrino flux at high energies in
$E_{\nu}\gtrsim 10^{4}~{}{\rm GeV}$ are suppressed. On the other hand, the
neutrino flux at lower energies in $E_{\nu}\lesssim 10^{4}~{}{\rm GeV}$, where
pion cooling is inefficient, is high, which is consistent with the observed
neutrino flux from NGC1068. Additionally, to reproduce the neutrino flux from
NGC1068, high $\epsilon_{\rm CR}$ and $L_{\rm j}$ are required. In the high
energy ranges around $\sim 10^{4}~{}{\rm GeV}$, where the atmospheric
backgrounds are much smaller, the observed flux is roughly consistent with the
prediction by model M3.

Although different values of $\epsilon_{\rm B}$ and $\epsilon_{\rm CR}$ (in
addition to $L_{\rm j}$) are required to explain the cosmic neutrino
background and the neutrino emission from NGC 1068, this is not unreasonable,
since we expect a broad distribution of $\epsilon_{\rm B}$ and $\epsilon_{\rm
CR}$. Multi-wavelength fits of gamma-ray burst afterglows revealed that the
distribution of $\epsilon_{\rm B}$ is indeed very broad (Panaitescu & Kumar,
2001; Santana et al., 2014) and $\epsilon_{\rm CR}$ is observationally less
constrained.111Theoretical models for magnetic amplification and cosmic-ray
acceleration in collisionless shocks are still incomplete with the current
computational resources (e.g. Caprioli & Spitkovsky, 2014; Caprioli et al.,
2020; Hu et al., 2022; Tomita et al., 2022), and we cannot predict the values
of $\epsilon_{\rm B}$ and $\epsilon_{\rm CR}$ in realistic astrophysical
environments from the first-principle calculations. If the peaks of the
$\epsilon_{\rm B}$ and $\epsilon_{\rm CR}$ distributions are as low as in
model M2 and the dispersions are large enough to contain a value as high as in
model M3, the neutrino emission from NC1068, and the background intensities of
gamma rays in $\sim{\rm MeV}$ bands, neutrinos in $\lesssim{\rm PeV}$ bands,
and UHECRs are all explained by a single unified model. This has not been
proposed so far, although we need to change the several parameters between
models M2 and M3. Also, the high value for $\epsilon_{\rm B}$, required to
explain the neutrino flux from NGC1068, may be difficult to be achieved. A
possible issue is that the jet produced by an sBH at $\sim 0.01~{}{\rm pc}$ in
NGC1068 needs to be directed to us. Since the SMBH mass and the accretion rate
of NGC1068 are high (Pier et al., 1994; Greenhill & Gwinn, 1997) and the
number of sBHs in an AGN disk is roughly proportional to the square root of
the SMBH mass and the accretion rate (Tagawa et al., 2020b, 2021), the number
of sBHs in an AGN disk in NGC1068 is higher by a factor of $\sim 14$ compared
to that in the fiducial model. Then, the probability that jets produced from
sBHs at $\sim 0.01~{}{\rm pc}$ in NGC1068 are directed to us is $\sim 0.3$,
which is viable. On the other hand, this also means that jets from $\sim 3$
sBHs at $\sim 1~{}{\rm pc}$ are directed to us, which may be problematic to be
consistent with the X-ray flux from NGC1068. This is because the X-ray
emission from corona in NGC1068 may be mostly absorbed at sub-parsec scales,
instead of several parsec scales, by considering the amount of mass estimated
in parsec scales (García-Burillo et al., 2016; Imanishi et al., 2018, 2020).
In that case, the X-ray emission from jets at $\sim 1~{}{\rm pc}$ may not be
fully absorbed and significantly exceed the observed X-ray flux. Hence,
dedicated estimates may be required to clarify whether or not this scenario is
consistent with hard X-ray observations by NuSTAR (Marinucci et al., 2016;
Zaino et al., 2020).

##  5\. Discussions

###  5.1. Fermi bubble

Two large gamma-ray bubbles, the so-called Fermi bubbles, have been discovered
above and below the center of our Galaxy by gamma-ray (Su et al., 2010), X-ray
(Bland-Hawthorn & Cohen, 2003), and radio telescopes (Finkbeiner, 2004). The
size of the bubble is $\sim 10~{}{\rm kpc}$, its shape is almost mirror-
symmetric with respect to the Galactic plane, and its origin is unknown (Yang
et al., 2018). We discuss whether the Fermi bubble may be related to sBHs in
an AGN disk. One of promising models for the origin of the Fermi bubble is the
leptonic jet model, in which a jet is estimated to be launched from the
central SMBH $1$–$3$ Myr ago with the duration of $\sim 0.1$–$0.5~{}{\rm Myr}$
(Yang et al., 2012; Guo & Mathews, 2012; Yang et al., 2022a), and the power of
the jet is roughly estimated to be $\sim 10^{43}$–$10^{44}~{}{\rm erg/s}$ (Guo
& Mathews, 2012; Yang et al., 2022a). We here propose that it is possible to
produce a similar bubble by jets launched from sBHs embedded in an AGN disk.
In our fiducial model, the power of the jet launched from an sBH is $L_{\rm
j}\sim 10^{42}$–$10^{43}~{}{\rm erg/s}$, and the active number of the jets is
$N_{\rm rotBH,AGN}\sim 10$ (Eq. 3.1). Thus, the total power of the jets
launched from the sBHs is consistent with the power required to explain the
properties of the Fermi bubble. Here we note that the symmetry of the Fermi
bubble is perhaps surprising, given that SMBH jets in general do not point
perpendicular to the galaxy disk (Kinney et al., 2000; Hopkins et al., 2012).
In the jets from the sBHs, the shocked gas is predicted to be symmetric with
respect to the galactic plane, if a pc-scale disk is aligned to the galactic
plane (Tagawa et al., 2020a). Since the angular momentum direction of the pc-
scale disks tend to be aligned to that of the nuclear star clusters due to
vector resonant relaxation (Kocsis & Tremaine, 2011; Impellizzeri et al.,
2019; Levin, 2022), the alignment of the pc-scale disk presumably occurs if
nuclear star clusters rotate in the same direction with the host galaxies,
which is often observed including our Galaxy (Levin & Beloborodov 2003; Yelda
et al. 2014; Do et al. 2020; Neumayer et al. 2020, but see also Kormendy & Ho
2013 for NGC4258 as a counter-example). Also, the duration of each jet from an
sBH is $\lesssim 10^{3}~{}{\rm yr}$ (Paper I), which may be consistent with
the uniform haze emission, as Yang et al. (2022a) suggested that multiple jets
may be required to explain the spatial distribution of the microwave emission
without being suppressed by magnetic pressure in the vicinity of the Galactic
center. Thus, jets from sBHs embedded in an AGN disk may be responsible for
producing such bubbles. Note that due to their intermittence and short
duration, the jets likely pass through high-pressure regions produced by
previous jets and could be significantly decelerated before interacting with
the inter-stellar medium. In this case, non-thermal emission from external
shocks of the jets might be inefficient. Such emission merits investigating in
the future.

###  5.2. Hadronic gamma-ray emission

When hadronuclear and photohadronic processes produce high-energy neutrinos in
astrophysical environments, these processes inevitably produce gamma rays
whose energy and luminosity are comparable to those of neutrinos. If these
gamma rays escape from the source, they are not absorbed during their
propagation to the Earth, and we should observe gamma rays of GeV–TeV
energies. However, the cosmic neutrino background intensity at $\sim 10$ TeV
is higher than the cosmic gamma-ray background intensity at $\sim 100$ GeV,
which implies that the cosmic neutrino source should be opaque to gamma rays
of $\gtrsim 100$ GeV (Murase et al., 2016; Kimura et al., 2021b). The same
arguments can be applied to the gamma-ray and neutrino fluxes from NGC1068,
and the neutrino emission region needs to be opaque to gamma rays of $\gtrsim
100$ MeV (Murase et al., 2020; Inoue et al., 2020).

In our models M2 and M3, gamma rays in $\gtrsim{\rm MeV}$ are significantly
suppressed by $\gamma\gamma$ annihilation, and thus, these models are
consistent with the gamma-ray and neutrino data for cosmic high-energy
backgrounds and NGC 1068. The typical energy of escaping gamma-rays should be
around MeV energies, and future MeV telescopes will be useful to probe the
hadron-induced electromagnetic cascade emission from these systems.

###  5.3. Caveats

In this section, we discuss caveats in our model. First, we fixed the values
of the parameters describing the AGN disks and the sBHs. However, in reality,
these parameter should vary from source to source, affecting their
contributions to the background. Additionally, we have not taken into account
the growth of sBHs due to gas accretion (Paper I, Section
$\S~{}\ref{sec:evolution}$) and mergers (Tagawa et al., 2021), which is a
promising pathway for sBH-sBH mergers reported by LIGO/Virgo/KAGRA (Abbott et
al., 2021). The evolution of the sBH mass influences the jet power ($L_{\rm
j}$) and $T_{\rm vari}$ (the dependence of the spectral energy distribution on
these quantities is presented in Appendix B). More precise estimates for the
background intensities considering these effects are desired to be conducted
in the future.

Meanwhile, in our model, the AGN disk’s role was considered to be just to feed
gas to sBHs – emission from shocks emerging from the AGN disk or from the
broad line regions is not considered unlike paper II, Tagawa et al. (2023b,
paper III), or Wang et al. (2021a). In this paper we considered persistent
high-energy emission, while in paper II and paper III, respectively, we
considered transient breakout emission from shocks emerging due to collision
between AGN disk gas and jets produced from merger remnants and solitary sBHs.
Note that we estimated that the transient emission from shocks emerging by
collision between the jets from sBHs and AGN disks cannot significantly
contribute to the gamma-ray and neutrino background intensities.

In model M1, a high Lorentz factor of 200 is assumed for fast shells (Appendix
A.2). On the other hand, it is unclear whether the shells in the jets can be
accelerated to such a high Lorentz factor via the Blandford-Znajek process
(e.g. Xiong & Zhang, 2014). Additionally, we assume a high contrast in the
Lorentz factors between fast and slow shells (e.g. 200 and 20 in model M1).
However, it is unclear whether such a high contrast is commonly realized (e.g.
Curd & Narayan, 2022). If either a high value or a high contrast in the
Lorentz factor is not achieved, the dissipation rate of the kinetic energy is
correspondingly overestimated in our model. Also, note that the equations in
Appendix A.2 and A.3 are approximated assuming that the Lorentz factor of the
shocked gas is much higher than 1. If the Lorentz factor of the shocked gas is
close to 1, the properties of electromagnetic emission and the jet structure
need to be appropriately modified. We note that the structures of the jets and
their dissipation are simply prescribed in this study. To quantitatively
discuss the dissipation of the jets and their emission, detailed numerical
simulations would be required.

The neutrino spectrum in our model might be affected by a few processes that
is not included in our calculations. First, we have not considered suppression
of neutrino production due to the Bethe-Heitler process, which is determined
by photon spectra. The photon spectra is also modified if we take into account
the electromagnetic cascade emission initiated by two-photon interactions and
hadronic processes. In models M2 and M3, such emission likely produces flatter
spectral energy distribution below the energy limited by $\gamma\gamma$
annihilation. We confirmed by performing numerical computations222We modify
the codes used in Kimura et al. (2019); Kimura & Toma (2020) to match the
physical condition in the current model.  that the suppression is not
significant in our models because of the flat photon spectral produced by the
cascade emission. The muon cooling is also not taken into account, which may
reduce a neutrino flux by a factor of 2. In model M3, this process is
important at $E_{\nu}\sim 1-10~{}{\rm TeV}$.

###  5.4. Evolution of sBHs

Here, we discuss how much sBHs grow by accretion and their possible fate (see
also Levin 2007; McKernan et al. 2012), as their growth is almost inevitable
if they are to explain the large power required to produce the various
background intensities discussed in this paper. In the fiducial model, the jet
power is set to $2\times 10^{42}~{}{\rm erg/s}$, corresponding to the
accretion rate of $3\times 10^{-4}~{}{\rm M}_{\odot}/{\rm yr}(\eta_{\rm
j}/0.1)^{-1}$. Assuming that each sBH accretes for $\sim 4~{}{\rm Myr}$ per
AGN phase and the accretion rate is constant, the mass of an sBH increases to
$m\sim 10^{3}~{}{\rm M}_{\odot}$. The timescale for these (modestly)
intermediate-mass BHs to migrate and inspiral to the central SMBH only by
gaseous and GW torques is longer than the active AGN phase (Tagawa et al.,
2022), but we consider migration by stellar torques during longer quiescent
phases. Since the relaxation timescale of the nuclear star cluster is $t_{\rm
rel}\lesssim 10^{11}~{}{\rm yr}$ (Merritt, 2010; Kocsis & Tremaine, 2011), and
the migration timescale of an sBH with the mass $m$ is roughly given by $\sim
t_{\rm rel}m_{\rm ave}/m\sim 10^{8}~{}{\rm yr}(t_{\rm rel}/10^{11}~{}{\rm
yr})[(m/m_{\rm ave})/10^{3}]^{-1}$, where $m_{\rm ave}$ is the average mass of
stars in a nuclear star cluster, the grown BHs migrate towards and accrete
onto the central SMBH during quiescent phases, whose duration is roughly given
as $\sim t_{\rm H}/N_{\rm AGN,H}\sim 10^{9}~{}{\rm yr}(t_{\rm
H}/10^{10}~{}{\rm yr})(N_{\rm AGN,H}/10)^{-1}$, where $N_{\rm AGN,H}$ is the
number of AGN phases per galaxy during the Hubble timescale, $t_{\rm H}$.
During the Salpeter timescale of $\sim 40~{}{\rm Myr}$, inspirals of a hundred
sBHs with the mass $10^{3}~{}{\rm M}_{\odot}$ enhances the SMBH mass
($10^{6}~{}{\rm M}_{\odot}$ in the fiducial model) by $\sim 10\%$. Thus, if
the high-energy phenomena are caused by sBHs in AGN disks and the number of
sBHs in an AGN disk is not significantly underestimated in our model, sBHs
grow to intermediate-mass BHs, and then, they merge with SMBHs during
quiescent phases, and can be observed as intermediate-mass ratio inspirals
(IMRIs) by the Laser Interferometer Space Antenna (Amaro-Seoane et al., 2022).

##  6\. Conclusions

In this paper, we considered high-energy EM, neutrino, and cosmic-ray
emissions arising from BZ jets launched from rapidly accreting and spinning
sBHs embedded in AGN disks. Our main results are summarized as follows:

  1. 1.

GeV gamma-ray emission from nearby radio-quiet Seyfert galaxies can be
explained by emission from the jets launched from sBHs if the Lorentz factor
of the jets is high ($\Gamma_{\rm j}\gtrsim 30$). This model can contribute
$\gtrsim 50\%$ of the gamma-ray emission from NGC1068 and $\sim 10-20$% of the
background gamma-ray intensity in the $\sim{\rm MeV}$–$\rm GeV$ bands.

  2. 2.

For low $\Gamma_{\rm j}$ ($\lesssim 4$), neutrino emission from the jets can
explain the cosmic neutrino background intensity at neutrino energies
$\lesssim 10^{6}~{}{\rm GeV}$. The background in the $\sim{\rm MeV}$ bands can
also be reproduced in this model if the efficiency of magnetic field
amplification ($\epsilon_{\rm B}$) is low. In this case, we predict that
future gamma-ray telescopes in $\sim{\rm MeV}$ bands can detect MeV gamma-rays
from nearby radio-quiet AGNs, providing a test of our model.

  3. 3.

Neutrino emission from NGC1068 can be explained by the jet model with moderate
$\Gamma_{\rm j}$ and high $\epsilon_{\rm B}$, $\epsilon_{\rm CR}$, and $L_{\rm
j}$. With this parameter set, jets by sBHs in AGNs can accelerate protons up
to energies of $\sim 10^{19}~{}{\rm eV}$ and account for the observed
intensity of UHECRs at these energies.

  4. 4.

If $\epsilon_{\rm B}$ and $\epsilon_{\rm CR}$ have a broad distribution from
source to source, our model can simultaneously explain the neutrino flux from
NGC1068 as well as the background intensities of gamma rays in $\sim{\rm MeV}$
bands, neutrinos in $\lesssim{\rm PeV}$ bands, and UHECRs, for the first time.

While AGN are known to have relativistic jets driven by their central SMBH,
our results in this paper suggest that the population of stellar-mass BHs,
embedded in the accretion disk fueling the central SMBH, can collectively
produce similar phenomena with energetics similar to that from the SMBH due to
the larger number of jets from sBHs. We found that the jets from sBHs in AGN
disks can explain various phenomena, which are observed to be unrelated to
jets from SMBHs. Our scenario is consistent with the observations that the
background neutrinos (Murase & Waxman, 2016) and UHECR (Takami et al., 2016)
emission are produced by the sources with a high local density.

We thank Kohta Murase for useful comments. This work was financially supported
by Japan Society for the Promotion of Science (JSPS) KAKENHI grant Number
JP21J00794 (HT) and 22K14028 (SSK). S.S.K. acknowledges the support by the
Tohoku Initiative for Fostering Global Researchers for Interdisciplinary
Sciences (TI-FRIS) of MEXT’s Strategic Professional Development Program for
Young Researchers. Z.H. was supported by NASA grant NNX15AB19G and NSF grants
AST-2006176 and AST-1715661.

##  Appendix A Emissions

We here describe the properties of the internal shocks in the jets and non-
thermal emission produced from them.

###  A.1. Accretion onto sBHs

We outline how the jet is launched from an accreting sBH embedded in an AGN
disk (see Paper I, for details). In the AGN disk, the gas accreting onto sBH
forms a circum-sBH disk (CsBD). When the CsBD is advection dominated, a
magnetically dominated state can be realized owing to the accumulation of the
magnetic flux in the vicinity of the sBH (e.g. Meier, 2001; Cao, 2011; Kimura
et al., 2021c). Even if the magnetic flux is initially weak, the outflow from
the disk converts the toroidal magnetic field generated by the shear motion
into a poloidal field (Liska et al., 2020). Such advection-dominated flows are
expected for super-Eddington accretion rates (Abramowicz et al., 1988) or low
accretion rates (e.g. Narayan & Yi, 1994; Blandford & Begelman, 1999). In
these cases, the jets from spinning sBHs can be launched through the BZ
process (Blandford & Znajek, 1977; Curd & Narayan, 2022).

Because super-Eddington accretion is predicted in an inner CsBD (Paper I), the
BZ jet is expected to be launched from rapidly accreting and spinning sBHs in
AGN disks. In the process, the jet power ($L_{\rm j}$) is proportional to the
mass accretion rate onto the sBH (${\dot{M}}_{\rm sBH}$) as

$\displaystyle L_{\rm j}=\eta_{\rm j}{\dot{M}}_{\rm sBH}c^{2},$ (A1)

where $\eta_{\rm j}$ is the jet conversion efficiency to the kinetic energy,
which is approximated as $\eta_{\rm j}\sim a_{\rm sBH}^{2}$ for a magnetically
dominated state (e.g. Tchekhovskoy et al., 2010; Narayan et al., 2022),
$a_{\rm sBH}$ is the dimensionless spin of the sBH, and $c$ is the speed of
light. We assume that the accretion rate onto sBHs in the AGN disk is given by
the Bondi-Hoyle-Lyttoleton (BHL) rate as used in Eq. (1) of Paper I. In the
formula, the accretion rate is determined by the AGN disk density and
temperature, and the Hill radius of the sBH.

###  A.2. Properties of jets

To dissipate the kinetic energy of jets, we assume multiple shells with
different Lorents factors collide with each other at some dissipation radius
($Z_{\rm diss}$), which is widely assumed to explain the prompt emission in
gamma-ray bursts. The relative Lorents factor between the slower (with the
Lorents factor $\Gamma_{\rm s}=(1-\beta_{\rm s}^{2})^{-1/2}$) and faster
shells (with the Lorents factor $\Gamma_{\rm r}=(1-\beta_{\rm r}^{2})^{-1/2}$)
is

$\displaystyle\Gamma_{\rm rel}=\Gamma_{\rm s}\Gamma_{\rm r}(1-\beta_{\rm
s}\beta_{\rm r})$ (A2)

The relative Lorents factor $\Gamma_{\rm rel}$ is also related to the relative
Lorents factors between the unshocked faster and slower shells in the rest
frame of the shocked fluid ($\Gamma_{12}$ and $\Gamma_{34}$) as

$\displaystyle\Gamma_{\rm rel}=\Gamma_{\rm 12}\Gamma_{\rm
34}+\sqrt{\Gamma_{\rm 12}^{2}-1}\sqrt{\Gamma_{\rm 34}^{2}-1},$ (A3)

and the ratio of the number densities of the faster ($n_{\rm r}^{\prime}$) and
slower ($n_{\rm s}^{\prime}$) shells in the fluid rest frame is related to
them as

$\displaystyle\frac{n_{\rm r}^{\prime}}{n_{\rm s}^{\prime}}=\frac{(\Gamma_{\rm
34}-1)(4\Gamma_{\rm 34}+3)}{(\Gamma_{\rm 12}-1)(4\Gamma_{\rm 12}+3)}$ (A4)

The Lorents factor of the shocked jet in the rest frame of the sBH is
approximately given as

$\displaystyle\Gamma_{\rm j}\approx\Gamma_{\rm r}(\Gamma_{\rm
12}-\sqrt{\Gamma_{\rm 12}^{2}-1}).$ (A5)

If we assume $\Gamma_{\rm r}=200$, $\Gamma_{\rm s}=20$, and $\frac{n_{\rm
r}^{\prime}}{n_{\rm s}^{\prime}}=\Gamma_{\rm s}^{2}/\Gamma_{\rm r}^{2}$, then
$\Gamma_{12}\sim 3.6$, and $\Gamma_{\rm j}\sim 28$, which we assume as
fiducial values. Since the results are sensitive to $\Gamma_{\rm j}$, we also
consider the models (M2, M3) with $\Gamma_{\rm r}=30$ and $\Gamma_{\rm s}=3$,
in which $\Gamma_{12}\sim 3.7$, and $\Gamma_{\rm j}=4.2$.

###  A.3. Non-thermal photons

We assume that the fraction $\epsilon_{e}$ of the kinetic energy of the shock
is used to accelerate electrons in a collisionless shock. Here, the plasma
and/or MHD instabilities are considered to amplify the magnetic field to
$\epsilon_{\rm B}\lesssim 10^{-3}$–$10^{-1}$ and electrons are accelerated via
the first-order Fermi process with the energy fraction of
$\epsilon_{e}\lesssim 10^{-2}$–$0.3$ (from observations, e.g., Waxman & Loeb
1999; Panaitescu & Kumar 2001; Frail et al. 2005; Uchiyama et al. 2007;
Santana et al. 2014, and theoretical studies, e.g., Medvedev & Loeb 1999;
Chang et al. 2008; Spitkovsky 2008; Martins et al. 2009; Keshet et al. 2009;
Sironi et al. 2013; Tomita et al. 2019). Assuming the fast cooling regime (Fan
& Piran, 2008), which is adequate for the internal shocks, the synchrotron
luminosity is calculated as

$\displaystyle L_{\rm syn}\sim\frac{L_{\rm
kin}\epsilon_{e}f_{\gamma\gamma}}{1+Y_{\rm SSC}},$ (A6)

(e.g. Fan & Piran, 2008), where $f_{\gamma\gamma}$ is the attenuation fraction
by the $\gamma\gamma$ annihilation given below, $Y_{\rm SSC}$ is the powers of
synchrotron self-Compton scattering compared to that of synchrotron emission,
and is calculated as $Y_{\rm SSC}=(-1+\sqrt{1+4\epsilon_{e}/\epsilon_{\rm
B}})/2$ (Fan & Piran, 2008).

We assume that electrons are accelerated in the shock to a power-law
distribution of Lorentz factor $\gamma^{\prime}_{e}$ as
$N(\gamma_{e}^{\prime})d\gamma_{e}^{\prime}\propto{\gamma_{e}^{\prime}}^{-p}d\gamma_{e}^{\prime}$
with a minimum ($\gamma_{\rm m}^{\prime}$) and maximum Lorentz factors
($\gamma_{\rm max}^{\prime}$). The minimum Lorentz factor $\gamma_{\rm
m}^{\prime}$ is

$\displaystyle\gamma_{\rm
m}^{\prime}\sim\epsilon_{e}\left(\frac{p-2}{p-1}\right)\frac{m_{p}}{m_{e}}(\Gamma_{\rm
12}-1)$ $\displaystyle\sim
30~{}\left(\frac{\epsilon_{e}}{0.02}\right)\left(\frac{\Gamma_{\rm
12}-1}{2.6}\right)$ (A7)

for $p=2.5$, where $m_{e}$ is the electron mass.

By comparing the cooling by the synchrotron radiation and the acceleration by
the first-order Fermi acceleration mechanism, the maximum Lorentz factor of
electrons is

$\displaystyle\gamma_{\rm max}^{\prime}=\left(\frac{6\pi e}{\sigma_{\rm
T}B^{\prime}_{\rm j}\xi}\right)^{1/2}\sim 1\times
10^{6}~{}\xi^{-1/2}\left(\frac{\Gamma_{\rm j}-1}{27}\right)^{-1/4}$
$\displaystyle\left(\frac{\epsilon_{\rm
B}}{0.01}\right)^{-1/4}\left(\frac{n^{\prime}_{p}}{9\times 10^{10}\,{\rm
cm^{-3}}}\right)^{-1/4}$ (A8)

where $\xi$ is the parameter representing the ratio of the mean free path to
the Larmor radius of electrons, which is adopted to be $\xi=1$ in this paper,

$\displaystyle B^{\prime}_{\rm j}=(8\pi\epsilon_{\rm B}e^{\prime}_{\rm
j})^{1/2}\sim 1\times 10^{4}~{}{\rm G}~{}\left(\frac{\Gamma_{\rm
12}-1}{2.6}\right)^{1/2}\left(\frac{\epsilon_{\rm B}}{0.01}\right)^{1/2}$
$\displaystyle\left(\frac{n^{\prime}_{p}}{9\times 10^{10}\,{\rm
cm^{-3}}}\right)^{1/2},$ (A9)

is the magnetic field,

$\displaystyle e^{\prime}_{\rm j}=(\Gamma_{\rm 12}-1)n^{\prime}_{p}m_{p}c^{2}$
(A10)

is the internal energy density of the shocked jet, and $n^{\prime}_{p}$ is the
proton number density of the shocked jet.

In the fiducial model, non-thermal emission is characterized by fast cooling
regimes (e.g. Sari et al., 1998; Fan & Piran, 2008). The cooling timescale for
electrons with $\gamma_{\rm{m}}^{\prime}$ is

$\displaystyle t_{\rm c}(\gamma_{\rm{m}}^{\prime})\sim
0.009\,{\rm{s}}~{}\left(\frac{\epsilon_{\rm{B}}}{0.01}\right)^{-1}\left(\frac{\gamma_{\rm
12}-1}{2.6}\right)^{-1}$ $\displaystyle\left(\frac{n^{\prime}_{p}}{9\times
10^{10}\,{\rm
cm^{-3}}}\right)^{-1}\left(\frac{\gamma_{\rm{m}}^{\prime}}{30}\right)^{-1}\left(\frac{\Gamma_{\rm
j}}{30}\right)^{-1}.$ (A11)

From $t_{\rm c}(\gamma_{\rm{m}}^{\prime})$, the typical shell width of
electrons with $\gamma_{\rm{m}}^{\prime}$ emitting the synchrotron photons is
approximated as $\Delta_{\rm shell}(\gamma_{\rm{m}}^{\prime})\sim t_{\rm
c}(\gamma_{\rm{m}}^{\prime})c\sim 3\times 10^{8}\,{\rm cm}~{}[t_{\rm
c}(\gamma_{\rm{m}}^{\prime})/0.009~{}{\rm{s}}]$.

The Lorentz factor at which self-absorption becomes effective is

$\displaystyle\gamma_{\rm a}^{\prime}=\gamma_{\rm
m}^{\prime}\times\left(\tau_{\rm q}C_{q+1}\right)^{1/(q+4)}$ (A12)

(Rybicki & Lightman, 1979; Fouka & Ouichaoui, 2009, 2011), where

$\displaystyle\tau_{q}=\frac{\pi}{3\sqrt{2}}\frac{(q^{2}+q-2)\gamma_{\rm
1}^{\prime-5}}{1-(\gamma_{\rm 2}^{\prime}/\gamma_{\rm
1}^{\prime})^{-q+1}}\frac{en^{\prime}_{\rm p}\Delta_{\rm
shell}(\gamma_{\rm{a}}^{\prime})}{\Gamma_{\rm j}B^{\prime}_{\rm j}},$ (A13)

$\gamma_{2}$ and $\gamma_{1}$ are the maximum and minimum Lorentz factors of
non-thermal electrons with the power-law index of $q$, respectively,

$\displaystyle
C_{q}=\frac{2^{(q+1)/2}}{q+1}\Gamma\left(\frac{q}{4}-\frac{1}{12}\right)\Gamma\left(\frac{q}{4}+\frac{19}{12}\right),$
(A14)

$\Gamma$ is the Gamma function, and we assume that electrons are randomly
oriented in the frame of the shocked fluid. We set
$\gamma_{2}^{\prime}=\gamma_{\rm max}^{\prime}$,
$\gamma_{1}^{\prime}=\gamma_{\rm m}^{\prime}$, and $q=p_{e}$ for $\gamma_{\rm
m}^{\prime}<\gamma_{\rm a}^{\prime}$, and $\gamma_{2}^{\prime}=\gamma_{\rm
m}^{\prime}$, $\gamma_{1}^{\prime}=\gamma_{\rm a}^{\prime}$, and $q=2$ for
$\gamma_{\rm a}^{\prime}<\gamma_{\rm m}^{\prime}$.

Using the variables derived above, we assume the spectra shapes for
synchrotron emission and synchrotron-self Compton scattering as prescribed in
Paper II, while we additionally take into account the bolometric correction as
prescribed by $R_{p}$ in Eq. (4.1). The spectra shapes for models M1–M3 are
presented in Fig. 2. In the fiducial model, we do not consider the second
order inverse Compton scattering as the Klein-Nishina effect is effective.

When the optical depth to the $\gamma\gamma$ annihilation
($\tau_{\gamma\gamma}$) exceeds $\gtrsim 1$, the high energy photons with
$h\nu^{\prime}>m_{e}c^{2}=511~{}{\rm keV}$ can be absorbed, where $h$ is the
Plank constant, and $\nu$ is the frequency of photon. To take into account
$\gamma\gamma$ attenuation, we refer to Eq. (6) of Kimura (2022) for
$\tau_{\gamma\gamma}$. We assume that the attenuation factor by the
$\gamma\gamma$ annihilation is given by $f_{\gamma\gamma}(\nu)={\rm
exp}[1-\tau_{\gamma\gamma}(h\nu)]/\tau_{\gamma\gamma}(h\nu)$. The peak
luminosity for attenuating photons is assumed to be Eq. (A6), and the peak
frequency to be the maximum of the minimum and absorption frequencies of the
synchrotron radiation.

###  A.4. Pion production efficiency

We calculate the efficiency of the $p\gamma$ reaction ($f_{p\gamma}$) and the
pion cooling suppression factor ($f_{\pi,sup}$) referring to Eq. (28) and two
asymptotes stated below Eq. (33) of Kimura (2022), respectively. Comparing the
dynamical timescale ($t^{\prime}_{\rm dyn}$) to the $pp$ cooling timescale
($t^{\prime}_{pp}$), the efficiency of $pp$ interaction is calculated as

$\displaystyle f_{pp}\approx n^{\prime}_{p}\kappa_{p}\sigma_{pp}Z_{\rm
diss}/\Gamma_{j}$ $\displaystyle\sim 3\times 10^{-6}\left(\frac{L_{\rm
j}/p_{\theta}}{1\times 10^{44}~{}{\rm erg}}\right)$
$\displaystyle\left(\frac{\Gamma_{j}}{30}\right)^{-3}\left(\frac{Z_{\rm
diss}}{5\times 10^{10}~{}{\rm cm}}\right)^{-1}$ (A15)

(e.g. Murase et al., 2014), where

$\displaystyle n^{\prime}_{p}=L_{\rm j}/4\pi\Gamma_{j}^{2}p_{\theta}Z_{\rm
diss}^{2}m_{p}c^{3}$ (A16)

is the proton number density in the jet, $\kappa_{p}\approx 0.5$ is the proton
inelasticity, and $\sigma_{pp}\approx 4\times 10^{-26}~{}{\rm cm}^{2}$ is the
cross section of the $pp$ interactions at $\sim 10$–$100~{}{\rm TeV}$.

To derive the efficiency of $p\gamma$ interaction, we adopt the delta function
approximation. We assume that most interactions occur via $\Delta$-resonance
process at the photon energy around the resonance peak (${\bar{\epsilon}}_{\rm
pk}\simeq 0.3~{}{\rm GeV}$), and approximate the cross section and
inelasticity to be

$\displaystyle\sigma_{p\gamma}\kappa_{p\gamma}\simeq\sigma_{\Delta}\kappa_{\Delta}\Delta{\bar{\epsilon}}_{\rm
pk}\delta({\bar{\epsilon}}_{\gamma}-{\bar{\epsilon}}_{\rm pk})$ (A17)

where $\sigma_{\Delta}\sim 5\times 10^{-28}~{}{\rm cm^{-2}}$,
$\kappa_{\Delta}\simeq 0.2$, and ${\bar{\epsilon}}_{\rm pk}\simeq 0.3~{}{\rm
GeV}$ are the cross section, inelasticity, and the photon energy at the
resonance peak, $\Delta{\bar{\epsilon}}_{\rm pk}\sim 0.2~{}{\rm GeV}$ is the
peak width, $\delta(x)$ is the Dirac delta function, and
${\bar{\epsilon}}_{\gamma}$ is the photon energy in the proton rest frame. For
$E_{p}<E_{p,\rm br}$, where

$\displaystyle E_{p,\rm br}=\frac{\Gamma_{j}^{2}{\bar{\epsilon}}_{\rm
pk}m_{p}c^{2}}{2E_{\gamma,\rm br}},$ (A18)

$E_{\gamma,\rm br}$ is the peak photon energy in the BH rest frame, the
fraction of cosmic-ray protons producing pions through the photomeson
production is approximated as

$\displaystyle
f_{p\gamma}\approx\frac{\sigma_{\Delta}\kappa_{\Delta}\Delta{\bar{\epsilon}}_{\rm
pk}L_{\gamma,\rm iso,br}}{2\pi c{\bar{\epsilon}}_{\rm pk}\Gamma_{j}^{2}Z_{\rm
diss}E_{\gamma,\rm br}(3-s)}(E_{p}/E_{p,\rm br})^{1-s}$ (A19)

(e.g. Kimura, 2022), where $L_{\gamma,\rm iso,br}$ is the isotropic
synchrotron luminosity at $E_{\gamma}=E_{\gamma,\rm br}$, and $s$ is the
power-law slope for the synchrotron spectrum for $E_{\gamma}>E_{\gamma,\rm
br}$. In model M2,

$\displaystyle f_{p\gamma}\sim 7\left(\frac{L_{\gamma,\rm iso,br}}{3\times
10^{41}\,{\rm erg/s}}\right)\left(\frac{Z_{\rm diss}}{10^{9}\,{\rm
cm}}\right)^{-1}$ $\displaystyle\left(\frac{\Gamma_{\rm
j}}{4}\right)^{-3}\left(\frac{E_{\gamma,\rm br}}{200~{}{\rm eV}}\right)^{-1}$
(A20)

at $E_{p}=E_{p,\rm br}$ and $E_{p,\rm br}\sim 3\times 10^{6}~{}{\rm
GeV}(\Gamma_{j}/4)^{2}(E_{\gamma,\rm br}/200~{}{\rm eV})^{-1}$. Since the
delta function approximation cannot predict $f_{p\gamma}$ above $E_{p,\rm
br}$, we do not present the neutrino luminosities in these regions in Fig. 5.
This approximation does not affect the detectability of neutrinos as the pion
cooling suppresses the neutrino fluence in the high-energy ranges. In models
M1–M3, the $pp$ reaction dominates the neutrino reaction in lower energy
ranges, while the ${p\gamma}$ reaction dominates over the $pp$ reaction in
middle energy ranges of $10^{4}~{}{\rm GeV}\lesssim E_{\nu}\lesssim
10^{6}~{}{\rm GeV}$ (Fig. 5).

Figure 6.—  The parameter range for efficient particle acceleration on
$\Gamma_{\rm j}$–$L_{\rm j,iso}$ plane for $\Gamma_{\rm 12}=4$ and $T_{\rm
vari}=10^{-3}~{}{\rm s}$. The parameters adopted for models M1, M2, and M3 are
presented by orange, cyan, and red circles, respectively.

###  A.5. Parameter space for particle acceleration

Here we discuss the parameter space in which electrons and protons can be
efficiently accelerated in internal shocks of the jets. If the upstream gas is
optically thick to incident photons, the shocks are mediated by interactions
with the photons (e.g. Ito et al., 2020). In this case, the velocity jump at
the shock is gradual, and particle acceleration should be inefficient. The
optical depth of the upstream shell is estimated as

$\displaystyle\tau_{\rm u}=n^{\prime}_{\rm r}\sigma_{\rm T}\frac{Z_{\rm
diss}}{\Gamma_{\rm r}}$ $\displaystyle=\frac{L_{\rm j,iso}\sigma_{\rm
T}}{4\pi\Gamma_{\rm r}^{3}Z_{\rm diss}m_{p}c^{3}}$ $\displaystyle\simeq
0.004\left(\frac{L_{\rm j,iso}}{10^{44}~{}{\rm
erg/s}}\right)\left(\frac{\Gamma_{\rm
12}}{4}\right)^{-3}\left(\frac{\Gamma_{\rm
j}}{4}\right)^{-5}\left(\frac{T_{\rm vari}}{10^{-3}~{}{\rm s}}\right)^{-1}$
(A21)

where $\sigma_{\rm T}$ is the Thomson scattering cross section, and we used
the approximation of $\Gamma_{\rm r}\sim 2\Gamma_{12}\Gamma_{\rm j}$ in the
last transformation. Hence, in the models we investigated (M1–M3), internal
shocks are collisionless and particles can be efficiently accelerated (see
Fig. 6).

###  A.6. Correction for redshift evolution

We here discuss the correction factor ($f_{z}$) due to the redshift evolution
of the number density of sources (Waxman & Bahcall, 1998), which needs to be
taken into account to estimate the diffuse background intensity. For instance,
$f_{z}\sim 3$ and $f_{z}\sim 0.6$, when the number density of sources is
proportional to $(z+1)^{3}$ and $(z+1)^{0}$, respectively. The value of
$f_{z}$ in our scenario should depend on two things: the number of sBHs
captured by AGN disks and the typical accretion rate onto the sBHs. The
accretion rate onto sBHs at a fixed radius is higher for highly-accreting
AGNs, but sBHs around heavier SMBHs tend to be captured in outer regions,
where the accretion rate onto sBHs is lower. These two effects likely
compensate each other, and thus, we assume that the accretion rate onto sBHs
is independent of the AGN luminosity. On the other hand, the number of sBHs is
higher for high-luminosity AGNs ($L_{X}\sim 10^{43}~{}{\rm erg/s}$) (Tagawa et
al., 2021), whose number density strongly depends on redshift,
$\sim(1+z)^{3}-(1+z)^{4}$ (Ueda et al., 2014). If we instead assume the number
of sBHs is independent of the AGN luminosity, sBHs in lower-luminosity AGNs
should also provide a significant contribution, whose number density is almost
independent of redshift, $\sim(1+z)^{0}$ (Ueda et al., 2014). Thus, the source
evolution in our scenario should be between $(1+z)^{3}$ and $(1+z)^{0}$. In
this study, we adopt $f_{z}=2$.

Figure 7.—  Same as Fig. 4, but with parameters varied from the fiducial model
(model M1) as indicated in the legend.  Figure 8.—  Same as Fig. 5, but with
parameters varied from model M2 as indicated in the legend.

##  Appendix B Parameter dependence

In this section, we present the parameter dependence of the contribution of
emission from internal shocks in jets to the gamma-ray (Fig. 7) and neutrino
background intensities (Fig. 8).

Fig. 7 shows the dependence of the gamma-ray flux from internal shocks on the
parameters, $p$, $\epsilon_{e}$, $T_{\rm vari}$, and $L_{\rm j}$. All of these
parameters influence the minimum synchrotron frequency ($\nu_{\rm
m}=\gamma_{\rm m}^{2}\nu_{\rm syn}$, where $\nu_{\rm syn}$ is the synchrotron
frequency). This is obvious that $\gamma_{\rm m}$ depends on $\epsilon_{e}$
and $p$ (Eq. A.3), while $\nu_{\rm syn}$ depends on $T_{\rm vari}$ and $L_{\rm
j}$ as $\nu_{\rm syn}\propto B_{\rm j}^{\prime}\propto n_{p}^{\prime
1/2}\propto L_{\rm j}^{1/2}T_{\rm vari}^{-1}$. In addition, $p$ influences the
power-law slope above $\nu_{\rm m}$ (dashed black line in Fig. 7a), and
$\epsilon_{e}$ and $L_{\rm j}$ influence the total radiation energy by non-
thermal emission (dotted black lines in Fig. 7). Note that $n_{\rm rotBH,acc}$
also linearly affects the intensity ($E_{\gamma}^{2}\Phi_{\gamma}\propto
n_{\rm rotBH,acc}$). Thus, the gamma-ray intensity depends on various
parameters. Nevertheless, gamma-ray emission is not eliminated by these
parameters, unlike the model with different $\Gamma_{\rm j}$ as shown in the
main manuscript (Fig. 4).

Fig. 8 shows the dependence of the neutrino flux from the internal shocks on
several parameters. The neutrino intensity in low energy ranges is low and
high for the long-$T_{\rm vari}$ or high-$L_{\rm j}$ models (dotted black and
solid orange lines) since $f_{pp}$ is proportional to the high jet density
($f_{pp}\propto n_{p}^{\prime}\propto L_{\rm j}T_{\rm vari}^{-2}$). Also, for
long $T_{\rm vari}$ or low $L_{\rm j}$, cooling by synchrotron radiation
during pion phases is inefficient, due to weak magnetic fields in shocks, and
then the neutrino intensities are high in high energy ranges. Similar to the
gamma-ray intensity, $n_{\rm rotBH,acc}$ linearly influences the neutrino
intensity ($E_{\nu}^{2}\Phi_{\nu}\propto n_{\rm rotBH,acc}$).

## References

  * Aab et al. (2014) Aab, A., Abreu, P., Aglietta, M., et al. 2014, Phys. Rev. D, 90, 122005, doi: 10.1103/PhysRevD.90.122005
  * Aartsen et al. (2013) Aartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013, Phys. Rev. Lett., 111, 021103, doi: 10.1103/PhysRevLett.111.021103
  * Aartsen et al. (2015) Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2015, ApJ, 809, 98, doi: 10.1088/0004-637X/809/1/98
  * Aartsen et al. (2020a) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020a, Phys. Rev. Lett., 125, 121104, doi: 10.1103/PhysRevLett.125.121104
  * Aartsen et al. (2020b) —. 2020b, Phys. Rev. Lett., 124, 051103, doi: 10.1103/PhysRevLett.124.051103
  * Abbasi et al. (2022) Abbasi, R., Ackermann, M., Adams, J., et al. 2022, Phys. Rev. D, 106, 022005, doi: 10.1103/PhysRevD.106.022005
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, Physical Review X, 11, 021053, doi: 10.1103/PhysRevX.11.021053
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
  * Abramowicz et al. (1988) Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 646, doi: 10.1086/166683
  * Acciari et al. (2019) Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2019, ApJ, 883, 135, doi: 10.3847/1538-4357/ab3a51
  * Ackermann et al. (2015) Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 799, 86, doi: 10.1088/0004-637X/799/1/86
  * Ajello et al. (2020) Ajello, M., Di Mauro, M., Paliya, V. S., & Garrappa, S. 2020, ApJ, 894, 88, doi: 10.3847/1538-4357/ab86a6
  * Ajello et al. (2008) Ajello, M., Greiner, J., Sato, G., et al. 2008, ApJ, 689, 666, doi: 10.1086/592595
  * Ajello et al. (2017) Ajello, M., Atwood, W. B., Baldini, L., et al. 2017, ApJS, 232, 18, doi: 10.3847/1538-4365/aa8221
  * Alves Batista et al. (2019) Alves Batista, R., Biteau, J., Bustamante, M., et al. 2019, Frontiers in Astronomy and Space Sciences, 6, 23, doi: 10.3389/fspas.2019.00023
  * Amaro-Seoane et al. (2022) Amaro-Seoane, P., Andrews, J., Arca Sedda, M., et al. 2022, arXiv e-prints, arXiv:2203.06016, doi: 10.48550/arXiv.2203.06016
  * Aramaki et al. (2020) Aramaki, T., Adrian, P. O. H., Karagiorgi, G., & Odaka, H. 2020, Astroparticle Physics, 114, 107, doi: 10.1016/j.astropartphys.2019.07.002
  * Artymowicz et al. (1993) Artymowicz, P., Lin, D. N. C., & Wampler, E. J. 1993, ApJ, 409, 592, doi: 10.1086/172690
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071, doi: 10.1088/0004-637X/697/2/1071
  * Bartos et al. (2017) Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017, ApJ, 835, 165\.  http://adsabs.harvard.edu/abs/2017ApJ...835..165B
  * Bartos et al. (2021) Bartos, I., Veske, D., Kowalski, M., Márka, Z., & Márka, S. 2021, ApJ, 921, 45, doi: 10.3847/1538-4357/ac1c7b
  * Bell (1978) Bell, A. R. 1978, MNRAS, 182, 147, doi: 10.1093/mnras/182.2.147
  * Berger (2014) Berger, E. 2014, ARA&A, 52, 43, doi: 10.1146/annurev-astro-081913-035926
  * Bland-Hawthorn & Cohen (2003) Bland-Hawthorn, J., & Cohen, M. 2003, ApJ, 582, 246, doi: 10.1086/344573
  * Blandford & Eichler (1987) Blandford, R., & Eichler, D. 1987, Phys. Rep., 154, 1, doi: 10.1016/0370-1573(87)90134-7
  * Blandford & Begelman (1999) Blandford, R. D., & Begelman, M. C. 1999, MNRAS, 303, L1, doi: 10.1046/j.1365-8711.1999.02358.x
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433, doi: 10.1093/mnras/179.3.433
  * Blasi (2013) Blasi, P. 2013, A&A Rev., 21, 70, doi: 10.1007/s00159-013-0070-7
  * Buitink et al. (2016) Buitink, S., Corstanje, A., Falcke, H., et al. 2016, Nature, 531, 70, doi: 10.1038/nature16976
  * Cao (2011) Cao, X. 2011, ApJ, 737, 94, doi: 10.1088/0004-637X/737/2/94
  * Caprioli (2015) Caprioli, D. 2015, ApJ, 811, L38, doi: 10.1088/2041-8205/811/2/L38
  * Caprioli et al. (2020) Caprioli, D., Haggerty, C. C., & Blasi, P. 2020, ApJ, 905, 2, doi: 10.3847/1538-4357/abbe05
  * Caprioli & Spitkovsky (2014) Caprioli, D., & Spitkovsky, A. 2014, ApJ, 783, 91, doi: 10.1088/0004-637X/783/2/91
  * Caputo et al. (2022) Caputo, R., Ajello, M., Kierans, C., et al. 2022, arXiv e-prints, arXiv:2208.04990.  https://arxiv.org/abs/2208.04990
  * Chang et al. (2008) Chang, P., Spitkovsky, A., & Arons, J. 2008, The Astrophysical Journal, 674, 378, doi: 10.1086/524764
  * Curd & Narayan (2022) Curd, B., & Narayan, R. 2022, arXiv e-prints, arXiv:2209.12081.  https://arxiv.org/abs/2209.12081
  * De Angelis et al. (2021) De Angelis, A., Tatischeff, V., Argan, A., et al. 2021, Experimental Astronomy, 51, 1225, doi: 10.1007/s10686-021-09706-y
  * Do et al. (2020) Do, T., David Martinez, G., Kerzendorf, W., et al. 2020, ApJ, 901, L28, doi: 10.3847/2041-8213/abb246
  * Eichmann & Becker Tjus (2016) Eichmann, B., & Becker Tjus, J. 2016, ApJ, 821, 87, doi: 10.3847/0004-637X/821/2/87
  * Eichmann et al. (2022) Eichmann, B., Oikonomou, F., Salvatore, S., Dettmar, R.-J., & Becker Tjus, J. 2022, arXiv e-prints, arXiv:2207.00102.  https://arxiv.org/abs/2207.00102
  * Fan & Piran (2008) Fan, Y.-Z., & Piran, T. 2008, Frontiers of Physics in China, 3, 306, doi: 10.1007/s11467-008-0033-z
  * Finkbeiner (2004) Finkbeiner, D. P. 2004, ApJ, 614, 186, doi: 10.1086/423482
  * Fouka & Ouichaoui (2009) Fouka, M., & Ouichaoui, S. 2009, ApJ, 707, 278, doi: 10.1088/0004-637X/707/1/278
  * Fouka & Ouichaoui (2011) —. 2011, ApJ, 743, 89, doi: 10.1088/0004-637X/743/1/89
  * Frail et al. (2005) Frail, D. A., Soderberg, A. M., Kulkarni, S. R., et al. 2005, ApJ, 619, 994, doi: 10.1086/426680
  * Funk et al. (2013) Funk, S., Hinton, J. A., & CTA Consortium. 2013, Astroparticle Physics, 43, 348, doi: 10.1016/j.astropartphys.2012.05.018
  * García-Burillo et al. (2016) García-Burillo, S., Combes, F., Ramos Almeida, C., et al. 2016, ApJ, 823, L12, doi: 10.3847/2041-8205/823/1/L12
  * Gilli et al. (2007) Gilli, R., Comastri, A., & Hasinger, G. 2007, A&A, 463, 79, doi: 10.1051/0004-6361:20066334
  * Giommi & Padovani (2015) Giommi, P., & Padovani, P. 2015, MNRAS, 450, 2404, doi: 10.1093/mnras/stv793
  * Goodman & Tan (2004) Goodman, J., & Tan, J. C. 2004, ApJ, 608, 108.  http://adsabs.harvard.edu/abs/2004ApJ...608..108G
  * Graham et al. (2020) Graham, M. J., Ford, K. E. S., McKernan, B., et al. 2020, Phys. Rev. Lett., 124, 251102, doi: 10.1103/PhysRevLett.124.251102
  * Greenhill & Gwinn (1997) Greenhill, L. J., & Gwinn, C. R. 1997, Ap&SS, 248, 261, doi: 10.1023/A:1000554317683
  * Grishin et al. (2021) Grishin, E., Bobrick, A., Hirai, R., Mandel, I., & Perets, H. B. 2021, MNRAS, 507, 156, doi: 10.1093/mnras/stab1957
  * Guo & Mathews (2012) Guo, F., & Mathews, W. G. 2012, The Astrophysical Journal, 756, 181, doi: 10.1088/0004-637x/756/2/181
  * Gutiérrez et al. (2021) Gutiérrez, E. M., Vieyro, F. L., & Romero, G. E. 2021, A&A, 649, A87, doi: 10.1051/0004-6361/202039671
  * Hada et al. (2013) Hada, K., Kino, M., Doi, A., et al. 2013, ApJ, 775, 70, doi: 10.1088/0004-637X/775/1/70
  * Hada et al. (2018) Hada, K., Doi, A., Wajima, K., et al. 2018, ApJ, 860, 141, doi: 10.3847/1538-4357/aac49f
  * Hayashida et al. (2013) Hayashida, M., Stawarz, Ł., Cheung, C. C., et al. 2013, ApJ, 779, 131, doi: 10.1088/0004-637X/779/2/131
  * Hopkins et al. (2012) Hopkins, P. F., Hernquist, L., Hayward, C. C., & Narayanan, D. 2012, MNRAS, 425, 1121, doi: 10.1111/j.1365-2966.2012.21449.x
  * Hu et al. (2022) Hu, Y., Xu, S., Stone, J. M., & Lazarian, A. 2022, ApJ, 941, 133, doi: 10.3847/1538-4357/ac9ebc
  * IceCube Collaboration et al. (2022) IceCube Collaboration, Abbasi, R., Ackermann, M., et al. 2022, Science, 378, 538, doi: 10.1126/science.abg3395
  * Imanishi et al. (2018) Imanishi, M., Nakanishi, K., Izumi, T., & Wada, K. 2018, ApJ, 853, L25, doi: 10.3847/2041-8213/aaa8df
  * Imanishi et al. (2020) Imanishi, M., Nguyen, D. D., Wada, K., et al. 2020, ApJ, 902, 99, doi: 10.3847/1538-4357/abaf50
  * Impellizzeri et al. (2019) Impellizzeri, C. M. V., Gallimore, J. F., Baum, S. A., et al. 2019, ApJ, 884, L28, doi: 10.3847/2041-8213/ab3c64
  * Inoue et al. (2022) Inoue, S., Cerruti, M., Murase, K., & Liu, R.-Y. 2022, arXiv e-prints, arXiv:2207.02097.  https://arxiv.org/abs/2207.02097
  * Inoue (2014) Inoue, Y. 2014, arXiv e-prints, arXiv:1412.3886.  https://arxiv.org/abs/1412.3886
  * Inoue et al. (2020) Inoue, Y., Khangulyan, D., & Doi, A. 2020, ApJ, 891, L33, doi: 10.3847/2041-8213/ab7661
  * Inoue et al. (2013) Inoue, Y., Murase, K., Madejski, G. M., & Uchiyama, Y. 2013, ApJ, 776, 33, doi: 10.1088/0004-637X/776/1/33
  * Ito et al. (2020) Ito, H., Levinson, A., & Nakar, E. 2020, MNRAS, 499, 4961, doi: 10.1093/mnras/staa3125
  * Keshet et al. (2009) Keshet, U., Katz, B., Spitkovsky, A., & Waxman, E. 2009, The Astrophysical Journal, 693, L127, doi: 10.1088/0004-637x/693/2/l127
  * Kimura (2022) Kimura, S. S. 2022, arXiv e-prints, arXiv:2202.06480.  https://arxiv.org/abs/2202.06480
  * Kimura et al. (2021a) Kimura, S. S., Murase, K., & Bartos, I. 2021a, ApJ, 916, 111, doi: 10.3847/1538-4357/ac0535
  * Kimura et al. (2019) Kimura, S. S., Murase, K., Ioka, K., et al. 2019, ApJ, 887, L16, doi: 10.3847/2041-8213/ab59e1
  * Kimura et al. (2021b) Kimura, S. S., Murase, K., & Mészáros, P. 2021b, Nature Communications, 12, 5615, doi: 10.1038/s41467-021-25111-7
  * Kimura et al. (2015) Kimura, S. S., Murase, K., & Toma, K. 2015, ApJ, 806, 159, doi: 10.1088/0004-637X/806/2/15910.48550/arXiv.1411.3588
  * Kimura et al. (2018) Kimura, S. S., Murase, K., & Zhang, B. T. 2018, Phys. Rev. D, 97, 023026, doi: 10.1103/PhysRevD.97.023026
  * Kimura et al. (2021c) Kimura, S. S., Sudoh, T., Kashiyama, K., & Kawanaka, N. 2021c, ApJ, 915, 31, doi: 10.3847/1538-4357/abff58
  * Kimura & Toma (2020) Kimura, S. S., & Toma, K. 2020, ApJ, 905, 178, doi: 10.3847/1538-4357/abc343
  * Kinney et al. (2000) Kinney, A. L., Schmitt, H. R., Clarke, C. J., et al. 2000, ApJ, 537, 152, doi: 10.1086/309016
  * Kocsis & Tremaine (2011) Kocsis, B., & Tremaine, S. 2011, MNRAS, 412, 187, doi: 10.1111/j.1365-2966.2010.17897.x
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Lacki et al. (2014) Lacki, B. C., Horiuchi, S., & Beacom, J. F. 2014, ApJ, 786, 40, doi: 10.1088/0004-637X/786/1/40
  * Levin (2007) Levin, Y. 2007, MNRAS, 374, 515, doi: 10.1111/j.1365-2966.2006.11155.x
  * Levin (2022) —. 2022, arXiv e-prints, arXiv:2211.12754, doi: 10.48550/arXiv.2211.12754
  * Levin & Beloborodov (2003) Levin, Y., & Beloborodov, A. M. 2003, ApJ, 590, L33, doi: 10.1086/376675
  * Liska et al. (2020) Liska, M., Tchekhovskoy, A., & Quataert, E. 2020, MNRAS, 494, 3656, doi: 10.1093/mnras/staa955
  * Lu et al. (2013) Lu, J. R., Do, T., Ghez, A. M., et al. 2013, ApJ, 764, 155.  http://adsabs.harvard.edu/abs/2013ApJ...764..155L
  * Marinucci et al. (2016) Marinucci, A., Bianchi, S., Matt, G., et al. 2016, MNRAS, 456, L94, doi: 10.1093/mnrasl/slv178
  * Martins et al. (2009) Martins, S. F., Fonseca, R. A., Silva, L. O., & Mori, W. B. 2009, ApJ, 695, L189, doi: 10.1088/0004-637X/695/2/L189
  * Matsumoto et al. (2020) Matsumoto, T., Kimura, S. S., Murase, K., & Mészáros, P. 2020, MNRAS, 493, 783, doi: 10.1093/mnras/staa305
  * McKernan et al. (2012) McKernan, B., Ford, K. E. S., Lyra, W., & Perets, H. B. 2012, MNRAS, 425, 460, doi: 10.1111/j.1365-2966.2012.21486.x
  * McKernan et al. (2020) McKernan, B., Ford, K. E. S., & O’Shaughnessy, R. 2020, MNRAS, 498, 4088, doi: 10.1093/mnras/staa2681
  * McKernan et al. (2018) McKernan, B., Ford, K. E. S., Bellovary, J., et al. 2018, ApJ, 866, 66, doi: 10.3847/1538-4357/aadae5
  * McKernan et al. (2019) McKernan, B., Ford, K. E. S., Bartos, I., et al. 2019, ApJ, 884, L50, doi: 10.3847/2041-8213/ab4886
  * Medvedev & Loeb (1999) Medvedev, M. V., & Loeb, A. 1999, ApJ, 526, 697, doi: 10.1086/308038
  * Meier (2001) Meier, D. L. 2001, ApJ, 548, L9, doi: 10.1086/318921
  * Merritt (2010) Merritt, D. 2010, ApJ, 718, 739, doi: 10.1088/0004-637X/718/2/739
  * Miller et al. (2019) Miller, R., Ajello, M., Beacom, J. F., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 123.  https://arxiv.org/abs/1907.07005
  * Milosavljević & Loeb (2004) Milosavljević, M., & Loeb, A. 2004, ApJ, 604, L45, doi: 10.1086/383467
  * Miralda-Escude & Gould (2000) Miralda-Escude, J., & Gould, A. 2000, ApJ, 545, 847.  http://adsabs.harvard.edu/abs/2000ApJ...545..847M
  * Moranchel-Basurto et al. (2021) Moranchel-Basurto, A., Sánchez-Salcedo, F. J., Chametla, R. O., & Velázquez, P. F. 2021, ApJ, 906, 15, doi: 10.3847/1538-4357/abca88
  * Murase et al. (2016) Murase, K., Guetta, D., & Ahlers, M. 2016, Phys. Rev. Lett., 116, 071101, doi: 10.1103/PhysRevLett.116.071101
  * Murase et al. (2014) Murase, K., Inoue, Y., & Dermer, C. D. 2014, Phys. Rev. D, 90, 023007, doi: 10.1103/PhysRevD.90.023007
  * Murase et al. (2020) Murase, K., Kimura, S. S., & Mészáros, P. 2020, Phys. Rev. Lett., 125, 011101, doi: 10.1103/PhysRevLett.125.011101
  * Murase & Waxman (2016) Murase, K., & Waxman, E. 2016, Phys. Rev. D, 94, 103006, doi: 10.1103/PhysRevD.94.103006
  * Narayan et al. (2022) Narayan, R., Chael, A., Chatterjee, K., Ricarte, A., & Curd, B. 2022, MNRAS, 511, 3795, doi: 10.1093/mnras/stac285
  * Narayan & Yi (1994) Narayan, R., & Yi, I. 1994, ApJ, 428, L13, doi: 10.1086/187381
  * Neumayer et al. (2020) Neumayer, N., Seth, A., & Böker, T. 2020, A&A Rev., 28, 4, doi: 10.1007/s00159-020-00125-0
  * Ostriker (1983) Ostriker, J. P. 1983, ApJ, 273, 99, doi: 10.1086/161351
  * Panaitescu & Kumar (2001) Panaitescu, A., & Kumar, P. 2001, ApJ, 560, L49, doi: 10.1086/324061
  * Pe’Er et al. (2009) Pe’Er, A., Murase, K., & Mészáros, P. 2009, Phys. Rev. D, 80, 123018, doi: 10.1103/PhysRevD.80.123018
  * Perna et al. (2021a) Perna, R., Lazzati, D., & Cantiello, M. 2021a, ApJ, 906, L7, doi: 10.3847/2041-8213/abd319
  * Perna et al. (2021b) Perna, R., Tagawa, H., Haiman, Z., & Bartos, I. 2021b, ApJ, 915, 10, doi: 10.3847/1538-4357/abfdb4
  * Pier et al. (1994) Pier, E. A., Antonucci, R., Hurt, T., Kriss, G., & Krolik, J. 1994, ApJ, 428, 124, doi: 10.1086/174225
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Pushkarev et al. (2009) Pushkarev, A. B., Kovalev, Y. Y., Lister, M. L., & Savolainen, T. 2009, A&A, 507, L33, doi: 10.1051/0004-6361/200913422
  * Razzaque et al. (2004) Razzaque, S., Mészáros, P., & Waxman, E. 2004, Phys. Rev. Lett., 93, 181101, doi: 10.1103/PhysRevLett.93.181101
  * Revnivtsev et al. (2003) Revnivtsev, M., Gilfanov, M., Sunyaev, R., Jahoda, K., & Markwardt, C. 2003, A&A, 411, 329, doi: 10.1051/0004-6361:20031386
  * Rybicki & Lightman (1979) Rybicki, G. B., & Lightman, A. P. 1979, Radiative processes in astrophysics (Wiley and Sons, Incorporated, John)
  * Santana et al. (2014) Santana, R., Barniol Duran, R., & Kumar, P. 2014, ApJ, 785, 29, doi: 10.1088/0004-637X/785/1/29
  * Sari et al. (1998) Sari, R., Piran, T., & Narayan, R. 1998, ApJ, 497, L17, doi: 10.1086/311269
  * Shields (1996) Shields, G. A. 1996, ApJ, 461, L9, doi: 10.1086/310001
  * Sironi et al. (2013) Sironi, L., Spitkovsky, A., & Arons, J. 2013, The Astrophysical Journal, 771, 54, doi: 10.1088/0004-637x/771/1/54
  * Spitkovsky (2008) Spitkovsky, A. 2008, ApJ, 673, L39, doi: 10.1086/527374
  * Sridhar et al. (2022) Sridhar, N., Metzger, B. D., & Fang, K. 2022, arXiv e-prints, arXiv:2212.11236.  https://arxiv.org/abs/2212.11236
  * Stone et al. (2017) Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464, 946.  http://adsabs.harvard.edu/abs/2017MNRAS.464..946S
  * Su et al. (2010) Su, M., Slatyer, T. R., & Finkbeiner, D. P. 2010, The Astrophysical Journal, 724, 1044, doi: 10.1088/0004-637x/724/2/1044
  * Syer et al. (1991) Syer, D., Clarke, C. J., & Rees, M. J. 1991, MNRAS, 250, 505, doi: 10.1093/mnras/250.3.505
  * Tagawa et al. (2020a) Tagawa, H., Haiman, Z., Bartos, I., & Kocsis, B. 2020a, ApJ, 899, 26, doi: 10.3847/1538-4357/aba2cc
  * Tagawa et al. (2020b) Tagawa, H., Haiman, Z., & Kocsis, B. 2020b, ApJ, 898, 25, doi: 10.3847/1538-4357/ab9b8c
  * Tagawa et al. (2023a) Tagawa, H., Kimura, S. S., Haiman, Z., Perna, R., & Bartos, I. 2023a, arXiv e-prints, arXiv:2301.07111, doi: 10.48550/arXiv.2301.07111
  * Tagawa et al. (2023b) —. 2023b, ApJ, 946, L3, doi: 10.3847/2041-8213/acc103
  * Tagawa et al. (2022) Tagawa, H., Kimura, S. S., Haiman, Z., et al. 2022, ApJ, 927, 41, doi: 10.3847/1538-4357/ac45f8
  * Tagawa et al. (2021) Tagawa, H., Kocsis, B., Haiman, Z., et al. 2021, ApJ, 908, 194, doi: 10.3847/1538-4357/abd555
  * Takami et al. (2016) Takami, H., Murase, K., & Dermer, C. D. 2016, ApJ, 817, 59, doi: 10.3847/0004-637X/817/1/59
  * Takami & Sato (2009) Takami, H., & Sato, K. 2009, Astroparticle Physics, 30, 306, doi: 10.1016/j.astropartphys.2008.10.004
  * Tchekhovskoy et al. (2010) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2010, ApJ, 711, 50, doi: 10.1088/0004-637X/711/1/50
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2021, arXiv e-prints, arXiv:2111.03606.  https://arxiv.org/abs/2111.03606
  * Thompson et al. (2005) Thompson, T. A., Quataert, E., & Murray, N. 2005, ApJ, 630, 167.  http://adsabs.harvard.edu/abs/2005ApJ...630..167T
  * Tomita et al. (2022) Tomita, S., Ohira, Y., Kimura, S. S., Tomida, K., & Toma, K. 2022, ApJ, 936, L9, doi: 10.3847/2041-8213/ac88be
  * Tomita et al. (2019) Tomita, S., Ohira, Y., & Yamazaki, R. 2019, ApJ, 886, 54, doi: 10.3847/1538-4357/ab4a10
  * Tomsick et al. (2019) Tomsick, J., Zoglauer, A., Sleator, C., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 98.  https://arxiv.org/abs/1908.04334
  * Troja et al. (2019) Troja, E., van Eerten, H., Ryan, G., et al. 2019, MNRAS, 489, 1919, doi: 10.1093/mnras/stz2248
  * Uchiyama et al. (2007) Uchiyama, Y., Aharonian, F. A., Tanaka, T., Takahashi, T., & Maeda, Y. 2007, Nature, 449, 576, doi: 10.1038/nature06210
  * Ueda et al. (2014) Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., & Watson, M. G. 2014, ApJ, 786, 104, doi: 10.1088/0004-637X/786/2/104
  * Wang et al. (2021a) Wang, J.-M., Liu, J.-R., Ho, L. C., & Du, P. 2021a, ApJ, 911, L14, doi: 10.3847/2041-8213/abee81
  * Wang et al. (2021b) Wang, J.-M., Liu, J.-R., Ho, L. C., Li, Y.-R., & Du, P. 2021b, ApJ, 916, L17, doi: 10.3847/2041-8213/ac0b46
  * Watanabe et al. (1997) Watanabe, K., Hartmann, D. H., Leising, M. D., et al. 1997, in American Institute of Physics Conference Series, Vol. 410, Proceedings of the Fourth Compton Symposium, ed. C. D. Dermer, M. S. Strickman, & J. D. Kurfess, 1223–1227, doi: 10.1063/1.53933
  * Waxman & Bahcall (1998) Waxman, E., & Bahcall, J. 1998, Phys. Rev. D, 59, 023002, doi: 10.1103/PhysRevD.59.023002
  * Waxman & Loeb (1999) Waxman, E., & Loeb, A. 1999, ApJ, 515, 721, doi: 10.1086/307066
  * Weidenspointner et al. (2000) Weidenspointner, G., Varendorff, M., Kappadath, S. C., et al. 2000, in American Institute of Physics Conference Series, Vol. 510, The Fifth Compton Symposium, ed. M. L. McConnell & J. M. Ryan, 467–470, doi: 10.1063/1.1307028
  * Wojaczyński & Niedźwiecki (2017) Wojaczyński, R., & Niedźwiecki, A. 2017, ApJ, 849, 97, doi: 10.3847/1538-4357/aa8f9d
  * Wojaczyński et al. (2015) Wojaczyński, R., Niedźwiecki, A., Xie, F.-G., & Szanecki, M. 2015, A&A, 584, A20, doi: 10.1051/0004-6361/201526621
  * Xiong & Zhang (2014) Xiong, D. R., & Zhang, X. 2014, MNRAS, 441, 3375, doi: 10.1093/mnras/stu755
  * Xu et al. (2018) Xu, F., Bian, F., Shen, Y., et al. 2018, MNRAS, 480, 345, doi: 10.1093/mnras/sty1763
  * Yang et al. (2018) Yang, H. Y., Ruszkowski, M., & Zweibel, E. 2018, Galaxies, 6, 29, doi: 10.3390/galaxies6010029
  * Yang et al. (2012) Yang, H. Y. K., Ruszkowski, M., Ricker, P. M., Zweibel, E., & Lee, D. 2012, ApJ, 761, 185, doi: 10.1088/0004-637X/761/2/185
  * Yang et al. (2022a) Yang, H. Y. K., Ruszkowski, M., & Zweibel, E. G. 2022a, Nature Astronomy, doi: 10.1038/s41550-022-01618-x
  * Yang et al. (2022b) Yang, Y., Bartos, I., Fragione, G., et al. 2022b, ApJ, 933, L28, doi: 10.3847/2041-8213/ac7c0b
  * Yang et al. (2019) Yang, Y., Bartos, I., Gayathri, V., et al. 2019, Phys. Rev. Lett., 123, 181101, doi: 10.1103/PhysRevLett.123.181101
  * Yelda et al. (2014) Yelda, S., Ghez, A. M., Lu, J. R., et al. 2014, ApJ, 783, 131.  http://adsabs.harvard.edu/abs/2014ApJ...783..131Y
  * Yuan et al. (2022) Yuan, C., Murase, K., Guetta, D., et al. 2022, ApJ, 932, 80, doi: 10.3847/1538-4357/ac6ddf
  * Zaino et al. (2020) Zaino, A., Bianchi, S., Marinucci, A., et al. 2020, MNRAS, 492, 3872, doi: 10.1093/mnras/staa107
  * Zhu et al. (2021a) Zhu, J.-P., Wang, K., & Zhang, B. 2021a, ApJ, 917, L28, doi: 10.3847/2041-8213/ac1a17
  * Zhu et al. (2021b) Zhu, J.-P., Yang, Y.-P., Zhang, B., et al. 2021b, ApJ, 914, L19, doi: 10.3847/2041-8213/abff5a
  * Zhu et al. (2021c) Zhu, J.-P., Zhang, B., Yu, Y.-W., & Gao, H. 2021c, ApJ, 906, L11, doi: 10.3847/2041-8213/abd412
  * Zirakashvili et al. (2022) Zirakashvili, V. N., Ptuskin, V. S., & Rogovaya, S. I. 2022, MNRAS, doi: 10.1093/mnrasl/slac142
  * Zweibel (2013) Zweibel, E. G. 2013, Physics of Plasmas, 20, 055501, doi: 10.1063/1.4807033