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

TL;DR

This paper develops a numerical model to study how ultrahigh-energy cosmic rays and their secondary particles produce gamma rays and neutrinos in radio galaxies, focusing on Centaurus A, and predicts detectable gamma-ray signatures from UHECR nuclei.

Contribution

It introduces a self-consistent numerical code for simulating nuclear and electromagnetic cascades from UHECRs in radio galaxies, specifically applied to Centaurus A.

Findings

VHE gamma-ray signatures from UHECR nuclei are detectable with current telescopes.

Intermediate-mass nuclei like oxygen can produce observable gamma-ray signals.

Interactions with photon fields lead to photodisintegration and pair production processes.

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

Tables

Table 1. Table 1: Physical parameters in the comoving frame used in the leptohadronic model for the kiloparsec-scale jet.

Parameter	Value
Viewing angle (${\theta }_{\mathrm{ob}}\left[\mathrm{deg}\right]$)	40
Lorentz factor (${\Gamma }_{\mathrm{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 (${\epsilon }_{e,\min }$ [eV])	$2.5\times {10}^7$
Primary electron maximum energy (${\epsilon }_{e,\max }$ [eV])	$2\times {10}^{12}$
Electron injection energy density ($u_{\mathrm{inj},e}[\mathrm{eV}{s}^{-1}$]	$3$
Oxygen minimum energy (${\epsilon }_{O,\min }$ [eV])	$1\times {10}^{17}$
Oxygen maximum energy (${\epsilon }_{O,\max }$ [eV])	$4\times {10}^{18}$
Oxygen spectral index ($s_{\mathrm{acc}}$)	2.6
Oxygen injection energy density ($u_{\mathrm{inj},O}[\mathrm{eV}{\mathrm{cm}}^{-3}$]	$300$
Iron minimum energy (${\epsilon }_{\mathrm{Fe},\min }$ [eV])	$1\times {10}^{17}$
Iron maximum energy (${\epsilon }_{\mathrm{Fe},\max }$ [eV])	$13\times {10}^{18}$
Iron spectral index ($s_{\mathrm{acc}}$)	2.6
Iron injection energy density ($u_{\mathrm{inj},\mathrm{Fe}}[\mathrm{eV}{\mathrm{cm}}^{-3}$]	$100$

Equations

$$q({\epsilon }_A,t)\equiv \frac{d{N}^A}{dtd{\epsilon }_A}={\dot{N}}_0^A{\left(\frac{{\epsilon }_A}{Z{\epsilon }_A^0}\right)}^{-s_{\mathrm{acc}}}\exp \left(-\frac{{\epsilon }_A}{Z{\epsilon }_p^{\max }}\right),$$

$$L_{\mathrm{CR}}=\frac{3}{2}{\Gamma }^2\sum _A\int d{\epsilon }_A\frac{d{N}^A}{dtd{\epsilon }_A},$$

$$\frac{\partial n_{{\epsilon }_a}^a}{\partial t}$$

$$+\int d{\epsilon }_a^{\ast }{n}_{{\epsilon }_a^{\ast }}^a{\mathcal{B}}_{a\to a}({\epsilon }_a,{\epsilon }_a^{\ast })$$

$$+\sum _b\int d{\epsilon }_b{n}_{{\epsilon }_b}^b{\mathcal{C}}_{b\to a}({\epsilon }_a,{\epsilon }_b)$$

$$+{\dot{n}}_a^{\mathrm{inj}}({\epsilon }_a),$$

$$t_{A\gamma }^{-1}$$

$$=\frac{c}{2}{\int }_{-1}^{1}d\mu {\int }_0^{\infty }d\epsilon (1-{\beta }_A\mu )\frac{dn}{d\epsilon }{\sigma }_{A\gamma }(\bar{\epsilon }),$$

$${\tau }_{A\gamma }\approx t_{\mathrm{esc}}/t_{A\gamma },$$

$$t_{\mathrm{esc}}\approx t_{\mathrm{adv}}\approx \frac{l_b}{V},$$

$$\frac{dn}{d\epsilon }=n_{{\epsilon }_b}\{\begin{array}{l}(\epsilon /{\epsilon }_b{)}^{-{\alpha }_l},\epsilon <{\epsilon }_b\\ (\epsilon /{\epsilon }_b{)}^{-{\alpha }_h},\epsilon >{\epsilon }_b\end{array},$$

$$A+\gamma \to A_d^{\star }+\mathrm{nucleons},$$

$${\tau }_{A\gamma }$$

$$\approx \frac{2}{\alpha +1}(c/V)l_b{\hat{\sigma }}_{\mathrm{dis}}{\epsilon }_b{n}_{{\epsilon }_b}{\left(\frac{E_A}{E_{A,b}^{\mathrm{dis}}}\right)}^{\alpha -1}$$

$$\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\mathrm{kpc}}\right){\left(\frac{E_A}{E_{A,b}^{\mathrm{dis}}}\right)}^{\alpha -1},$$

$$f_{\mathrm{dis}}$$

$$\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\mathrm{kpc}}\right){\left(\frac{E_A}{E_{A,b}^{\mathrm{dis}}}\right)}^{\alpha -1},$$

$$A_d^{\star }\to A_d+\mathrm{photons},$$

$$f_{\mathrm{deex}}$$

$$\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\mathrm{kpc}}\right){\left(\frac{E_A}{E_{A,b}^{\mathrm{dis}}}\right)}^{\alpha -1},$$

$$f_{\mathrm{mes}}$$

$$\sim 0.3\times 10^{-3}\frac{2}{\alpha +1}{\left(\frac{V}{0.5c}\right)}^{-1}\left(\frac{l_b}{1\mathrm{kpc}}\right){\left(\frac{E_A}{E_{A,b}^{\mathrm{mes}}}\right)}^{\alpha -1},$$

$$f_{\mathrm{BH}}$$

$$\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\mathrm{kpc}}\right){\left(\frac{E_A}{E_{A,b}^{\mathrm{BH}}}\right)}^{\alpha -1},$$

