#  Nuclear and electromagnetic cascades induced by ultrahigh-energy cosmic
rays in radio galaxies: implications for Centaurus A

B. Theodore Zhang (张兵),1 Kohta Murase2,3,4,5,1
1Center for Gravitational Physics and Quantum Information, Yukawa Institute
for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan
2Department of Physics, The Pennsylvania State University, University Park, PA
16802, USA
3Department of Astronomy & Astrophysics, The Pennsylvania State University,
University Park, PA 16802, USA
4Center for Multimessenger Astrophysics, Institute for Gravitation and the
Cosmos, The Pennsylvania State University, University Park, PA 16802, USA
5School of Natural Sciences, Institute for Advanced Study, Princeton, NJ
08540, USA
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Very-high-energy (VHE) $\gamma$-rays ($\gtrsim 0.1\rm~{}TeV$) and neutrinos
are crucial for identifying accelerators of ultrahigh-energy cosmic rays
(UHECRs), but this is challenging especially for UHECR nuclei. In this work,
we develop a numerical code to solve the transport equation for UHECRs and
their secondaries, where both nuclear and electromagnetic cascades are taken
into account self-consistently, considering steady UHECR accelerators such as
radio galaxies. In particular, we focus on Centaurus A, which has been
proposed as one of the most promising UHECR sources in the local universe.
Motivated by observations of extended VHE $\gamma$-ray emission from its
kiloparsec-scale jet by the H.E.S.S. telescope, we study interactions between
UHECRs accelerated in the large-scale jet and various target photon fields
including blazar-like beamed core emission, and present a quantitative study
on VHE $\gamma$-ray signatures of UHECR nuclei, including the
photodisintegration and Bethe-Heitler pair-production processes. We show that
VHE $\gamma$-rays from UHECR nuclei could be detected by the ground-based
$\gamma$-ray telescopes given that the dominant composition of UHECRs consists
of intermediate-mass (such as oxygen) nuclei.

###### keywords:

galaxies: jets - gamma-rays: galaxies - radio continuum: galaxies -
acceleration of particles - radiation mechanisms: non-thermal.

††pubyear: 2022††pagerange:  Nuclear and electromagnetic cascades induced by
ultrahigh-energy cosmic rays in radio galaxies: implications for Centaurus
A–A.6

##  1 Introduction

Very-high-energy (VHE) $\gamma$-rays ($\gtrsim 0.1\rm~{}TeV$) from
extragalactic sources have been detected by ground-based Imaging Atmospheric
Cherenkov Telescopes (IACTs) and surface detectors (e.g., Hinton & Hofmann,
2009; Sitarek, 2022). Based on the unification model of radio-loud AGNs,
blazars including BL Lac objects and flat-spectrum radio quasars (FSRQs) have
relativistic jets pointing along the line-of-sight, while radio galaxies
including Fanaroff-Riley (FR) I and II galaxies are viewed with an inclination
angle to the jet axis where the effect of Doppler boosting is modest (e.g.,
Urry & Padovani, 1995). In the local universe, several radio-bright galaxies
have been detected in the VHE band, including Centaurus A (Cen A), M87, NGC
1275, IC 310, 3C 264, and PKS 0625-35 (Rieger & Levinson, 2018; Rulten, 2022).
Cen A, the closest radio galaxy located at $\sim 3.7\rm~{}Mpc$ from Earth, has
been extensively studied due to its potential as a source of ultrahigh-energy
cosmic rays (UHECRs, $\gtrsim 1\rm~{}EeV$) detected (e.g., Kotera & Olinto,
2011; Anchordoqui, 2019; Alves Batista et al., 2019; Rieger, 2022, for
reviews). Additionally, the Pierre Auger Collaboration observed a $\sim
4\sigma$ deviation from isotropy at the intermediate angular scales in the
Centaurus region for UHECRs with energies beyond $\sim 40\rm~{}EeV$ (Biteau et
al., 2021). Cen A was also the first extragalactic extended source detected in
the GeV sky by the Fermi-LAT (The Fermi-LAT Collaboration et al., 2010). H. E.
S. S. Collaboration et al. (2020) has provided strong evidence for VHE
$\gamma$-ray production from components beyond the inner core of the radio
galaxy (Aharonian et al., 2009; Abdalla et al., 2018). However, current
instruments still face limitations in resolving the emission regions from
other VHE radio galaxies (Rulten, 2022).

The production of VHE $\gamma$-rays from radio galaxies, such as Cen A, is
thought to be the result of inverse-Compton (IC) emission by non-thermal
electrons accelerated in the energy dissipation region within the jet. This
mechanism is consistent with the available evidence (e.g., H. E. S. S.
Collaboration et al., 2020). However, hadronic scenarios, in which VHE
$\gamma$-rays are produced through the interaction of high-energy cosmic rays,
are still a viable possibility (e.g., Petropoulou et al., 2014; Fraija, 2014;
Fraija & Marinelli, 2016). The hadronic components of the VHE $\gamma$-rays
include electromagnetic cascade emission from photomeson production,
hadronuclear interaction, Bethe-Heitler electron-position production, and
photodisintegration accompanied by de-excitation $\gamma$-rays. In magnetized
environments, synchrotron radiation from UHECR protons and nuclei could also
contribute significantly to the observed VHE $\gamma$-rays (e.g., Mücke &
Protheroe, 2001; Murase et al., 2008). The detection of de-excitation VHE
$\gamma$-rays from nearby UHECR sources may provide direct evidence of the
acceleration of the heavier nuclei components of UHECRs (Murase & Beacom,
2010b). Based on observations showing that the fraction of heavier nuclei
increases beyond $\sim 4$ EeV (e.g., Guido et al., 2021), it is likely that
the highest energy range of UHECRs accelerated in nearby sources is dominated
by heavier nuclei.

The photodisintegration of UHECR nuclei is a crucial process in determining
their fate, especially when target photon energy in the nuclear rest frame
exceeds the nuclear binding energy of approximately 10 MeV. This leads to the
ejection of one or several nucleons, as described by Stecker (1969). Even
fragmentation may occur when the target photon energy is high, and the
photomeson production may occur when it exceeds the pion production threshold
of approximately 140 MeV. As a result of photodisintegration, the nuclear
fragments are left in an excited state, which quickly de-excites through the
emission of one or several photons with energies around MeV in the nuclear
rest frame. In the observer frame, these de-excitation $\gamma$-rays are
boosted to the VHE range for ultra-relativistic cosmic ray nuclei. The process
was proposed for not only Galactic point sources (Karakula et al., 1994;
Anchordoqui et al., 2007a, b) but also extragalactic sources (Murase & Beacom,
2010b). It has been shown that UHECR nuclei can survive in radio galaxies like
Cen A (Murase et al., 2012), implying that photodisintegration should not be
efficient in regions where UHECRs are produced. For a given target photon
energy the efficiency of de-excitation is also lower than those of the
photomeson and Bethe-Heitler pair production, but TeV gamma rays can still be
dominated by de-excitation $\gamma$-rays (see, e.g., Fig. 1 of Murase &
Beacom, 2010b). The potential production of de-excitation VHE $\gamma$-rays
from the core of Cen A via photodisintegration of heavy nuclei has been
discussed (Murase & Beacom, 2010b; Kundu & Gupta, 2014; Morejon et al., 2021).

In this study, we present the numerical framework implemented in the
Astrophysical Multimessenger Emission Simulator (AMES). The framework can
simultaneously treat both the nuclear cascade and the electromagnetic cascade
by solving the coupled transport equations. The fate of UHECR nuclei in
various astrophysical sources has been widely explored in previous studies
(e.g., Rodrigues et al., 2018; Zhang et al., 2018; Biehl et al., 2018;
Boncioli et al., 2019; Zhang & Murase, 2019). Our goal is to examine the
impact of electromagnetic cascades on the multi-wavelength spectral energy
distribution (SED), considering the injection of UHECR nuclei. We apply our
code to model leptohadronic processes in the large-scale jet of Cen A, where
we investigate the detectability of de-excitation $\gamma$-rays. In contrast
to the previous work on blazars, we consider the acceleration zone located in
the large-scale jet, motivated by recent studies from the H.E.S.S.
Collaboration (H. E. S. S. Collaboration et al., 2020). Additionally, we take
into account the beamed photons from the inner core as the dominant target
photons (e.g., Bednarek, 2019; Sudoh et al., 2020).

The structure of the paper is as follows: In Sec. 2, we present an in-depth
explanation of the physical processes related to the modeling of nuclear and
electromagnetic cascades. In Sec. 3, we use our numerical framework to
investigate the hadronic origin of the VHE $\gamma$-rays detected from Cen A
and the feasibility of detecting de-excitation $\gamma$-rays with present and
future ground-based $\gamma$-ray detectors. Sec. 4 explores the implications
of our results. Finally, in Sec. 5, we summarize the work.

##  2 Physical processes related to nuclear and electromagnetic cascade

To begin, let us consider a basic physical model where all physical processes
occur within a uniform spherical emission region with comoving radius $l_{b}$
and radius from the black hole $R$. This emission region is encompassed by
tangled magnetic fields with magnetic field strength $B$, and Doppler factor
$\delta_{D}$. In the following, we adopt the notation
$E=\delta_{D}\varepsilon$, where $E$ represents the particle energy measured
in the observer frame and $\varepsilon$ represents the particle energy
measured in the comoving frame. Note the redshift evolution factor $(1+z)$ is
not included considering the source is nearby.

We assume the injection of CR nuclei following a power-law distribution with
an exponential cutoff,

$q(\varepsilon_{A},t)\equiv\frac{dN^{A}}{dtd\varepsilon_{A}}=\dot{N}_{0}^{A}\left(\frac{\varepsilon_{A}}{Z\varepsilon_{A}^{0}}\right)^{-s_{\rm
acc}}{{\rm exp}\left(-\frac{\varepsilon_{A}}{Z\varepsilon_{p}^{\rm
max}}\right)},$ (1)

where $A$ is the CR nuclear mass number, $Z$ is the charge number, $s_{\rm
acc}$ is the acceleration spectral index, $\varepsilon_{p}^{\rm max}$ is the
proton maximum acceleration energy measured in the comoving frame,
$\dot{N}_{0}^{A}$ is the normalization constant with units [$\rm eV^{-1}\
s^{-1}$], which is determined by

$L_{\rm CR}=\frac{3}{2}\Gamma^{2}\sum_{A}\int
d\varepsilon_{A}\frac{dN^{A}}{dtd\varepsilon_{A}},$ (2)

where $L_{\rm CR}$ is the total CR injection luminosity measured in the source
frame.

We model the leptohadronic processes by solving a series of coupled transport
equations for various particles including photons, electrons, neutrinos,
neutrons, protons, and nuclei,

$\displaystyle\frac{\partial n_{\varepsilon_{a}}^{a}}{\partial t}$
$\displaystyle=-n_{\varepsilon_{a}}^{a}\mathcal{A}_{a}(\varepsilon_{a})$
$\displaystyle+\int
d\varepsilon_{a}^{*}n_{\varepsilon_{a}^{*}}^{a}\mathcal{B}_{a\to
a}(\varepsilon_{a},\varepsilon_{a}^{*})$ $\displaystyle+\sum_{b}\int
d\varepsilon_{b}n_{\varepsilon_{b}}^{b}\mathcal{C}_{b\to
a}(\varepsilon_{a},\varepsilon_{b})$ $\displaystyle+\dot{n}_{a}^{\rm
inj}(\varepsilon_{a}),$ (3)

where $a$ is particle type, $n_{\varepsilon_{a}}^{a}$ is particle differential
number density at energy $\varepsilon_{a}$, $\mathcal{A}_{a}(\varepsilon_{a})$
is the total interaction rate at energy $\varepsilon_{a}$ including particle
escape, $\mathcal{B}_{a\to a}(\varepsilon_{a},\varepsilon_{a}^{*})$ is the
self-production rate of particles with energy $\varepsilon_{a}$ generated from
the same type of particles with energy $\varepsilon_{a}^{*}$, and
$\mathcal{C}_{b\to a(\varepsilon_{a},\varepsilon_{b})}$ is the generation rate
of particles with energy $\varepsilon_{a}$ from other types of particles with
energy $\varepsilon_{b}$, and $\dot{n}_{a}^{\rm inj}$ is the source injection
rate at energy $\varepsilon_{a}$. In Appendix A, we provide the details of the
numerical scheme and the related physical processes described in the above
sections for each particle species, respectively. We here provide a part of
the Astrophysical Multimessenger Emission Simulator (AMES) code, which
integrates those for hadronic and leptonic emissions from various
astrophysical objects. Some of the earlier calculations without nuclear
cascades are found in, e.g., Murase (2012); Murase & Beacom (2012); Murase et
al. (2015); Murase (2018, 2022).

The photonuclear interaction rate for CR nuclei can be calculated by

$\displaystyle t_{A\gamma}^{-1}$ $\displaystyle=c\oint
d\Omega\int_{0}^{\infty}d\varepsilon(1-\beta_{A}\mu)\frac{dn}{d\varepsilon
d\Omega}\sigma_{A\gamma}(\bar{\varepsilon})$
$\displaystyle=\frac{c}{2}\int_{-1}^{1}d\mu\int_{0}^{\infty}d\varepsilon(1-\beta_{A}\mu)\frac{dn}{d\varepsilon}\sigma_{A\gamma}(\bar{\varepsilon}),$
(4)

where $c$ is the speed of light, the solid angle averaged target photon
distribution $dn/d\varepsilon\equiv(1/4\pi)\oint d\Omega(dn/d\varepsilon
d\Omega)$, $\beta_{A}$ is the particle velocity, $\mu=\rm cos\theta$ is the
angle between nuclei and incident target photon, $\sigma_{A\gamma}$ is the
cross section, and $\bar{\varepsilon}=\gamma_{A}(1-\beta_{A}\mu)\varepsilon$
is the target photon energy measured in the nuclear rest frame. We can define
optical depth as

$\tau_{A\gamma}\approx t_{\rm esc}/t_{A\gamma},$ (5)

where $t_{\rm esc}$ is the escape time and in the limit that it is dominated
by advection we have

$t_{\rm esc}\approx t_{\rm adv}\approx\frac{l_{b}}{V},$ (6)

where $t_{\rm adv}$ is the advection time scale.

Note that the energy loss rate, which is useful in analytical estimates, can
be calculated similarly, by considering the inelasticity $\kappa_{A\gamma}$.

For the purpose of understanding physical processes analytically, let us
consider target photon fields with a broken power law,

$\frac{dn}{d\varepsilon}=n_{\varepsilon_{b}}\begin{cases}(\varepsilon/\varepsilon_{b})^{-\alpha_{l}},\varepsilon<\varepsilon_{b}\\\
(\varepsilon/\varepsilon_{b})^{-\alpha_{h}},\varepsilon>\varepsilon_{b}\end{cases},$
(7)

in the comoving frame of the blob, where $n_{\varepsilon_{b}}$ is the
differential photon energy density at $\varepsilon_{b}$ with units [$\rm
eV^{-1}~{}cm^{-3}$], $\varepsilon_{b}$ is the break energy, $\alpha_{l}$ and
$\alpha_{h}$ are spectral indices.

In Fig. 1, we show a schematic picture of interactions by CR nuclei and
electrons, including nuclear and electromagnetic cascades induced by CR
nuclei, in the large-scale jet of radio galaxies. The relevant physical
processes include photodisintegration of nuclei, photomeson production
process, Bethe-Heitler pair production process, and accompanied
electromagnetic cascades.

Figure 1: Schematic picture of neutrino and gamma-ray production processes
caused by protons, nuclei and electrons, including nuclear and electromagnetic
cascades induced by CR nuclei in radio galaxies.

###  2.1 Photodisintegration process

The photodisintegration process is inelastic interaction between CR nuclei and
target photons, resulting in the emission of one or more nucleons or lighter
nuclei from a parent nuclei, as

$A+\gamma\to A_{d}^{\star}+{\rm nucleons},$ (8)

where $A$ represents a primary nuclei, $\gamma$ is an incoming target photon
and $A_{d}^{\star}$ is a daughter nuclei. This process occurs when the
incoming target photons in the nuclear rest frame have energy greater than the
nuclear binding energy, typically $\bar{\varepsilon}_{\gamma}\gtrsim
10\rm~{}MeV$ (e.g., Rachen, 1996). The photodisintegration process is
typically dominated by the giant dipole resonance (GDR), which may be
approximated by the $\delta$-function, $\sigma_{\rm phdis}\approx\sigma_{\rm
GDR}\delta(\bar{\varepsilon}_{\gamma}-\bar{\varepsilon}_{\rm GDR})$, where
$\bar{\varepsilon}_{\rm GDR}$ is the typical photon energy in the nuclear rest
frame, and the width of the GDR process is $\Delta\bar{\varepsilon}_{\rm
GDR}$. The optical depth to the photodisintegration process is estimated by
(e.g., Murase & Beacom, 2010b),

$\displaystyle\tau_{A\gamma}$ $\displaystyle\approx t_{\rm esc}/t_{\rm phdis}$
$\displaystyle\approx\frac{2}{\alpha+1}(c/V)l_{b}\hat{\sigma}_{\rm
dis}\varepsilon_{b}n_{\varepsilon_{b}}\left(\frac{E_{A}}{E_{A,b}^{\rm
dis}}\right)^{\alpha-1}$ $\displaystyle\sim
0.05\frac{2}{\alpha+1}\left(\frac{A}{16}\right)^{1.21}\left(\frac{V}{0.5c}\right)^{-1}\left(\frac{l_{b}}{1\rm~{}kpc}\right)\left(\frac{E_{A}}{E_{A,b}^{\rm
dis}}\right)^{\alpha-1},$ (9)

where $t_{\rm phdis}$ is the photodisintegration interaction time scale,
$\varepsilon_{b}n_{\varepsilon_{b}}=10^{3}\rm~{}cm^{-3}$,
$\varepsilon_{b}=10\rm~{}eV$, $\hat{\sigma}_{\rm dis}=\kappa_{\rm
GDR}\sigma_{\rm GDR}\Delta\bar{\varepsilon}_{\rm GDR}/\bar{\varepsilon}_{\rm
GDR}$ is the effective photo disintegration cross section, $\kappa_{\rm GDR}$
is the inelasticity, $\sigma_{\rm GDR}\approx 1.45\times
10^{-27}A\rm~{}cm^{2}$, $\bar{\varepsilon}_{\rm GDR}\approx
42.65A^{-0.21}\rm~{}MeV$ ($A>4$), $\Delta\bar{\varepsilon}_{\rm GDR}\approx
8\rm~{}MeV$ and $E_{A,b}^{\rm dis}\approx
0.5\delta_{D}m_{A}c^{2}\bar{\varepsilon}_{\rm GDR}/\varepsilon_{b}\simeq
7.2\times
10^{16}~{}\delta_{D}(A/16)^{0.79}(\varepsilon_{b}/10\rm~{}eV)^{-1}\rm~{}eV$ is
the typical energy of the nuclei that interacts with target photons with
energy $\varepsilon_{b}$. For the photodisintegration process with the
ejection of one nucleon, we have $\kappa_{\rm GDR}=1/A$. However, in general,
the mean inelasticity depends on the contributions from all the dominant
channels, see the details in Appendix A of Zhang et al. (2017).

The corresponding effective optical depth is

$\displaystyle f_{\rm dis}$ $\displaystyle\approx t_{\rm esc}/t_{\rm
dis}\sim\tau_{A\gamma}/A$ $\displaystyle\sim 3\times
10^{-3}\frac{2}{\alpha+1}\left(\frac{A}{16}\right)^{0.21}\left(\frac{V}{0.5c}\right)^{-1}\left(\frac{l_{b}}{1\rm~{}kpc}\right)\left(\frac{E_{A}}{E_{A,b}^{\rm
dis}}\right)^{\alpha-1},$ (10)

where $t_{\rm dis}$ is the energy loss time scale.

The remaining daughter nuclei are in an excited state, which undergoes a
subsequent de-excitation process by emitting one or multiple photons,

$A_{d}^{\star}\to A_{d}+{\rm photons},$ (11)

where $\bar{\varepsilon}_{\gamma,\rm deex}$ is the typical photon energy in
the nuclear rest frame. The value of $\bar{\varepsilon}_{\gamma,\rm deex}$
extends from hundreds keV to a few MeV, which depends on the specific de-
excitation channel (Anchordoqui et al., 2007b; Murase & Beacom, 2010b; Kundu &
Gupta, 2014). The observed de-excited $\gamma$-ray have boosted energy
$E_{\gamma,\rm deex}=\delta_{D}\gamma_{A}\bar{E}_{\gamma,\rm deex}\sim
100\delta_{D}(\gamma_{A}/10^{8})(\bar{E}_{\gamma,\rm
deex}/1\rm~{}MeV)\rm~{}TeV$, where $\gamma_{A}$ is the Lorentz factor of the
daughter nuclei $A_{d}^{\star}$. The detailed inclusive cross section for the
production of de-excitation photons can be calculated with numerical code
Talys (Goriely et al., 2008). However, the direct output of the total cross-
section of the photodisintegration process for light and intermediate-mass
nuclei is inconsistent with experimental data (Batista et al., 2016; Tamii et
al., 2022). For this purpose, we directly adopt the same data files used in
CRPropa 3 which had been prepared with Talys 1.8 using adjusted GDR parameters
in order to better match the experimental data (Batista et al., 2016).
However, considering uncertainties and the discrepancies in the
photodisintegration process for the inclusive cross section for photon
production of Talys with other numerical code, e.g., Fluka (Böhlen et al.,
2014), we also adopt the method used in previous works (e.g., Anchordoqui et
al. (2007b); Murase & Beacom (2010b); Kundu & Gupta (2014)). For numerical
calculations, the injection rate of de-excitation $\gamma$-rays from the
photodisintegration process is calculated using Eq. 16 of Anchordoqui et al.
(2007b). We assume the average energy of emitted photons is
$\bar{\varepsilon}_{\gamma,\rm deex}=2\rm~{}MeV$ when measured in the nuclear
rest frame and the multiplicity is $\bar{n}_{\gamma}=3$ (Morejon, 2021). The
effective optical depth is estimated to be (Murase & Beacom, 2010b)

$\displaystyle f_{\rm deex}$ $\displaystyle\approx t_{\rm esc}/t_{\rm
deex}\approx(\kappa_{\rm deex}/\kappa_{\rm GDR})f_{\rm dis}$
$\displaystyle\sim 2\times
10^{-5}\frac{2}{\alpha+1}\left(\frac{A}{16}\right)^{0.21}\left(\frac{V}{0.5c}\right)^{-1}\left(\frac{l_{b}}{1\rm~{}kpc}\right)\left(\frac{E_{A}}{E_{A,b}^{\rm
dis}}\right)^{\alpha-1},$ (12)

where $\kappa_{\rm deex}\approx\bar{n}_{\gamma}\bar{\varepsilon}_{\gamma,\rm
deex}/m_{A}c^{2}$ is the energy carried by the de-excitation $\gamma$-rays.
Furthermore, we ignore the de-excitation photons generated from the excited
fragments via photomeson production.

###  2.2 Photomeson production process

The photomeson production process occurs when the energy of target photons in
the nuclear rest frame exceeds the pion production threshold,
$\varepsilon_{\rm th}\approx m_{\pi}c^{2}(1+m_{\pi}/2m_{N})\sim 140\rm~{}MeV$.
The effective optical depth of the photomeson production process can be
estimated using the $\Delta$-resonance as

$\displaystyle f_{\rm mes}$ $\displaystyle\approx t_{\rm esc}/t_{\rm
mes}\sim\frac{2}{\alpha+1}(c/V)l_{b}\hat{\sigma}_{\rm
mes}\varepsilon_{b}n_{\varepsilon_{b}}\left(\frac{E_{A}}{E_{A,b}^{\rm
mes}}\right)^{\alpha-1}$ $\displaystyle\sim 0.3\times
10^{-3}\frac{2}{\alpha+1}\left(\frac{V}{0.5c}\right)^{-1}\left(\frac{l_{b}}{1\rm~{}kpc}\right)\left(\frac{E_{A}}{E_{A,b}^{\rm
mes}}\right)^{\alpha-1},$ (13)

where $t_{\rm mes}$ is the photomeson production energy loss time scale,
$\hat{\sigma}_{\rm mes}\sim A\hat{\sigma}_{p\gamma}\approx A\kappa_{\rm
mes}\sigma_{\Delta}\Delta\bar{\varepsilon}_{\Delta}/\bar{\varepsilon}_{\Delta}$,
$\kappa_{\rm mes}\sim\kappa_{p\gamma}/A$ is the nuclear inelasticity,
$\kappa_{p\gamma}\sim 0.2$ is the proton inelasticity, $\sigma_{\Delta}\approx
4.4\times 10^{-28}\rm~{}cm^{2}$, $\bar{\varepsilon}_{\Delta}\simeq
0.34\rm~{}GeV$, $\Delta\bar{\varepsilon}_{\Delta}\simeq 0.2\rm~{}GeV$ and
$E_{A,b}^{\rm mes}\approx
0.5\delta_{D}m_{A}c^{2}\bar{\varepsilon}_{\Delta}/\varepsilon_{b}\simeq
2.5\times 10^{17}\delta_{D}(A/16)(\varepsilon_{b}/10\rm~{}eV)^{-1}\rm~{}eV$
(e.g., Murase & Beacom, 2010b). The photomeson production interaction time
scale $t_{\rm phmes}$ can be estimated without considering the inelasticity.
In our numerical approach, we employ the Monte Carlo event generator SOPHIA to
determine the differential cross-section of all stable secondary particles
(Mücke et al., 2000). We adopt the superposition model, where the photomeson
cross section is $\sigma_{\rm mes}\sim A\sigma_{p\gamma}$, where
$\sigma_{p\gamma}$ is the photomeson cross-section of protons (Zhang & Murase,
2019). The superposition model is employed in CRPropa 3 with a slightly
different scaling law as described in Eq. 3 in Kampert et al. (2013), which
assumes the emission of one proton or neutron in each interaction process. In
accordance with Batista et al. (2016), we assume that the cross-section for
nuclei with the mass number $A>8$ is $0.85(Z+N)\sigma_{p\gamma}$ times higher,
where $Z$ and $N$ are the numbers of protons and neutrons, respectively. On
the other hand, for nuclei with the mass number $1<A\leqslant 8$, the cross
section is estimated to be $0.85(Z^{2/3}+N^{2/3})\sigma_{p\gamma}$. In Morejon
et al. (2019), a comprehensive examination of the photomeson production
process that incorporates the impact of the nuclear medium and the
fragmentation of the primary nucleus is presented. In the context of radio
galaxies under consideration in this work, the target photons typically have
lower energies and the photodisintegration process of primary nuclei is
dominant, thus we opt for the simplified superposition model for ease of
calculation. However, the photomeson production process and nuclear
fragmentation become significant when the target photons are produced by the
prompt emission of relativistic jets, such as in $\gamma$-ray bursts and tidal
disruption events (Murase et al., 2008; Murase & Beacom, 2010b; Morejon et
al., 2019).

###  2.3 Bethe-Heitler pair production process

The Behte-Heitler pair production process will occur once the target photons
have energy beyond $\varepsilon_{\rm th}\approx 2m_{e}c^{2}\sim 1\rm~{}MeV$ in
the nuclear rest frame (e.g., Blumenthal, 1970a). The effective optical depth
to the Bethe-Heitler pair production process is estimated to be