$$t_{\mathrm{syn}}^{-1}\approx \frac{4{\sigma }_T{Z}^4{m}_e^2}{3m^4{c}^3}\frac{E_A}{{\delta }_D}{U}_B,$$

$$t_{\mathrm{IC}}^{-1}\approx \frac{4{\sigma }_{T}c}{3m_e^2{c}^4}\frac{E_e}{{\delta }_D}{U}_{\mathrm{ph}}{F}_{\mathrm{KN}},$$

$${\tau }_{\gamma \gamma }$$

$$\sim 0.2\left(\frac{l_b}{1\mathrm{kpc}}\right){\left(\frac{E_{\gamma }}{E_{\gamma ,b}}\right)}^{\alpha -1},$$

$$E_{\max }$$

$$\simeq 110(Z/26){\Gamma }_{\mathrm{kpc}}{\left(\frac{\eta }{7}\right)}^{-1}\left(\frac{B}{10^{-4}\mathrm{G}}\right)\left(\frac{\beta }{0.3}\right)\left(\frac{l_b}{1\mathrm{kpc}}\right)\mathrm{EeV},$$

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

$${\epsilon }^2\frac{dn}{d\epsilon }=\frac{d_L^2}{R^{2}c}\frac{{\delta }_{D,\mathrm{core},\mathrm{kpc}}^4}{{\delta }_{D,\mathrm{core},\mathrm{ob}}^4}{\delta }_{D,\mathrm{deboost}}^{4}E{F}_E,$$

$$E{F}_E=A{x}^{-{\alpha }_1}{\left(\frac{1+x^{1/{\alpha }_3}}{2}\right)}^{({\alpha }_1-{\alpha }_2){\alpha }_3},$$

$$D(E)\approx \frac{c{l}_{\mathrm{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],$$

$$E_c\approx ZeB{l}_{\mathrm{coh}}\simeq 9.6Z\frac{B}{10^{-4}\mathrm{G}}\frac{l_{\mathrm{coh}}}{0.1\mathrm{kpc}}\mathrm{EeV},$$

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

Full text

††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,

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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.

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

[TABLE]

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),

[TABLE]

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

[TABLE]

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,

[TABLE]

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)

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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),

[TABLE]

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)

[TABLE]

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

[TABLE]

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,

[TABLE]

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

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

[TABLE]

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

[TABLE]

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)

[TABLE]

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

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.

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

[TABLE]

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),

[TABLE]

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

[TABLE]

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

[TABLE]

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)

[TABLE]

where

[TABLE]

$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)

[TABLE]

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,

[TABLE]

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

[TABLE]

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.

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.

Appendix A Numerical method and related physical processes

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

[TABLE]

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),

[TABLE]

and we have

[TABLE]

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).

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)

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where the discretization formula is

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

The coefficient $\mathcal{B}$ is

[TABLE]

The coefficient $\mathcal{C}$ is

[TABLE]

A.4 Neutron

The kinetic equation for neutrons is

[TABLE]

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

[TABLE]

The coefficient $\mathcal{B}$ is,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where the discretization form is

[TABLE]

The coefficient $\mathcal{B}$ is

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

A.6 Nuclei

The kinetic equation for nuclei is

[TABLE]

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

[TABLE]

where the discretization form is

[TABLE]

The coefficient $\mathcal{B}$ is

[TABLE]

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

[TABLE]

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