$\displaystyle f_{\rm BH}$ $\displaystyle\approx t_{\rm esc}/t_{\rm
BH}\sim\frac{2}{\alpha+1}(c/V)l_{b}\hat{\sigma}_{\rm
BH}\varepsilon_{b}n_{\varepsilon_{b}}\left(\frac{E_{A}}{E_{A,b}^{\rm
BH}}\right)^{\alpha-1}$ $\displaystyle\sim 2\times
10^{-5}\frac{2}{\alpha+1}\left(\frac{Z}{8}\right)^{2}\left(\frac{A}{16}\right)^{-1}\left(\frac{V}{0.5c}\right)^{-1}\left(\frac{l_{b}}{1\rm~{}kpc}\right)\left(\frac{E_{A}}{E_{A,b}^{\rm
BH}}\right)^{\alpha-1},$ (14)

where $\hat{\sigma}_{\rm BH}\sim 8\times 10^{-31}(Z^{2}/A)\rm~{}cm^{2}$ is the
pair production cross section of nuclei, $E_{A,b}^{\rm BH}\approx
0.5\delta_{D}\bar{\varepsilon}_{\rm BH}m_{A}c^{2}/\varepsilon_{b}\simeq
7.5\times 10^{15}\delta_{D}(A/16)(\varepsilon_{b}/10\rm~{}eV)^{-1}\rm~{}eV$
and $\bar{\varepsilon}_{\rm BH}\approx 10\rm~{}MeV$. To calculate the energy
loss rate of the photopair production process for relativistic nuclei, we
employ Equation 3.11 in Chodorowski et al. (1992). The energy spectrum of
secondary electron-positron pairs can be determined by utilizing Eq. 62 in
Kelner & Aharonian (2008). This equation calculates the double-differential
cross-section of emitted electrons (and positrons) as a function of energy and
emission angle in the nuclear rest frame, which is obtained from Eq. 10 in
Blumenthal (1970b). We note that a factor of 2 should be multiplied when using
Eq. 62 in Kelner & Aharonian (2008) to accurately account for both electrons
and positrons.

Note that the relative contribution of de-excitation VHE $\gamma$-rays and the
Bethe-Heitler pair production process to the observed flux at the TeV energy
range depends on several factors, such as the source magnetic field strength,
target photon fields, and the composition of UHECR nuclei, as noted in studies
by Murase & Beacom (2010b); Aharonian & Taylor (2010). The energy loss rate
due to de-excitation $\gamma$-rays can be lower than that of the Bethe-Heitler
pair production process, but the relative contribution depends on details of
electromagnetic cascades from the Bethe-Heitler electron-positron pairs (e.g.,
Murase & Beacom, 2010b).

###  2.4 Electromagnetic cascade

High-energy electrons and positrons will lose energy through processes such as
synchrotron emission and inverse-Compton (IC) scattering, which are influenced
by the strength of the magnetic fields and the density of the ambient target
photon fields, respectively. The synchrotron energy loss time scale for high-
energy particles is

$t_{\rm
syn}^{-1}\approx\frac{4\sigma_{T}Z^{4}m_{e}^{2}}{3m^{4}c^{3}}\frac{E_{A}}{\delta_{D}}U_{B},$
(15)

where $Z$ is the particle charge, $m$ is the particle mass, $\sigma_{T}$ is
the Thomson cross section, and $U_{B}=B^{2}/8\pi$ is the magnetic energy
density (e.g., Rybicki & Lightman, 1986). The IC energy loss time scale for
high-energy electrons is given by

$t_{\rm
IC}^{-1}\approx\frac{4\sigma_{T}c}{3m_{e}^{2}c^{4}}\frac{E_{e}}{\delta_{D}}U_{\rm
ph}F_{\rm KN},$ (16)

where $U_{\rm ph}$ is the energy density of target photons and $F_{\rm KN}$ is
a factor to take into account Klein-Nishina effect (Jones, 1968). In our
numerical calculations, we use the known results of the total and differential
cross sections (e.g., Rybicki & Lightman, 1986). In this work, we also adopt
the continuous energy loss approximation, which is further detailed in
Appendix A.

High-energy $\gamma$-rays may interact with target photons, leading to the
creation of electron-positron pairs. The two-photon annihilation optical depth
can be estimated by

$\displaystyle\tau_{\gamma\gamma}$ $\displaystyle\approx t_{\rm
lc}/t_{\rm\gamma\gamma}\sim\eta_{\gamma\gamma}\sigma_{T}l_{b}\varepsilon_{b}n_{\varepsilon_{b}}\left(\frac{E_{\gamma}}{E_{\gamma,b}}\right)^{\alpha-1}$
$\displaystyle\sim
0.2\left(\frac{l_{b}}{1\rm~{}kpc}\right)\left(\frac{E_{\gamma}}{E_{\gamma,b}}\right)^{\alpha-1},$
(17)

where $t_{\gamma\gamma}$ is the two-photon annihilation time scale,
$\eta_{\gamma\gamma}\sim 0.1$ is a numerical factor that depends on the target
photon spectral index (e.g., Svensson, 1987),
$E_{\gamma,b}\approx\delta_{D}(m_{e}c^{2})^{2}/\varepsilon_{b}\simeq 2.6\times
10^{10}\delta_{D}(\varepsilon_{b}/10)^{-1}\rm~{}eV$. Note that the validity of
Eq. 2.4 is based on the assumption that the target photon spectrum is soft,
with $\alpha\gtrsim 1$ (e.g., Dermer & Menon, 2009). In this work for
simplicity we also assume that electrons and positrons share the photon energy
(but the detailed distribution may be used as in Murase & Beacom, 2012). The
injection of high-energy $\gamma$-rays and electrons will trigger an
electromagnetic cascade process, which is essential for predicting the
observed multi-wavelength energy spectrum.

In addition to the isotropic IC scattering case, we also consider the effect
of anisotropic IC scattering process (e.g., Brunetti, 2000). It had been
proposed that the observed flux from the anisotropic IC scattering process is
sensitive to the observation angle, especially for nearby radio galaxies
(e.g., Brunetti, 2000; Bednarek, 2020).

##  3 Application to the nearest radio galaxy Centaurus A

###  3.1 Physical model

In this section, we apply the above method to the nearest radio galaxy, Cen A.
Motivated by the possible connection between Cen A and the observed UHECRs
(e.g., Biteau et al., 2021), the hadronic origin of the VHE $\gamma$-rays from
Cen A has been explored by various authors based on the inner core model
(e.g., Fraija, 2014; Petropoulou et al., 2014; Fraija et al., 2018; Joshi et
al., 2018; Fraija et al., 2018; Banik et al., 2020). However, due to the
smaller distance of the inner core to the central black hole, these models
contradict the observed morphology of VHE $\gamma$-rays, which is consistent
with the origin from the kiloparsec-scale jet (H. E. S. S. Collaboration et
al., 2020).

It is natural to expect that charged particles, including electrons and
nuclei, could accelerate in the kiloparsec-scale jet via stochastic
acceleration or shear acceleration. The accelerated high-energy electrons
could up-scatter the surrounding target photon fields to the VHE energy range,
such as infrared photons from dust torus (H. E. S. S. Collaboration et al.,
2020; Liu et al., 2017; Wang et al., 2021), optical and ultra-violet emission
from disk starlight (Hardcastle & Croston, 2011; Tanada et al., 2019) and
broad-band non-thermal emission from the inner core or “hidden” core
(Bednarek, 2019, 2020).

Shear acceleration could operate when charged particles are scattered off the
magnetic field inhomogeneities from different layers of the shearing flow
(e.g., Rieger & Duffy, 2019). In the steady-state, the energy spectrum of
accelerated particles typically may follow a power-law distribution with an
exponential cutoff $dN^{A}/d\varepsilon\propto\varepsilon^{-s_{\rm acc}}{\rm
exp}(\varepsilon/\varepsilon_{\rm max})$, where $s_{\rm acc}$ is the spectral
index. The spectral index $s_{\rm acc}$ of accelerated particles without
radiation energy losses could be much steeper for non-relativistic flow speeds
because of the efficient diffusive escape process (Rieger & Duffy, 2022). On
the other hand, the spectral index is harder for trans-relativistic and
relativistic flows and the shear acceleration has been proposed as an
effective mechanism to accelerate charged nuclei to the UHE energy range in
the large-scale jet (e.g., Kimura et al., 2018; Wang et al., 2022; Rieger,
2022; Seo et al., 2023). For large-scale relativistic jets, not only the shear
reacceleration but also the one-shot reacceleration can be efficient,
especially for inner jets and powerful FR-II jets (Mbarek & Caprioli, 2021;
Mbarek et al., 2023). The maximum available energy of the accelerated CR
nuclei can be estimated under the Hillas-type confinement condition (Hillas,
1984),

$\displaystyle E_{\rm max}$ $\displaystyle\approx\frac{1}{\eta}\Gamma_{\rm
kpc}ZeB\beta l_{b}$ $\displaystyle\simeq 110(Z/26){\Gamma_{\rm
kpc}}\left(\frac{\eta}{7}\right)^{-1}\left(\frac{B}{10^{-4}\rm~{}G}\right)\left(\frac{\beta}{0.3}\right)\left(\frac{l_{b}}{1\rm~{}kpc}\right)\rm~{}EeV,$
(18)

where $\eta$ represents a prefactor which is $\sim$a few in the Bohm limit
(e.g., Drury, 1983), $B$ is the magnetic field strength, $\beta$ is the shock
velocity. From Eq. 3.1, we can see the required magnetic luminosity is (Pe’Er
et al., 2009; Murase et al., 2012)

$L_{B}\sim 4\times
10^{44}{(Z/26)}^{-2}\left(\frac{\eta}{7}\right)^{2}\left(\frac{E_{\rm
max}}{100\rm~{}EeV}\right)^{2}\left(\frac{\beta}{0.3}\right)^{-1}\rm~{}erg~{}s^{-1}.$
(19)

The Eq. 19 gives the condition for a source capable of acceleration particles
to 100 EeV. Note that the composition of the accelerated UHECR nuclei could be
different from the typical composition of the interstellar medium (e.g.,
Kimura et al., 2018). In this work, for the demonstration, we only consider
two typical elements of heavier nuclei, oxygen, and iron nuclei, to study the
detectability of de-excitation VHE $\gamma$-rays.

For simplicity, we assume the emission region in the kiloparsec-scale jet is
modeled as a spherical blob moving at sub-relativistic speed towards us with
the same inclination angle $\theta_{\rm ob}$ as the inner core (H. E. S. S.
Collaboration et al., 2020). The viewing angle of the jet is $\theta_{\rm
ob}\sim 20^{\circ}-50^{\circ}$ (Tingay et al., 1998; Hardcastle et al., 2003).
The jet velocity can be derived from observed apparent velocity $\beta_{\rm
app}$ combining with the viewing angle $\theta_{\rm ob}$, $\beta=\beta_{\rm
app}/(\beta_{\rm app}{\rm cos}\theta_{\rm ob}+{\rm sin}\theta_{\rm ob})$.
Similar to Bednarek (2019, 2020), we consider the broad-band non-thermal
emission originating from the inner core as the dominant target photon field
in the kiloparsec-scale jet. The comoving frame photon energy density of the
inner core emission “observed” in the kpc-scale jet is

${\varepsilon}^{2}\frac{dn}{d\varepsilon}=\frac{d_{L}^{2}}{R^{2}c}\frac{\delta_{D,\rm
core,kpc}^{4}}{\delta_{D,\rm core,ob}^{4}}\delta_{D,\rm deboost}^{4}EF_{E},$
(20)

where $d_{L}$ is the luminosity distance to the observer at Earth,
$R=1\rm~{}kpc$ is the distance of the kiloparsec-scale jet to the inner core,
$\delta_{D,\rm core,ob}=1/(\Gamma_{\rm core}(1-\beta_{\rm core}{{\rm
cos}\theta_{\rm ob}}))$ is the Doppler factor of the inner core when observed
at Earth, $\delta_{D,\rm core,kpc}=1/(\Gamma_{\rm core}(1-\beta_{\rm core}))$
is the Doppler factor of the inner core when observed at the distance of kpc-
scale jet in the black hole rest frame, $\delta_{D,\rm deboost}=1/(\Gamma_{\rm
kpc}(1+\beta_{\rm kpc}))$ is the de-boosting factor when converting the inner
core emission from black hole rest frame to the comoving frame, and $EF_{E}$
is the observed SED of the inner core. In this work, we adopt $\Gamma_{\rm
core}=5$ and $\Gamma_{\rm kpc}=1.05$. If we adopt the jet viewing angle as
$\theta_{\rm ob}=40^{\circ}$, the corresponding Doppler factor are
$\delta_{D,\rm core,ob}\simeq 1$, $\delta_{D,\rm core,kpc}\simeq 8$ and
$\delta_{D,\rm deboost}\simeq 0.7$. We fit the observed SED of inner core
emission with the following formula,

$EF_{E}=Ax^{-\alpha_{1}}\left(\frac{1+x^{1/\alpha_{3}}}{2}\right)^{(\alpha_{1}-\alpha_{2})\alpha_{3}},$
(21)

where $x=E/E_{b}$. The adopted values for the low-energy bump are $A=1.5\times
10^{-10}\rm~{}erg~{}cm^{-2}~{}s^{-1}$, $E_{b}=1\rm~{}eV$, $\alpha_{1}=-0.8$,
$\alpha_{2}=3.2$, and $\alpha_{3}=1$. The values for the high-energy bump are
$A=4.5\times 10^{-10}\rm~{}erg~{}cm^{-2}~{}s^{-1}$, $E_{b}=10^{6}\rm~{}eV$,
$\alpha_{1}=-0.6$, $\alpha_{2}=1$, and $\alpha_{3}=3$.

###  3.2 Time Scales

Figure 2: Different time scales of various energy loss processes in the
kiloparsec-scale jet, with the light crossing time scale indicated as a dotted
line, the advection time scale represented by the black dashed line and the
diffusion escape time scales shown in grey lines.

In Fig. 2, we show various time scales in the kiloparsec-scale jet. We model
the kiloparsec-scale jet as a spherical blob located at a distance $r_{\rm
kpc}=1\rm~{}kpc$ from the inner core. The blob has a comoving radius of
$l_{b}=0.5\rm~{}kpc$ and moves with a velocity of $V=\beta_{\rm kpc}c=0.3$
towards the observer with a viewing angle of $\theta_{\rm ob}\sim
40\rm~{}deg$. The Doppler factor is $\delta_{D,\rm kpc}$. To simplify
calculations, we assume isotropic target photon fields are uniformly
distributed throughout the kiloparsec-scale jet region. However, in reality,
the emission from the inner core is highly beamed and concentrated along the
jet axis, varying depending on the activity of the central supermassive black
holes. The mean-field magnetic field strength is set to $B=1\times
10^{-4}\rm~{}G$.

The black-dashed line represents the advection time scale $t_{\rm
adv}=l_{b}/V$ for charged particles. We also consider the diffusive escape
process for oxygen nuclei (grey dashed line) and iron nuclei (grey dot-dashed
line). The diffusive escape time scales can be estimated by $t_{\rm
diff}\approx l_{b}^{2}/6D$ assuming spherical geometry, where $D$ is the
diffusion coefficient. We adopt the following form of the diffusion
coefficient

$D(E)\approx\frac{cl_{\rm
coh}}{3}\left[4\left(\frac{E}{E_{c}}\right)^{2}+a_{I}\left(\frac{E}{E_{c}}\right)+a_{L}\left(\frac{E}{E_{c}}\right)^{1/3}\right],$
(22)

where $a_{I}=0.9$ and $a_{L}=0.23$ for a Kolmogorov spectrum with $m=5/3$
(Harari et al., 2014). The critical energy $E_{c}$ is given by

$E_{c}\approx ZeBl_{\rm coh}\simeq 9.6Z\frac{B}{10^{-4}\rm~{}G}\frac{l_{\rm
coh}}{\rm 0.1~{}kpc}\rm~{}EeV,$ (23)

where $l_{\rm coh}=0.1\rm~{}kpc$ is the coherence length. Note the escape time
scale should be larger than the light crossing time scale in order to avoid
the superluminal escape. For charged particles, we use the effective
confinement time $t_{\rm conf}={\rm max}[t_{\rm diff},t_{\rm lc}]$, while for
neutral particles, the escape time scale is equivalent to $t_{\rm
lc}=l_{b}/c$, where $t_{\rm lc}$ represents the light crossing time scale.

We also show the Bethe-Heitler pair production energy loss rate, photomeson
production energy loss rate, and photodisintegration energy loss rate for
UHECR oxygen and iron nuclei, respectively. The fate of CR nuclei inside the
kiloparsec-scale jet can be characterized by the effective optical depth,
$f_{\rm dis}\approx t_{\rm esc}/t_{\rm dis}$, given that the photonuclear
reaction is dominated by photodisintegration process. As indicated in Fig. 2,
both oxygen and iron nuclei could survive up to the highest energy. We
calculate the de-excitation energy loss rate for UHECR oxygen and iron nuclei,
respectively. The ratio between the effective optical depth of the de-
excitation process and the Bethe-Heitler pair production process above
energies of nuclei interacting with photons at $\varepsilon_{b}$ is (e.g.,
Murase & Beacom, 2010b; Aharonian & Taylor, 2010)

$f_{\rm deex}/f_{\rm BH}\sim
1(0.1)^{\alpha-1}\left(\frac{Z}{8}\right)^{-2}\left(\frac{A}{16}\right)^{0.21\alpha+1},$
(24)

This is consistent with the numerical values displayed in Figure 2 considering
$\alpha\sim 1$, where $f_{\rm deex}/f_{\rm BH}\sim 1$ for oxygen nuclei and
$f_{\rm deex}/f_{\rm BH}\sim 0.4$ for iron nuclei. Our results reveal the
feasibility of detecting de-excitation $\gamma$-rays emitted by light and
intermediate-mass nuclei groups, despite the low-energy loss efficiency.

By examining the blue dotted-dash line, which signifies the two-photon
annihilation time scale, and comparing it to the light crossing time scale, as
shown in Equation 2.4, it is apparent that $\gamma$-rays with high energy
levels of up to approximately 100 TeV may be able to flee from the source.

###  3.3 Results

Figure 3: Predicted multi-wavelength SED from the kiloparsec-scale jet after the injection of primary high-energy electrons (solid blue lines) and UHECR nuclei (solid red lines). The solid black lines represent the inner core’s SED. The upper panel displays the results from the injection of UHECR oxygen nuclei, while the lower panel shows the results from the injection of UHECR iron nuclei. The green squares and red circles represent data points from the inner core obtained from (Abdalla et al., 2018) and cover a range from low-energy radio to optical and high-energy X-ray to GeV bands. The cyan circles represent observed data points from the radio to X-ray band, which correspond to the total flux in the inner region of the large-scale jet and were taken from (Hardcastle et al., 2006). The cyan triangles represent the total X-ray flux, which includes an even larger region of the large-scale jet, as reported by (Kataoka et al., 2006). The cyan ‘butterfly’ represents the observed VHE $\gamma$-ray flux by H.E.S.S. (H. E. S. S. Collaboration et al., 2020). In addition, this figure displays the expected VHE $\gamma$-ray sensitivities of CTA (50 hours, grey solid lines) (The CTA Consortium, 2019), LHAASO (1 yr, grey dashed lines) (Addazi et al., 2022), and SWGO (5 yr, grey dotted lines) (Albert et al., 2019).  Table 1: Physical parameters in the comoving frame used in the leptohadronic model for the kiloparsec-scale jet. Parameter | Value
---|---
Viewing angle ($\theta_{\rm ob}[{\rm deg}]$) | 40
Lorentz factor ($\Gamma_{\rm kpc}$) | 1.05
Radius ($l_{b}$ [kpc]) | 0.5
Magnetic field strength ($B$ [G]) | $1\times 10^{-4}$
Primary electron index | 2.3
Primary electron minimum energy ($\varepsilon_{e,\rm min}$ [eV]) | $2.5\times 10^{7}$
Primary electron maximum energy ($\varepsilon_{e,\rm max}$ [eV]) | $2\times 10^{12}$
Electron injection energy density ($u_{{\rm inj},e}[\rm eV~{}s^{-1}$] | $3$
Oxygen minimum energy ($\varepsilon_{\rm O,min}$ [eV]) | $1\times 10^{17}$
Oxygen maximum energy ($\varepsilon_{\rm O,max}$ [eV]) | $4\times 10^{18}$
Oxygen spectral index ($s_{\rm acc}$) | 2.6
Oxygen injection energy density ($u_{\rm inj,O}[\rm eV~{}cm^{-3}$] | $300$
Iron minimum energy ($\varepsilon_{\rm Fe,min}$ [eV]) | $1\times 10^{17}$
Iron maximum energy ($\varepsilon_{\rm Fe,max}$ [eV]) | $13\times 10^{18}$
Iron spectral index ($s_{\rm acc}$) | 2.6
Iron injection energy density ($u_{\rm inj,Fe}[\rm eV~{}cm^{-3}$] | $100$

In Fig. 3, we show the results from the injection of UHECR oxygen and iron
nuclei, respectively. The relevant physical parameters are summarized in Table
1. We note that both oxygen and iron nuclei require a minimum injection energy
of $10^{17}\rm~{}eV$. In this study, we assume that the energy spectrum of CRs
accelerated in the kiloparsec-scale jet undergoes a break at $\sim
10^{17}\rm~{}eV$, where the spectral index $s_{\rm acc}$ may be harder below
this energy compared to the value of $2.6$ adopted in this work. This change
in the spectral index could be attributed to different acceleration
mechanisms. In this study, we assume that UHECRs are more likely to be
accelerated via the shear acceleration mechanism, while the low-energy CRs are
predominantly contributed by the diffusive shock acceleration or stochastic
acceleration mechanism. Introducing a low-energy cutoff at $\sim
10^{17}\rm~{}eV$ is crucial to ensure that the total CR energy does not exceed
the jet power.

The solid black lines represent the inner core’s SED fitted with Eq. 21, while
the solid blue and red lines represent the results of the injection of primary
electrons and UHECR nuclei, respectively. We can see that there are four
distinct peaks on the predicted multi-wavelength SED from the injection of
UHECR nuclei. These peaks will be thoroughly discussed in the following
paragraphs.

  1. 1.

The first peak observed in the keV energy range is attributed to the
synchrotron emission generated by electron-positron pairs produced through the
Bethe-Heitler pair production process. The peak frequency of the synchrotron
emission can be roughly estimated using the formula $E_{\rm
syn,pk}\approx\delta_{D,\rm kpc}\gamma_{e}^{2}\hbar eB/(m_{e}c)\simeq
1.1\times 10^{2}\delta_{D,\rm
kpc}(B/10^{-4}{\rm~{}G})(\gamma_{e}/10^{7})^{2}\rm~{}eV$ for electrons with
typical energy $E_{e}\sim 5\delta_{D,\rm kpc}\rm~{}TeV$.

  2. 2.

The second peak observed in the GeV energy range is dominated by the
synchrotron emission from electron-positron pairs produced through the
photomeson production process. These electrons have higher energies, leading
to a peak synchrotron frequency estimated to be $E_{\rm syn,pk}\simeq
0.1\delta_{D,\rm kpc}(B/10^{-4}{\rm~{}G})(\gamma_{e}/10^{10})^{2}\rm~{}GeV$,
assuming typical electron energies of around $E_{e}\sim 5\delta_{D,\rm
kpc}\rm~{}PeV$. As seen in Fig. 3, the leptonic component dominates the flux
observed in the GeV band, compared to the hadronic component.

  3. 3.

The third peak, located in the TeV energy range, is dominated by the IC
emission via electron-positron pairs produced from the Bethe-Heitler pair
production process. The typical energy of IC scattered photons can be
estimated to be $E_{\rm IC,pk}\sim\delta_{D,\rm
kpc}\gamma_{e}^{2}\varepsilon\simeq 10\delta_{D,\rm
kpc}(\gamma_{e}/10^{7})^{2}(\varepsilon/0.1{\rm~{}eV})\rm~{}TeV$, where
$\varepsilon$ is the typical energy of the seed photons in the comoving frame.
As seen in the figure, the hadronic component can contribute to, and in some
cases even dominate, the observed VHE $\gamma$-ray flux above $\sim
1~{}\rm~{}TeV$.

  4. 4.

The last peak in the energy range $\sim 100\rm~{}TeV$ is predominantly caused
by the de-excitation $\gamma$-rays produced during the photodisintegration
process of UHECR nuclei. These de-excitation $\gamma$-rays can be estimated to
have an energy of $E_{\gamma,\rm deex}\approx\delta_{D,\rm
kpc}\gamma_{A}\bar{\varepsilon}_{\gamma,\rm deex}\sim 200\delta_{D,\rm
kpc}(\gamma_{A}/10^{8})(\bar{\varepsilon}_{\gamma,\rm
deex}/2\rm~{}MeV)\rm~{}TeV$, where $\gamma_{A}=10^{8}$ is the Lorentz factor
of nuclei. The uncertainty in the flux of de-excitation $\gamma$-rays, as
discussed in Sec. 2, is represented by the shaded area between the two curves
at energies beyond $\sim 10\rm~{}TeV$. The light-shaded region represents the
spectrum without extra-background light (EBL) absorption, while the dark-
shaded region accounts for absorption by EBL (Gilmore et al., 2012). Our
results indicate that despite strong absorption by the EBL, the detection of
de-excitation $\gamma$-rays by current and future ground-based VHE
$\gamma$-ray telescopes, such as Cherenkov Telescope Array (CTA) (The CTA
Consortium, 2019), the Southern Wide-field $\gamma$-ray Observatory (SWGO)
(Albert et al., 2019), and Large High Altitude Air Shower Observatory (LHAASO)
(Addazi et al., 2022), is possible. Note that LHAASO is included just for
comparison because Cen A is not located in its field of view.

##  4 Discussion and implications

###  4.1 Impact of target photon fields

Note that the detectability of the de-excitation $\gamma$-rays is sensitive to
both the target photon energy density and magnetic field strength. The
observed flux in the X-ray band, dominated by the synchrotron emission from
pair-induced electrons, increases with a stronger magnetic field. However,
suppose the magnetic field energy density is much lower than the target photon
energy density. In that case, the IC emission from pair-induced electrons can
surpass the de-excitation $\gamma$-ray flux near its peak energy. Similarly, a
low target photon energy density reduces the energy loss efficiency of the
photodisintegration process. When the target photon field is dense, the escape
of de-excitation $\gamma$-rays becomes impossible due to two-photon
annihilation. According to Eq. 2.1 and Eq. 2.4, the ratio between the optical
depth of the two-photon annihilation process and the photodisintegration
process is

$\tau_{\gamma\gamma}/\tau_{\rm dis}\sim
2\left(\frac{V}{0.5c}\right)2.7^{\alpha-1}\left(\frac{A}{16}\right)^{0.79\alpha-2}\left(\frac{E_{\gamma}}{E_{A}}\right)^{\alpha-1}.$
(25)

We can see $\tau_{\gamma\gamma}/\tau_{\rm dis}\sim 2$ for $\alpha\sim 1$. As
pointed out by Murase et al. (2008), the sources where UHECR nuclei can
survive would be optically thin to high-energy $\gamma$-rays (e.g., Murase et
al., 2008). In conclusion, the parameter space allowing for the detection of
de-excitation $\gamma$-rays is a rather limited, and further study of the
available parameter space is necessary.

In our calculations, we initially assume an isotropic distribution of target
photons in the comoving frame of the emission region. However, we now consider
the impact of the anisotropically distributed beamed target photon fields on
the observed IC spectrum (e.g., Bednarek, 2019, 2020). When viewing an
approaching jet, the IC radiation is reduced compared to the isotropic case
due to the small scattering angle between the line-of-sight and the target
photon beam direction. Conversely, the large scattering angle between the
line-of-sight and the target photon beam direction enhances the IC radiation
from the counter jet. However, the Doppler beaming effect causes the emission
from the counter jet to be reduced by a few factors compared to the emission
from the approaching jet. One significant outcome of the anisotropic IC
scattering effect is that the radiation flux from Bethe-Heitler electron-
positron pairs can be significantly lower than the flux of de-excitation
$\gamma$-rays when viewing an approaching jet, making it easier to identify
de-excitation $\gamma$-rays.

To perform our calculations within the framework of the one-zone model, we
assume that the photons from the inner core are uniformly distributed in the
comoving frame of the kiloparsec-scale jet. However, future research should be
conducted to study the non-uniform distribution of target photons and the
diffusion of UHECR nuclei inside the kiloparsec-scale jet, potentially through
Monte Carlo simulations.

###  4.2 Implications for UHECRs

The injection luminosity of UHECR nuclei can be estimated as (e.g., Dermer et
al., 2012),

$\displaystyle L_{\rm inj}^{\rm UHECR}$ $\displaystyle\approx
2\Omega_{j}\beta_{\rm kpc}cR^{2}\Gamma_{\rm kpc}^{2}u_{\rm inj},$ (26)

where $\Omega_{j}\approx\pi l_{b}^{2}/R^{2}$ is the opening solid-angle and
$u_{\rm inj}$ is the comoving frame UHECR injection energy density. The
injection luminosity of UHECR nuclei is $L_{\rm inj}^{\rm O}\sim 7.2\times
10^{43}\rm~{}erg~{}s^{-1}$ oxygen nuclei component and $L_{\rm inj}^{\rm
Fe}\sim 2.4\times 10^{43}\rm~{}erg~{}s^{-1}$ iron nuclei component,
respectively. The total jet power should be larger than the value estimated in
Eq. 26, considering additional contributions from thermal particles, low-
energy cosmic rays, radiation fields, and magnetic fields. The mean jet power
of Cen A inferred from the observed enthalpy and age of the southern inner
blob is $L_{\rm jet}\sim 10^{43}\rm~{}erg~{}s^{-1}$ (e.g., Croston et al.,
2009; Wykes et al., 2013). It is apparent that the total jet power estimated
in this work is larger than the mean jet power of Cen A. However, higher jet
power is still allowed if the jet activity of Cen A has been intermittent. An
upper limit of the jet power is the Eddington luminosity $L_{\rm Edd}=4\pi
Gm_{p}M_{\bullet}/\sigma_{T}\sim 7\times 10^{45}\rm~{}erg~{}s^{-1}$ with a
black hole mass $M_{\bullet}\sim 5.5\times 10^{7}M_{\odot}$ (e.g., Neumayer et
al., 2007; Cappellari et al., 2009). If we assume that the energy density in
thermal particles and low-energy cosmic rays are not far away from magnetic
energy density, we find the total jet power is still less than the Eddington
luminosity.

Cen A, along with other radio galaxies, has been proposed as a candidate
source of UHECRs detected on Earth (e.g., Kimura et al., 2018; Eichmann et
al., 2018; Matthews et al., 2018; Bell & Matthews, 2022; Taylor et al., 2023).
However, recent studies by Eichmann et al. (2022) indicate that the CR power
of Cen A is about an order of magnitude smaller than its jet power, i.e.,
$\sim 10^{42}\rm~{}erg~{}s^{-1}$, which reduces its likelihood as the main
source of UHECRs. Moreover, this value is about two orders of magnitude
smaller than the amount of CR nuclei required in this work. Nevertheless, if a
heavy composition of UHECR nuclei is considered, Cen A could still contribute
significantly to the observed UHECRs, as heavy nuclei are less likely to
violate the strong quadruple anisotropy constraint (Eichmann et al., 2022). To
account for the intermediate-scale anisotropies observed in UHECRs, it has
been suggested that the scattering of UHECRs emitted from Cen A by the local
structure may be a contributing factor (Bell & Matthews, 2022; Taylor et al.,
2023).

The escaped flux of charged particles from sources depends on the details of
magnetic fields. The confinement time scale of charged particles could be
estimated as $t_{\rm conf}\approx{\rm max}[t_{\rm diff},t_{\rm lc}]$. Assuming
that UHECRs are isotropized in lobes and/or large-scale structures, the
luminosity of escaping cosmic rays are

$L_{\rm esc}^{A}\approx\frac{1}{t_{\rm
conf}}\varepsilon_{A}^{2}\frac{dn^{A}}{d{\varepsilon_{A}}}\Gamma_{\rm
kpc}^{2}2\pi l_{b}^{3},$ (27)

where $dn^{A}/d\varepsilon_{A}$ is the steady-state differential energy
density of charged particles. In Fig. 4, we show the predicted fluxes of CRs
and neutrinos that have reached Earth after escaping from their sources. Our
calculations, assuming the rectilinear propagation of CRs without energy
losses during propagation, indicated that the expected UHECR flux can be below
the observed values, as indicated by black lines. However, the propagation of
UHECRs is strongly impacted by the strength of the intergalactic magnetic
field, leading to the magnetic horizon effect, which limits the arrival of CRs
at the Earth to only the highest energy particles.

Neglecting energy losses during propagation and the effect of Galactic
magnetic fields, the observed flux of CRs at Earth can be expressed by

$E_{A}F_{E_{A}}=\frac{L_{\rm esc}^{A}}{4\pi d_{L}^{2}}\xi(E_{A},d_{L},t_{\rm
act}),$ (28)

where $d_{L}$ is the source luminosity distance and $t_{\rm act}$ is the
activity time. The emission from Cen A is assumed to have been continuous but
with a recent burst of activity, with $t_{\rm act}=1\rm~{}Myr$, where $t_{\rm
act}$ is the burst lifetime. The enhancement factor $\xi(E_{A},d_{L},t_{\rm
act})$ can be estimated as (e.g., Harari et al., 2021; Eichmann et al., 2022)

$\xi(E_{A},d_{L},t_{\rm act})\approx\frac{1}{\mathcal{C}(E_{A},d_{L})}{\rm
exp}\left[-\left(\frac{d_{L}^{2}}{0.6l_{D}(ct_{\rm
act}+d_{L})}\right)^{0.8}\right],$ (29)

where

$\mathcal{C}(E_{A},d_{L})=\frac{l_{D}}{3d_{L}}\left[1-{\rm
exp}\left(-3\left(\frac{d_{L}}{l_{D}}\right)-3.5\left(\frac{d_{L}}{l_{D}}\right)^{2}\right)\right],$
(30)

$l_{D}\equiv 3D/c$ is the diffusion length and D is the diffusion coefficient
(see Eq. 22). Note $ct_{\rm act}+d_{L}$ corresponding to the maximum distance
traveled by the observed CRs, where the distance for rectilinear propagation
is $d_{L}$. Our expected UHECR flux is in agreement with the observed results
when taking into account the average intergalactic magnetic field strength of
$B=10^{-9}\rm~{}G$ and the typical coherence length of $l_{\rm
coh}=100\rm~{}kpc$ (e.g., Harari et al., 2021). Note that the burst lifetime
is much smaller than the typical source lifetime of Cen A, $t_{\rm act}\ll
t_{\rm source}\sim 100\rm~{}Myr$ (e.g., Taylor et al., 2023). The typical time
delay $t_{\rm delay}$ between the arrival time of UHECRs and photons that are
emitted from sources simultaneously is (e.g., Miralda-Escudé & Waxman, 1996;
Murase & Takami, 2009)

$t_{\rm delay}(E_{A})\sim
10\left(\frac{E_{A}/Z}{10^{18}\rm~{}eV}\right)^{-2}\left(\frac{d_{L}}{3\rm~{}Mpc}\right)^{2}\left(\frac{B}{10^{-8}\rm~{}G}\right)^{2}\left(\frac{l_{\rm
coh}}{0.1\rm~{}Mpc}\right)\rm~{}Myr.$ (31)

The time-profile spread is comparable to time delay, $\sigma_{d}(E_{A})\sim
t_{\rm delay}(E_{A})$ (Takami & Murase, 2012).

###  4.3 Implications for Neutrinos

In Fig. 4, we also show the predicted all-flavor energy spectrum of high-
energy neutrinos from pion decay and neutron $\beta$-decay during the
photodisintegration process. Although we adopt the numerical approach, the
observed all-flavor energy spectrum of high-energy neutrinos produced from the
pion decay process can be analytically estimated with the following formula,

$\displaystyle E_{\nu}F_{E_{\nu}}^{(\rm mes)}$
$\displaystyle\approx\delta_{D,\rm kpc}^{4}\frac{1}{4\pi
d_{L}^{2}}\frac{3}{8}f_{\rm
mes}(\varepsilon_{A})\varepsilon_{A}L^{\prime}_{\varepsilon_{A}},$ (32)

where $\varepsilon_{A}L^{\prime}_{\varepsilon_{A}}$ is the comoving frame
cosmic-ray injection luminosity, $E_{\nu}\approx\delta_{D,\rm
kpc}0.05(\varepsilon_{A}/A)\sim\delta_{D,\rm kpc}10^{16}\rm~{}eV$ is the
typical neutrino energy in the observer frame for $\varepsilon_{A}\sim 3\times
10^{18}\rm~{}eV$ and $A=16$. Note both direct neutrino contribution by the
photomeson production process on nuclei and indirect neutrino contribution by
the photomeson production process on secondary neutrons and protons are
included (Murase & Beacom, 2010a; Zhang & Murase, 2019). The observed energy
spectrum of anti-electron neutrinos from neutron $\beta$-decay can be written
as

$\displaystyle E_{\nu}F_{E_{\nu}}^{(\beta_{\rm dec})}$
$\displaystyle\approx\delta_{D,\rm kpc}^{4}\frac{1}{4\pi
d_{L}^{2}}\kappa_{\beta_{\rm dec}}\frac{t_{\rm
esc}^{n}}{\gamma_{n}\tau_{n}}\xi_{n}f_{\rm
dis}(\varepsilon_{A})\varepsilon_{A}L^{\prime}_{\varepsilon_{A}},$ (33)

where $\gamma_{n}$ is the Lorentz factor of neutrons, $\tau_{n}\simeq
879.6\rm~{}s$ is the mean lifetime of free neutrons (Workman et al., 2022),
$\kappa_{\beta_{\rm dec}}\approx\langle\varepsilon_{\nu}\rangle/m_{n}c^{2}$ is
the inelasticity, $m_{n}$ is the neutron mass in the neutron rest frame,
$\xi_{n}\sim 1/2$ is the fraction of neutrons in the emitted nucleons. The
average electron kinetic energy in the neutron rest frame is
$\langle\varepsilon_{e}\rangle-m_{e}c^{2}\approx 0.30\rm~{}MeV$. The typical
neutrino energy in the neutron rest frame can be estimated as
$\langle\varepsilon_{\nu}\rangle\approx
Q_{\beta}-(\langle\varepsilon_{e}\rangle-
m_{e}c^{2})-(\langle\varepsilon_{p}\rangle-m_{p}c^{2})\simeq 0.48\rm~{}MeV$,
when measured in the neutron rest frame, where $Q_{\beta}\approx
0.78\rm~{}MeV$ is the $Q$-value representing the difference between the
initial and final mass energies, $\langle\varepsilon_{p}\rangle-m_{p}c^{2}\sim
0.3\rm~{}keV$, is the proton recoil kinetic energy. The typical neutrino
energy in the observer frame is $E_{\nu}\approx\delta_{D,\rm
kpc}\langle\varepsilon_{\nu}\rangle\gamma_{n}\sim\delta_{D,\rm kpc}3\times
10^{13}\rm~{}eV$ for $\varepsilon_{A}\sim 1\times 10^{18}\rm~{}eV$ and $A=16$.
Note that $\beta$-decay from nuclei is not included in this study.

We can see the neutrino flux lies in the $\sim 10\rm~{}PeV$ energy range
dominated by the pion decay process, while the neutron $\beta$-decay process
dominates the neutrino flux in the lower energy range, $\sim 0.01\rm~{}PeV$,
High-energy neutrino emission from Cen A can be searched for by the next-
generation neutrino telescopes, such as KM3Net (Adrian-Martinez et al., 2016)
and IceCube-Gen2 (Aartsen et al., 2021). As Cen A located in the Southern sky,
KM3Net is more sensitive than IceCube-Gen2 for detecting high-energy neutrinos
from Cen A. However, the flux of high-energy neutrinos predicted in this work
is $\sim 1-2$ orders smaller than the prediction of the magnetically-powered
corona model for Cen A (Kheirandish et al., 2021), and the detection of high-
energy neutrinos seem challenging.

Figure 4: Predicted fluxes of high-energy $\gamma$-rays, neutrinos, and cosmic
rays at Earth from Cen A from the kiloparsec-scale jet after the injection of
UHECR oxygen nuclei (solid lines) and iron nuclei (dashed lines). For high-
energy $\gamma$-rays, we also show the unattenuated spectrum without
considering EBL absorption (thin lines).

The grey triangles are the diffuse flux of UHECRs measured by the Pierre Auger
Collaboration (Aab et al., 2020).

##  5 Summary

In this study, we evaluated the feasibility of detecting de-excitation VHE
$\gamma$-rays through our numerical code that self-consistently considers both
nuclear and electromagnetic cascades. The accuracy of the code has been
validated through comparisons with Monte Carlo simulations using a modified
version of CRPropa 3. Our code can be used to explore the behavior of UHECR
nuclei in astrophysical sources under the assumption of a one-zone model.

We then applied our numerical code to the closest radio galaxy, Cen A, which
we considered a potential accelerator of UHECR nuclei in its kiloparsec-scale
jet. In our model, the primary target photons are the beamed inner core
emission that illuminates the jet when viewed along its axis, but this
emission is greatly reduced when viewed from Earth due to the Doppler beaming
effect.

Our results, assuming the dominant injection of UHECR nuclei consisting of
oxygen and/or iron, indicate that the de-excitation VHE $\gamma$-rays will be
the dominant contributor to the multi-wavelength spectrum at $\gtrsim
10-100\rm~{}TeV$ if the UHECR nuclei are dominated by oxygen-group components.
The de-excitation VHE $\gamma$-rays from Cen A could be detected by current
and future ground-based VHE $\gamma$-ray telescopes. The results obtained in
this work provide valuable insight into the composition of UHECR nuclei in
nearby extragalactic sources.

## Acknowledgements

K.M. acknowledges John Beacom, Charles Dermer, and Asaf Pe’er for early
discussions in 2011-2012. The work of B.T.Z. is supported by KAKENHI No.
20H01901. The work of K.M. is supported by the NSF Grant No. AST-1908689, No.
AST-2108466 and No. AST-2108467, and KAKENHI No. 20H01901 and No. 20H05852.

## Data Availability

The data developed for the calculation in this work is available upon request.
The code used will be made public in the future as a part of the AMES.

## References

  * Aab et al. (2020) Aab A., et al., 2020, Physical Review D, 102, 062005
  * Aartsen et al. (2021) Aartsen M. G., et al., 2021, Journal of Physics G: Nuclear and Particle Physics, 48, 060501
  * Abdalla et al. (2018) Abdalla H., et al., 2018, Astron. Astrophys., 619, A71
  * Addazi et al. (2022) Addazi A., et al., 2022, Chin. Phys. C, 46, 035001
  * Adrian-Martinez et al. (2016) Adrian-Martinez S., et al., 2016, J. Phys. G, 43, 084001
  * Aharonian & Taylor (2010) Aharonian F., Taylor A. M., 2010, Astroparticle Physics, 34, 258
  * Aharonian et al. (2009) Aharonian F., et al., 2009, The Astrophysical Journal, 695, L40
  * Albert et al. (2019) Albert A., et al., 2019, arXiv e-prints, p. arXiv:1902.08429
  * Alves Batista et al. (2019) Alves Batista R., et al., 2019, Frontiers in Astronomy and Space Sciences, 6, 23
  * Anchordoqui (2019) Anchordoqui L. A., 2019, Physics Reports, 801, 1
  * Anchordoqui et al. (2007a) Anchordoqui L. A., Beacom J. F., Goldberg H., Palomares-Ruiz S., Weiler T. J., 2007a, Physical Review D, 75, 063001
  * Anchordoqui et al. (2007b) Anchordoqui L. A., Beacom J. F., Goldberg H., Palomares-Ruiz S., Weiler T. J., 2007b, Physical Review Letters, 98, 121101
  * Banik et al. (2020) Banik P., Bhadra A., Bhattacharyya A., 2020, Monthly Notices of the Royal Astronomical Society, 500, 1087
  * Batista et al. (2016) Batista R. A., et al., 2016, Journal of Cosmology and Astroparticle Physics, 2016, 038
  * Bednarek (2019) Bednarek W., 2019, Monthly Notices of the Royal Astronomical Society, 483, 1003
  * Bednarek (2020) Bednarek W., 2020, The Astrophysical Journal, 891, 145
  * Bell & Matthews (2022) Bell A. R., Matthews J. H., 2022, Monthly Notices of the Royal Astronomical Society, 511, 448
  * Biehl et al. (2018) Biehl D., Boncioli D., Fedynitch A., Winter W., 2018, Astronomy & Astrophysics, 611, A101
  * Biteau et al. (2021) Biteau J., et al., 2021, in Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021). Sissa Medialab, Berlin, Germany - Online, p. 307, doi:10.22323/1.395.0307, https://pos.sissa.it/395/307
  * Blumenthal (1970a) Blumenthal G. R., 1970a, Physical Review D, 1, 1596
  * Blumenthal (1970b) Blumenthal G. R., 1970b, Physical Review D, 1, 1596
  * Boncioli et al. (2019) Boncioli D., Biehl D., Winter W., 2019, The Astrophysical Journal, 872, 110
  * Brunetti (2000) Brunetti G., 2000, Astroparticle Physics, 13, 107
  * Böhlen et al. (2014) Böhlen T., et al., 2014, Nuclear Data Sheets, 120, 211
  * Cappellari et al. (2009) Cappellari M., Neumayer N., Reunanen J., van der Werf P. P., de Zeeuw P. T., Rix H.-W., 2009, Mon. Not. Roy. Astron. Soc., 394, 660
  * Chodorowski et al. (1992) Chodorowski M. J., Zdziarski A. A., Sikora M., 1992, The Astrophysical Journal, 400, 181
  * Croston et al. (2009) Croston J. H., et al., 2009, Monthly Notices of the Royal Astronomical Society, 395, 1999
  * Dermer & Menon (2009) Dermer C. D., Menon G., 2009, High Energy Radiation from Black Holes: Gamma Rays, Cosmic Rays, and Neutrinos. https://ui.adsabs.harvard.edu/abs/2009herb.book.....D
  * Dermer et al. (2012) Dermer C. D., Murase K., Takami H., 2012, Astrophys. J., 755, 147
  * Drury (1983) Drury L. O., 1983, Reports on Progress in Physics, 46, 973
  * Eichmann et al. (2018) Eichmann B., Rachen J., Merten L., Vliet A. v., Tjus J. B., 2018, Journal of Cosmology and Astroparticle Physics, 2018, 036
  * Eichmann et al. (2022) Eichmann B., Kachelrieß M., Oikonomou F., 2022, JCAP, 07, 006
  * Fraija (2014) Fraija N., 2014, Monthly Notices of the Royal Astronomical Society, 441, 1209
  * Fraija & Marinelli (2016) Fraija N., Marinelli A., 2016, The Astrophysical Journal, 830, 81
  * Fraija et al. (2018) Fraija N., Aguilar-Ruiz E., Galván-Gámez A., Marinelli A., de Diego J. A., 2018, Monthly Notices of the Royal Astronomical Society, 481, 4461
  * Gilmore et al. (2012) Gilmore R. C., Somerville R. S., Primack J. R., Domínguez A., 2012, Monthly Notices of the Royal Astronomical Society, 422, 3189
  * Goriely et al. (2008) Goriely S., Hilaire S., Koning A. J., 2008, Astronomy & Astrophysics, 487, 767
  * Guido et al. (2021) Guido E., et al., 2021, in Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021). Sissa Medialab, Berlin, Germany - Online, p. 311, doi:10.22323/1.395.0311, https://pos.sissa.it/395/311
  * H. E. S. S. Collaboration et al. (2020) H. E. S. S. Collaboration et al., 2020, Nature, 582, 356
  * Harari et al. (2014) Harari D., Mollerach S., Roulet E., 2014, Physical Review D, 89, 123001
  * Harari et al. (2021) Harari D., Mollerach S., Roulet E., 2021, Physical Review D, 103, 023012
  * Hardcastle & Croston (2011) Hardcastle M. J., Croston J. H., 2011, Monthly Notices of the Royal Astronomical Society, 415, 133
  * Hardcastle et al. (2003) Hardcastle M. J., Worrall D. M., Kraft R. P., Forman W. R., Jones C., Murray S. S., 2003, The Astrophysical Journal, 593, 169
  * Hardcastle et al. (2006) Hardcastle M. J., Kraft R. P., Worrall D. M., 2006, Monthly Notices of the Royal Astronomical Society: Letters, 368, L15
  * Hillas (1984) Hillas A. M., 1984, Annual Review of Astronomy and Astrophysics, 22, 425
  * Hinton & Hofmann (2009) Hinton J. A., Hofmann W., 2009, Ann. Rev. Astron. Astrophys., 47, 523
  * Jones (1968) Jones F. C., 1968, Physical Review, 167, 1159
  * Joshi et al. (2018) Joshi J. C., Miranda L. S., Razzaque S., Yang L., 2018, Monthly Notices of the Royal Astronomical Society: Letters, 478, L1
  * Kalashev & Kido (2015) Kalashev O. E., Kido E., 2015, Journal of Experimental and Theoretical Physics, 120, 790
  * Kampert et al. (2013) Kampert K.-H., Kulbartz J., Maccione L., Nierstenhoefer N., Schiffer P., Sigl G., van Vliet A. R., 2013, Astroparticle Physics, 42, 41
  * Karakula et al. (1994) Karakula S., Kociolek G., Moskalenko I. V., Tkaczyk W., 1994, The Astrophysical Journal Supplement Series, 92, 481
  * Kataoka et al. (2006) Kataoka J., Stawarz L., Aharonian F., Takahara F., Ostrowski M., Edwards P. G., 2006, The Astrophysical Journal, 641, 158
  * Kelner & Aharonian (2008) Kelner S. R., Aharonian F. A., 2008, Physical Review D, 78, 034013
  * Kheirandish et al. (2021) Kheirandish A., Murase K., Kimura S. S., 2021, The Astrophysical Journal, 922, 45
  * Kimura et al. (2018) Kimura S. S., Murase K., Zhang B. T., 2018, Physical Review D, 97, 023026
  * Kotera & Olinto (2011) Kotera K., Olinto A. V., 2011, Annual Review of Astronomy and Astrophysics, 49, 119
  * Kundu & Gupta (2014) Kundu E., Gupta N., 2014, Journal of Cosmology and Astroparticle Physics, 2014, 030
  * Lee (1998) Lee S., 1998, Physical Review D, 58, 043004
  * Liu et al. (2017) Liu R.-Y., Rieger F. M., Aharonian F. A., 2017, The Astrophysical Journal, 842, 39
  * Matthews et al. (2018) Matthews J. H., Bell A. R., Blundell K. M., Araudo A. T., 2018, Monthly Notices of the Royal Astronomical Society: Letters, 479, L76
  * Mbarek & Caprioli (2021) Mbarek R., Caprioli D., 2021, Astrophys. J., 921, 85
  * Mbarek et al. (2023) Mbarek R., Caprioli D., Murase K., 2023, The Astrophysical Journal, 942, 37
  * Miralda-Escudé & Waxman (1996) Miralda-Escudé J., Waxman E., 1996, The Astrophysical Journal, 462
  * Morejon (2021) Morejon L., 2021, PhD thesis, Humboldt U., Berlin, Humboldt U., Berlin, doi:10.18452/22585
  * Morejon et al. (2019) Morejon L., Fedynitch A., Boncioli D., Biehl D., Winter W., 2019, Journal of Cosmology and Astroparticle Physics, 2019, 007
  * Morejon et al. (2021) Morejon L. R., Rodrigues X., Rudolph A., Gao S., Winter W., 2021, in Proceedings of 37th International Cosmic Ray Conference - PoS(ICRC2021). Sissa Medialab, https://video.desy.de/video/Excited-isomer-photons-and-the-VHE-emission-from-Centaurus-A/839c870c60c47465926b0a672f91c3b7
  * Mücke & Protheroe (2001) Mücke A., Protheroe R., 2001, Astroparticle Physics, 15, 121
  * Mücke et al. (2000) Mücke A., Engel R., Rachen J., Protheroe R., Stanev T., 2000, Computer Physics Communications, 124, 290
  * Murase (2009) Murase K., 2009, Phys. Rev. Lett., 103, 081102
  * Murase (2012) Murase K., 2012, ApJ, 745, L16
  * Murase (2018) Murase K., 2018, Physical Review D, 97, 081301
  * Murase (2022) Murase K., 2022, Astrophys. J. Lett., 941, L17
  * Murase & Beacom (2010a) Murase K., Beacom J. F., 2010a, Physical Review D, 81, 123001
  * Murase & Beacom (2010b) Murase K., Beacom J. F., 2010b, Physical Review D, 82, 043008
  * Murase & Beacom (2012) Murase K., Beacom J. F., 2012, J. Cosmology Astropart. Phys., 2012, 043
  * Murase & Takami (2009) Murase K., Takami H., 2009, Astrophys. J. Lett., 690, L14
  * Murase et al. (2008) Murase K., Ioka K., Nagataki S., Nakamura T., 2008, Phys. Rev. D, 78, 023005
  * Murase et al. (2012) Murase K., Dermer C. D., Takami H., Migliori G., 2012, Astrophys. J., 749, 63
  * Murase et al. (2015) Murase K., Kashiyama K., Kiuchi K., Bartos I., 2015, Astrophys. J., 805, 82
  * Neumayer et al. (2007) Neumayer N., Cappellari M., Reunanen J., Rix H., van der Werf P. P., de Zeeuw P. T., Davies R. I., 2007, Astrophys. J., 671, 1329
  * Pe’Er et al. (2009) Pe’Er A., Murase K., Mészáros P., 2009, Phys. Rev. D, 80, 123018
  * Petropoulou et al. (2014) Petropoulou M., Lefa E., Dimitrakoudis S., Mastichiadis A., 2014, Astronomy & Astrophysics, 562, A12
  * Rachen (1996) Rachen J. P., 1996, PhD thesis, Max-Planck-Institute for Radioastronomy, Bonn
  * Rieger (2022) Rieger F. M., 2022, Universe, 8, 607
  * Rieger & Duffy (2019) Rieger F. M., Duffy P., 2019, The Astrophysical Journal, 886, L26
  * Rieger & Duffy (2022) Rieger F. M., Duffy P., 2022, The Astrophysical Journal, 933, 149
  * Rieger & Levinson (2018) Rieger F., Levinson A., 2018, Galaxies, 6, 116
  * Rodrigues et al. (2018) Rodrigues X., Fedynitch A., Gao S., Boncioli D., Winter W., 2018, The Astrophysical Journal, 854, 54
  * Rulten (2022) Rulten C., 2022, Galaxies, 10, 61
  * Rybicki & Lightman (1986) Rybicki G. B., Lightman A. P., 1986, Radiative Processes in Astrophysics. https://ui.adsabs.harvard.edu/abs/1986rpa..book.....R
  * Seo et al. (2023) Seo J., Ryu D., Kang H., 2023, The Astrophysical Journal, 944, 199
  * Sitarek (2022) Sitarek J., 2022, Galaxies, 10, 21
  * Stecker (1969) Stecker F. W., 1969, Physical Review, 180, 1264
  * Sudoh et al. (2020) Sudoh T., Khangulyan D., Inoue Y., 2020, The Astrophysical Journal, 901, L27
  * Svensson (1987) Svensson R., 1987, Monthly Notices of the Royal Astronomical Society, 227, 403
  * Takami & Murase (2012) Takami H., Murase K., 2012, The Astrophysical Journal, 748, 9
  * Tamii et al. (2022) Tamii A., et al., 2022, PANDORA project: photo-nuclear reactions below $A=60$, doi:10.48550/arXiv.2211.03986, http://arxiv.org/abs/2211.03986
  * Tanada et al. (2019) Tanada K., Kataoka J., Inoue Y., 2019, The Astrophysical Journal, 878, 139
  * Taylor et al. (2023) Taylor A. M., Matthews J. H., Bell A. R., 2023, UHECR Echoes from the Council of Giants (arXiv:2302.06489), doi:10.48550/arXiv.2302.06489
  * The CTA Consortium (2019) The CTA Consortium 2019, Science with the Cherenkov Telescope Array. WORLD SCIENTIFIC, doi:10.1142/10986
  * The Fermi-LAT Collaboration et al. (2010) The Fermi-LAT Collaboration et al., 2010, Science, 328, 725
  * Tingay et al. (1998) Tingay S. J., et al., 1998, The Astronomical Journal, 115, 960
  * Urry & Padovani (1995) Urry C. M., Padovani P., 1995, Publications of the Astronomical Society of the Pacific, 107, 803
  * Wang et al. (2021) Wang J.-S., Reville B., Liu R.-Y., Rieger F. M., Aharonian F. A., 2021, Monthly Notices of the Royal Astronomical Society, 505, 1334
  * Wang et al. (2022) Wang J.-S., Reville B., Mizuno Y., Rieger F. M., Aharonian F. A., 2022, Monthly Notices of the Royal Astronomical Society, 519, 1872
  * Workman et al. (2022) Workman R. L., et al., 2022, PTEP, 2022, 083C01
  * Wykes et al. (2013) Wykes S., et al., 2013, Astronomy & Astrophysics, 558, A19
  * Zhang & Murase (2019) Zhang B. T., Murase K., 2019, Physical Review D, 100, 103004
  * Zhang et al. (2017) Zhang B. T., Murase K., Oikonomou F., Li Z., 2017, Physical Review D, 96, 063007
  * Zhang et al. (2018) Zhang B. T., Murase K., Kimura S. S., Horiuchi S., Mészáros P., 2018, Physical Review D, 97, 083010
  * Zhang et al. (2020) Zhang B. T., Petropoulou M., Murase K., Oikonomou F., 2020, The Astrophysical Journal, 889, 118

##  Appendix A Numerical method and related physical processes

The coupled transport equation given by Eq. 2 can be discretized as

$\displaystyle\frac{d}{dt}n_{a,i}$
$\displaystyle=-n_{a,i}\mathcal{A}_{a,i}+\sum_{j\geqslant
i}n_{a,j}\mathcal{B}_{a,j\to i}$ $\displaystyle+\sum_{b}\sum_{j\geqslant
i}n_{b,j}\mathcal{C}_{b,j\to a,i}+\dot{n}_{a,i}^{\rm inj},$ (34)

where the indices $i$ and $j$ represent different energy bins. The above
equations can be numerically solved with the first-order implicit scheme
(e.g., Lee, 1998),

$\displaystyle\frac{n_{a,i}^{m+1}-n_{a,i}^{m}}{\Delta t}$
$\displaystyle=-n_{a,i}^{m+1}\mathcal{A}_{a,i}+\sum_{j\geqslant
i}n_{a,j}^{m+1}\mathcal{B}_{a,j\to i}$ $\displaystyle+\sum_{b}\sum_{j\geqslant
i}n_{b,j}^{m+1}\mathcal{}_{b,j\to a,i}+\dot{n}_{a,i}^{{\rm inj},m},$ (35)

and we have

$\displaystyle n_{a,i}^{m+1}=\frac{\frac{n_{a,i}^{m}}{\Delta
t}+\sum_{j>i}n_{a,j}^{m+1}\mathcal{B}_{a,j\to i}+\sum_{b}\sum_{j\geqslant
i}n_{b,j}^{m+1}\mathcal{C}_{b,j\to a,i}+\dot{n}_{a,i}^{{\rm
inj},m}}{\frac{1}{\Delta t}+\mathcal{A}_{a,i}-\mathcal{B}_{a,i\to
i}},\,\,\,\,\,\,\,\,$ (36)

where $\Delta t$ is the time step, the index $m$ represents the current
particle number density at time $t$, and the index $m+1$ represents the
particle number density at time $t+\Delta t$. In order to keep the accuracy,
we adopt the method used in Murase (2009, 2012); Murase & Beacom (2012), where
at each time step the solution of Eq. 36 is found when the number density of
each species converges.

To check the availability of our transport code, which is a part of the
Astrophysical Multimessenger Emission Simulator (AMES), we compared our
results with Monte Carlo simulations with a modified version of CRPropa 3. For
the Monte Carlo simulations, instead of implementing the injection term
$\dot{n}_{a}^{\rm inj}$, we choose $\dot{n}_{a}^{\rm inj}=0$ when solving Eq.
2 just for the comparison purpose. We also neglect the escape term. We
consider a spherical blob with a radius of $l_{b}=1\rm~{}pc$ moving toward the
observer with a Lorentz factor of $\Gamma=10$. The target photon fields in the
comoving frame of the blob can be described by Eq. 7, which are isotropically
distributed inside the blob with $\varepsilon_{b}=1\rm~{}eV$, $\alpha_{l}=1$,
$\alpha_{h}=2.5$ and $u_{\rm ph}=10^{7}\rm~{}eV$ is the comoving frame target
photon energy density. The injected energy spectrum of oxygen nuclei follows a
power-law distribution with an exponential cutoff, where $\varepsilon_{\rm
max}=8\times 10^{\rm 17}\rm~{}eV$ and $s_{\rm acc}=2$. The dynamical time
scale of the system is given by $t_{\rm dyn}=l_{b}/V$, where $V$ is the
characteristic velocity (e.g., shock crossing time during which photons are
generated).

Figure 5: Spectral energy distributions of primary nuclei and secondary
particles, including daughter nuclei and other particles, which are calculated
with our kinetic code (solid curves) and Monte Carlo code (dots), where we
consider the Bethe-Heitler pair production process (upper panel), photomeson
production process (middle panel), and photodisintegration process (lower
panel), respectively.

In Fig. 5, we compare the output from our kinetic code with what derived from
Monte Carlo simulations using a modified version of CRPropa 3 (Batista et al.,
2016; Zhang et al., 2020). We compare the spectrum of surviving nuclei and
generated photons, electrons, and neutrinos, considering three main hadronic
processes including the Bethe-Heitler pair production process (upper panel),
the photomeson production process (middle panel), and the photodisintegration
process (lower panel). The results from the two calculation methods are
consistent with each other, except for stochastic fluctuations. Indeed, the
difference in the high-energy part of Bethe-Heilter pairs is due to the
limited sampling of target photons for CRPropa 3 in handling the Bethe-Heilter
process.

The details of the value of the coefficients $\mathcal{A}$, $\mathcal{B}$, and
$\mathcal{C}$ for various particles are discussed below. Note all the
quantities are in the comoving frame.

###  A.1 Photon

The kinetic equation for $\gamma$-rays is (see also Supplemental Material of
Murase, 2018)

$\displaystyle\frac{\partial n_{\varepsilon_{\gamma}}^{\gamma}}{\partial t}$
$\displaystyle=-n_{\varepsilon_{\gamma}}^{\gamma}\left(\frac{1}{t_{\rm
esc}^{\gamma}}+\frac{1}{t_{\gamma\gamma}}\right)$
$\displaystyle+\frac{\partial}{\partial t}\left(n_{\varepsilon_{\gamma}}^{\rm
syn}+n_{\varepsilon_{\gamma}}^{\rm IC}+n_{\varepsilon_{\gamma}}^{\rm
phmes}+n_{\varepsilon_{\gamma}}^{\rm
deex}\right)+\dot{n}_{\varepsilon_{\gamma}}^{\rm inj},$ (37)

where $t_{\rm esc}^{\gamma}$ is the photon escape time scale,
$t_{\rm\gamma\gamma}$ is the interaction time scale of high-energy
$\gamma$-rays with ambient target photon fields, $\partial
n_{\varepsilon_{\gamma}}^{\rm syn}/\partial t$ is the $\gamma$-ray generation
rate from synchrotron emission process of all the charged particles, $\partial
n_{\varepsilon_{\gamma}}^{\rm IC}/\partial t$ is the $\gamma$-ray generation
rate from the inverse-Compton process of electrons, $\partial
n_{\varepsilon_{\gamma}}^{\rm phmes}/\partial t$ is the $\gamma$-ray
generation rate from the photomeson production process of both nucleons and
nuclei, $\partial n_{\varepsilon_{\gamma}}^{\rm phdis}/\partial t$ is the
$\gamma$-ray generation rate from the photodisintegration process via de-
excitation process and $\dot{n}_{\varepsilon_{\gamma}}^{\rm inj}$ is the
primary $\gamma$-ray injection rate (e.g., from external regions).

Based on Eq. A.1, the coefficient $\mathcal{A}$ in Eq. 2 corresponds to

$\mathcal{A}_{\gamma}=\frac{1}{t_{\rm
esc}^{\gamma}}+\frac{1}{t_{\rm\gamma\gamma}}.$ (38)

The coefficient $\mathcal{B}$ is set to be zero

$\mathcal{B}_{\gamma\to\gamma}=0,$ (39)

which means that incident high-energy $\gamma$-rays disappear once they
annihilate with target photons. The coefficient $\mathcal{C}$ is

$\mathcal{C}_{\gamma}=\frac{\partial}{\partial
t}\left(n_{\varepsilon_{\gamma}}^{\rm syn}+n_{\varepsilon_{\gamma}}^{\rm
IC}+n_{\varepsilon_{\gamma}}^{\rm phmes}+n_{\varepsilon_{\gamma}}^{\rm
phdis}\right).$ (40)

The escape time scale for photons is set to $t_{\rm esc}^{\gamma}=t_{\rm
lc}=R/c$, where $t_{\rm lc}$ is the light crossing time scale and $R$ is the
radius of the emission region.

###  A.2 Electron

The kinetic equation for high-energy electrons (and positrons) is

$\displaystyle\frac{\partial n_{\varepsilon_{e}}^{e}}{\partial t}$
$\displaystyle=-\frac{n_{\varepsilon_{e}}^{e}}{t_{\rm
esc}^{e}}-\frac{\partial}{\partial\varepsilon_{e}}\left[(P_{\rm ad}+P_{\rm
syn}^{e}+P_{\rm IC}^{e})n_{\varepsilon_{e}}^{e}\right]$
$\displaystyle+\frac{\partial}{\partial
t}\left(n_{\varepsilon_{e}}^{\gamma\gamma}+n_{\varepsilon_{e}}^{\rm
BH}+n_{\varepsilon_{e}}^{\rm phmes}+n_{\varepsilon_{e}}^{\beta_{\rm
dec}}\right)+\dot{n}_{\varepsilon_{e}}^{\rm inj},$ (41)

where $t_{\rm esc}^{e}$ is the electron escape time scale, $P_{\rm ad}$ is the
adiabatic energy loss rate, $P_{\rm syn}^{e}$ is the synchrotron energy loss
rate, $P_{\rm IC}^{e}$ is the inverse-Compton energy loss rate, $\partial
n_{\varepsilon_{e}}^{\rm BH}/\partial t$ is the electron-positron generation
rate from the Bethe-Heitler pair production process, $\partial
n_{\varepsilon_{e}}^{\rm phmes}/\partial t$ is the electron-positron
generation rate from the photomeson production process,
$n_{\varepsilon_{e}}^{\beta_{\rm dec}}$ is the electron generation rate from
$\beta$ decay, $\dot{n}_{\varepsilon_{e}}^{\rm inj}$ is the primary electron
injection rate.

For electromagnetic cascades inside the source, we adopt the continuous
energy-loss approximation, as in Murase et al. (2015), where the second term
in the right hand of Eq. A.2 can be rewritten as
$\frac{\partial}{\partial\varepsilon_{e}}\left[(P_{\rm ad}+P_{\rm
syn}^{e}+P_{\rm
IC}^{e})n_{\varepsilon_{e}}^{e}\right]=\frac{\partial}{\partial\varepsilon_{e}}\left[P_{\rm
ad}+P_{\rm syn}^{e}+P_{\rm
IC}^{e}\right]n_{\varepsilon_{e}}^{e}+\frac{\partial}{\partial\varepsilon_{e}}\left[n_{\varepsilon_{e}}^{e}\right](P_{\rm
ad}+P_{\rm syn}^{e}+P_{\rm IC}^{e})$. We can express the coefficient
$\mathcal{A}$ as

$\mathcal{A}_{e}=\frac{1}{t_{\rm
esc}^{e}}+\frac{\partial}{\partial\varepsilon_{e}}\left[P_{\rm ad}+P_{\rm
syn}^{e}+P_{\rm IC}^{e}\right],$ (42)

where the discretization formula is

$\mathcal{A}_{e,i}=\frac{1}{t_{\rm
esc}^{e}}\bigg{|}_{\varepsilon=\varepsilon_{i}}+\frac{1}{\varepsilon_{i+1/2}-\varepsilon_{i-1/2}}\left(P_{\rm
ad}+P_{\rm syn}^{e}+P_{\rm IC}^{e}\right)|_{\varepsilon=\varepsilon_{i}}.$
(43)

The coefficient $\mathcal{B}$ can be written as

$\mathcal{B}_{e,j\to
e,i}=\delta_{j}^{i+1}\frac{1}{\varepsilon_{j+1/2}-\varepsilon_{j-1/2}}\left(P_{\rm
ad}+P_{\rm syn}^{e}+P_{\rm IC}^{e}\right)|_{\varepsilon=\varepsilon_{j}},$
(44)

where $\delta$ is the Kronecker delta function, the index $i$ and $j$
represent electron energy index (e.g., Murase & Beacom, 2012; Kalashev & Kido,
2015). The coefficient $\mathcal{C}$ can be written as

$\mathcal{C}_{e}=\frac{\partial}{\partial
t}\left(n_{\varepsilon_{\gamma}}^{\gamma\gamma}+n_{\varepsilon_{\gamma}}^{\rm
BH}+n_{\varepsilon_{\gamma}}^{\rm phmes}+n_{\varepsilon_{\gamma}}^{\beta_{\rm
dec}}\right).$ (45)

In general, the escape time depends on the diffusion coefficient or details of
magnetic fields. In the limit that the magnetic confinement is sufficiently
long, which is valid for electrons in the energy range of interest, one may
approximate the escape time scale by

$t_{\rm esc}^{e}\approx t_{\rm adv}=\frac{l_{b}}{V},$ (46)

where $t_{\rm adv}$ is the advection time scale and $V$ is the expansion speed
of the emission region. The adiabatic energy loss time scale can be estimated,

$t_{\rm
ad}^{-1}\equiv\frac{1}{\varepsilon_{e}}\frac{d\varepsilon_{e}}{dt}\approx\left(\frac{l_{b}}{V_{\rm
adv}}\right)^{-1},$ (47)

where the adiabatic energy loss rate is $P_{\rm ad}=\varepsilon_{e}t_{\rm
ad}^{-1}$, $l_{b}$ is the comoving size of the blob, and $V_{\rm adv}\sim V$
is the advection velocity at the dissipation region.

###  A.3 Neutrino

The kinetic equation for neutrinos is

$\displaystyle\frac{\partial n_{\varepsilon_{\nu}}^{\nu}}{\partial t}$
$\displaystyle=-\frac{n_{\varepsilon_{\nu}}^{\nu}}{t_{\rm
esc}}+\frac{\partial}{\partial t}\left(n_{\varepsilon_{\nu}}^{\rm
phmes}+n_{\varepsilon_{\nu}}^{\beta_{\rm dec}}\right),$ (48)

where $t_{\rm esc}^{\nu}=t_{\rm lc}$ is the neutrino escape time scale,
$\partial n_{\varepsilon_{\nu}}^{\rm phmes}/\partial t$ is the neutrino
generation rate from the photomeson production process, and $\partial
n_{\varepsilon_{\nu}}^{\beta_{\rm dec}}/\partial t$ is the (anti-)electron
neutrino generation rate from $\beta$ decay.

The coefficient $\mathcal{A}$ is

$\mathcal{A}_{\nu}=\frac{1}{t_{\rm esc}^{\nu}}.$ (49)

The coefficient $\mathcal{B}$ is

$\mathcal{B}_{\nu}=0.$ (50)

The coefficient $\mathcal{C}$ is

$\mathcal{C}_{\nu}=\frac{\partial}{\partial t}\left(n_{\varepsilon_{\nu}}^{\rm
phmes}+n_{\varepsilon_{\nu}}^{\beta_{\rm dec}}\right).$ (51)

###  A.4 Neutron

The kinetic equation for neutrons is

$\displaystyle\frac{\partial n_{\varepsilon_{n}}^{n}}{\partial
t}=-\frac{n_{\varepsilon_{n}}^{n}}{t_{\rm
esc}^{n}}-\frac{n_{\varepsilon_{n}}^{n}}{t_{\rm
phmes}}-\frac{n_{\varepsilon_{n}}^{n}}{t_{\beta_{\rm
dec}}}+\frac{\partial}{\partial t}\left(n_{\varepsilon_{n}}^{\rm
phmes}+n_{\varepsilon_{n}}^{\rm phdis}\right),$ (52)

where $t_{\rm esc}^{n}=t_{\rm lc}$ is the neutron escape time scale, $t_{\rm
phmes}$ is the photomeson production time scale, $t_{\beta_{\rm dec}}$ is the
neutron lifetime, $\partial n_{\varepsilon_{n}}^{\rm phmes}/\partial t$ is the
neutron generation rate from the photomeson production process, and $\partial
n_{\varepsilon_{n}}^{\rm phdis}/\partial t$ is the neutron generation rate
from nuclear photodisintegration.

The coefficient $\mathcal{A}$ for neutron is

$\mathcal{A}_{n}=\frac{1}{t_{\rm esc}^{n}}+\frac{1}{t_{\rm
phmes}}+\frac{1}{t_{\beta_{\rm dec}}}.$ (53)

The coefficient $\mathcal{B}$ is,

$\mathcal{B}_{n}=\frac{\partial n_{\varepsilon_{n}}^{\rm phmes}}{\partial t},$
(54)

where $\partial n_{\varepsilon_{n}}^{\rm phmes}/\partial t$ is the generation
rate of the neutron from the photomeson production process when the primary
particle is a neutron. The coefficient $\mathcal{C}$ is

$\mathcal{C}_{n}=\frac{\partial}{\partial t}\left(n_{\varepsilon_{n}}^{\rm
phmes}+n_{\varepsilon_{n}}^{\rm phdis}\right),$ (55)

where $\partial n_{\varepsilon_{n}}^{\rm phmes}/\partial t$ is the generation
rate of neutrons from the photomeson production process when the primary
particle is a proton or a nucleus.

###  A.5 Proton

The kinetic equation for protons is

$\displaystyle\frac{\partial n_{\varepsilon_{p}}^{p}}{\partial t}$
$\displaystyle=-\frac{n_{\varepsilon_{p}}^{p}}{t_{\rm
esc}}-\frac{n_{\varepsilon_{p}}^{p}}{t_{\rm
phmes}}+\frac{\partial}{\partial\varepsilon_{p}}\left[(P_{\rm ad}+P_{\rm
BH}^{p}+P_{\rm syn}^{p})n_{\varepsilon_{p}}^{p}\right]$
$\displaystyle+\frac{\partial}{\partial t}\left(n_{\varepsilon_{p}}^{\rm
phmes}+n_{\varepsilon_{p}}^{\beta_{\rm dec}}+n_{\varepsilon_{p}}^{\rm
phdis}\right)+\dot{n}_{\varepsilon_{p}}^{\rm inj},$ (56)

where $t_{\rm esc}^{p}$ is the proton escape time scale, $t_{\rm phmes}$ is
the photomeson production interaction time scale, $P_{\rm BH}$ is the energy
loss rate due to the Bethe-Heitler pair production process, $P_{\rm syn}^{p}$
is the proton synchrotron energy loss rate, $\partial n_{\varepsilon_{p}}^{\rm
phmes}/\partial t$ is the proton generation rate from the photomeson
production process, $n_{\varepsilon_{p}}^{\beta_{\rm dec}}$ is the proton
generation rate from $\beta$ decay, $\partial n_{\varepsilon_{p}}^{\rm
phdis}/\partial t$ is the proton generation rate from nuclear
photodisintegration, $\dot{n}_{\varepsilon_{p}}^{\rm inj}$ is the primary
proton injection rate.

For the Bethe-Heitler pair production process, we adopt the continuous energy
loss approximation. The efficient $\mathcal{A}$ is

$\mathcal{A}_{p}=\frac{1}{t_{\rm esc}^{p}}+\frac{1}{t_{\rm
phmes}}+\frac{\partial}{\partial\varepsilon_{p}}(P_{\rm ad}+P_{\rm
BH}^{p}+P_{\rm syn}^{p}),$ (57)

where the discretization form is

$\mathcal{A}_{p,i}=\frac{1}{t_{\rm esc}^{p}}\bigg{|}_{i}+\frac{1}{t_{\rm
phmes}}\bigg{|}_{i}+\frac{1}{\varepsilon_{i+1/2}-\varepsilon_{i-1/2}}(P_{\rm
ad}+P_{\rm BH}^{p}+P_{\rm syn}^{p})|_{\varepsilon=\varepsilon_{i}}.$ (58)

The coefficient $\mathcal{B}$ is

$\mathcal{B}_{p,j\to p,i}=\frac{\partial n_{\varepsilon_{p}}^{\rm
phmes}}{\partial
t}+\delta_{j}^{i+1}\frac{1}{\varepsilon_{j+1/2}-\varepsilon_{j-1/2}}(P_{\rm
ad}+P_{\rm BH}^{p}+P_{\rm syn}^{p})|_{\varepsilon=\varepsilon_{j}},$ (59)

where $\partial n_{\varepsilon_{p}}^{\rm phmes}/\partial t$ is the generation
rate of protons from the photomeson production process when the primary
particle is a proton. The coefficient $\mathcal{C}$ is written as

$\mathcal{C}_{p}=\frac{\partial}{\partial
t}\left(n_{\varepsilon_{p}}^{\beta_{\rm dec}}+n_{\varepsilon_{p}}^{\rm
phmes}+n_{\varepsilon_{p}}^{\rm phdis}\right).$ (60)

For the escape processes, we consider both diffusion and advection. The
confinement time scale is given by

$t_{\rm conf}={\rm max}[t_{\rm diff},t_{\rm lc}],$ (61)

where

$t_{\rm diff}\approx\frac{l_{b}^{2}}{6D},$ (62)

is the diffusion time scale for a spherical blob geometry and $D$ is the
diffusion coefficient, see Eq. 22. The escape time scale can be estimated as

$t_{\rm esc}^{p}={\rm min}[t_{\rm conf},t_{\rm adv}],$ (63)

where $t_{\rm adv}=l_{b}/V$ is the advection escape time scale.

###  A.6 Nuclei

The kinetic equation for nuclei is

$\displaystyle\frac{\partial n_{\varepsilon_{A}}^{A}}{\partial t}$
$\displaystyle=-\frac{n_{\varepsilon_{A}}^{A}}{t_{\rm
esc}^{A}}-\frac{n_{\varepsilon_{A}}^{A}}{t_{\rm
phmes}}-\frac{n_{\varepsilon_{A}}^{A}}{t_{\rm
phdis}}+\frac{\partial}{\partial\varepsilon_{A}}\left[(P_{\rm ad}+P_{\rm
BH}^{A}+P_{\rm syn}^{A})n_{\varepsilon_{A}}^{A}\right]$
$\displaystyle+\frac{\partial}{\partial t}\left(n_{\varepsilon_{p}}^{\rm
phmes}+n_{\varepsilon_{A}}^{\rm phdis}\right)+\dot{n}_{\varepsilon_{A}}^{\rm
inj},$ (64)

where $t_{\rm esc}^{A}$ is the nuclear escape time scale, $t_{\rm phmes}$ is
the photomeson production interaction time scale, $t_{\rm phdis}$ is the
photodisintegration production interaction time scale, $P_{\rm BH}$ is the
energy loss rate due to the Bethe-Heitler pair production process, $P_{\rm
syn}^{A}$ is the nuclei synchrotron energy loss rate, $\partial
n_{\varepsilon_{A}}^{\rm phmes}/\partial t$ is the nuclear generation rate
from the photomeson production process, $\partial n_{\varepsilon_{A}}^{\rm
phdis}/\partial t$ is the nuclear generation rate from nuclear
photodisintegration process, $\dot{n}_{\varepsilon_{A}}^{\rm inj}$ is the
injection rate of primary nuclei.

For the Bethe-Heitler pair production process, we adopt the continuous energy
loss approximation similar to the proton case. The coefficient $\mathcal{A}$
for high-energy nuclei is

$\mathcal{A}_{A}=\frac{1}{t_{\rm esc}^{A}}+\frac{1}{t_{\rm
phmes}}+\frac{1}{t_{\rm phdis}}+\frac{\partial(P_{\rm ad}+P_{\rm
BH}^{A}+P_{\rm syn}^{A})}{\partial\varepsilon_{p}}.$ (65)

where the discretization form is

$\mathcal{A}_{A,i}=\frac{1}{t_{\rm esc}^{A}}\bigg{|}_{i}+\frac{1}{t_{\rm
phmes}}\bigg{|}_{i}+\frac{1}{t_{\rm
phdis}}\bigg{|}_{i}+\frac{1}{\varepsilon_{i+1/2}-\varepsilon_{i-1/2}}(P_{\rm
ad}+P_{\rm BH}^{A}+P_{\rm syn}^{A})|_{\varepsilon=\varepsilon_{i}}.$ (66)

The coefficient $\mathcal{B}$ is

$\mathcal{B}_{A,j\to
A,i}=\delta_{j}^{i+1}\frac{1}{\varepsilon_{j+1/2}-\varepsilon_{j-1/2}}(P_{\rm
ad}+P_{\rm BH}^{A}+P_{\rm syn}^{A})|_{\varepsilon=\varepsilon_{j}}.$ (67)

The coefficient $\mathcal{C}$ is written as

$\mathcal{C}_{A}=\frac{\partial}{\partial t}\left(n_{\varepsilon_{A}}^{\rm
phmes}+n_{\varepsilon_{A}}^{\rm phdis}\right).$ (68)

Similar to protons, the escape time scale of nuclei depends on both diffusion
and advection.