High-Energy Neutrino Emission from Espresso-Reaccelerated Ions in Jets of Active Galactic Nuclei

We present a bottom-up calculation of the flux of ultra-high energy cosmic rays (UHECRs) and high-energy neutrinos produced by powerful jets of active galactic nuclei (AGNs). By propagating test particles in 3D relativistic magnetohydrodynamic jet simulations, including a Monte Carlo treatment of sub-grid pitch-angle scattering and attenuation losses due to realistic photon fields, we study the spectrum and composition of the accelerated UHECRs and estimate the amount of neutrinos produced in such sources. We find that UHECRs may not be significantly affected by photodisintegration in AGN jets where the espresso mechanism efficiently accelerates particles, consistent with Auger's results that favor a heavy composition at the highest energies. Moreover, we present estimates and upper bounds for the flux of high-energy neutrinos expected from AGN jets. In particular, we find that: i) source neutrinos may account for a sizable fraction, or even dominate, the expected flux of cosmogenic neutrinos; ii) neutrinos from the beta-decay of secondary neutrons produced in nucleus photodisintegration end up in the TeV-PeV band observed by IceCube, but can hardly account for the observed flux; iii) UHECRs accelerated via the espresso mechanism lead to nearly isotropic neutrino emission, which suggests that nearby radio galaxies may be more promising as potential sources. We discuss our results in the light of multimessenger astronomy and current/future neutrino experiments.


INTRODUCTION
A comprehensive theory that accounts for particle injection, particle acceleration, and spectra of ultrahigh energy cosmic rays (UHECRs) with energies above 10 18 eV is still lacking. The origin of these particles, however, has been the subject of many theoretical studies that rely mostly on estimates of the maximum energy that particles can achieve in specific environments (Cavallo 1978;Hillas 1984;Alves Batista et al. 2019). In that regard, one of the most promising astrophysical sources of UHECRs are AGN jets, which satisfy the Hillas criterion up to 10 20 eV, especially if the highestenergy CRs are heavy ions. Additionally, their luminosities can explain the energy injection rate necessary to sustain the UHECR flux (e.g., Katz et al. 2009;Murase & Fukugita 2019;Jiang et al. 2021). Other sources have also been suggested such as newly-born millisecond pul-rmbarek@uchicago.edu sars (e.g., Fang et al. 2012Fang et al. , 2014, γ-ray bursts (e.g., Vietri 1995;Waxman 1995), engine-driven transrelativistic supernovae including low-luminosity γ-ray bursts and hypernovae (e.g., Murase et al. 2006;Wang et al. 2008;Zhang & Murase 2019), and tidal disruption events (e.g., Farrar & Piran 2014;Zhang et al. 2017). However, the exact acceleration mechanism boosting UHECRs to their energies in these environments is not pinned down.
This paper is the third in a series of projects aimed to analyze the promising espresso model (Caprioli 2015), checking if it may satisfy all the requirements of a particle acceleration theory and outlining its observational predictions. In essence, the espresso framework suggests that UHECRs can be produced in relativistic AGN jets via the reacceleration of galactic CR seeds. Such seeds, accelerated in supernova remnants up to a few PV in rigidity, penetrate in the highly relativistic regions of the jets and tap in their radial electric field to receive one, or even multiple, ∼ Γ 2 boosts in energy. If the jet is sufficiently powerful with Γ ∼ 20 − 30, a single shot would allow them to reach UHECR energies. In Mbarek & Caprioli (2019); Mbarek & Caprioli (2021), the espresso mechanism has been tested by propagating particles in high-resolution magnetohydrodynamical (MHD) simulations of AGN jets (Mignone et al. 2007;Rossi et al. 2008). The reacceleration of galactic CR seeds is also promising even for subrelativistic AGN jets that are seen at kpc scales, and the model can fit the Auger spectrum and composition data simultaneously (Kimura et al. 2018).
In Mbarek & Caprioli (2019), hereafter MC19, we found that the spectra, chemical composition, and anisotropy of the reaccelerated particles are qualitatively consistent with UHECR phenomenology. Then, in Mbarek & Caprioli (2021), hereafter MC21, we included sub-grid scattering (SGS) to model small-scale magnetic turbulence that cannot be resolved by MHD simulations, constraining for the first time one potentially crucial but hard-to-model ingredient. We established the relative importance of espresso and stochastic shear acceleration in relativistic jets, finding that strong SGS, on one hand, can promote the injection and acceleration of lower-energy UHECRs, but on the other hand, is irrelevant for the acceleration of the highestenergy CRs, which are invariably espresso-accelerated in the inner regions of relativistic jets. However, we should also keep in mind that shear acceleration may play a more important role in accelerating particles in subrelativistic jets that are usually seen in kiloparsec scales especially for Fanaroff-Riley (FR) I galaxies. Overall, we expect the neutrino production to be relatively important in the relativisitic spine of the jet, in particular powerful jets seen in FR II galaxies and flat spectrum radio quasars (FSRQs), where particles are promoted to the highest energies.
In this paper, we evaluate the effects of photodisintegration and high-energy neutrino production in AGN jets within our self-consistent particle acceleration framework. In particular, we investigate how the intense radiation fields of the blazar zone, the broadline region, and the dusty torus may affect the chemical composition of the accelerated particles. Moreover, modeling UHECR attenuation in a realistic jet environment allows us to calculate the spectrum of high-energy neutrinos produced in these sources.
Within our bottom-up approach we aim to address, with as few assumptions as possible, some key open questions such as: • What are the effects of losses on espressoaccelerated particles in sub-kiloparsec-scale AGN jets?
• What is the expected spectrum of UHE (Ultra-High-Energy) neutrinos produced by a typical AGN jet?
• If UHECRs are accelerated in AGN jets via the espresso mechanism, can they be responsible for the observed IceCube flux, too?
This work would be particularly important to unravel questions associated with Ultra-High-Energy (UHE) neutrinos. UHE neutrinos are created through interactions of UHECRs and are pivotal tools to advance our knowledge of extreme astrophysical environments. Many current and proposed experiments, such as the balloon-borne interferometer ANITA (Gorham et al. 2018a,b), its successor the Payload for Ultrahigh Energy Observations (PUEO) (Abarr et al. 2021), the Askaryan Radio Array (ARA) (Allison et al. 2012(Allison et al. , 2016, the In-Ice Radio Array ARIANNA (Barwick et al. 2017), the exciting next-generation IceCube detector IceCube-Gen2 (IceCube-Gen2 Collaboration et al. 2014), the Giant Radio Array for Neutrino Detection (GRAND) (Álvarez-Muñiz et al. 2020), the Radio Neutrino Observatory in Greenland (RNO-G) (Aguilar et al. 2021), the Beamforming Elevated Array for Cosmic Neutrinos (BEACON) (Wissel et al. 2020), and the proposed PO-EMMA mission (Olinto et al. 2017), all aim to detect EeV neutrinos for the first time. There are many studies that focus on cosmogenic neutrinos, i.e., the neutrinos created by UHECRs through interactions with the extragalactic photon background during intergalactic propagation (e.g. Beresinsky & Zatsepin 1969;Yoshida & Teshima 1993;Takami et al. 2009;Heinze et al. 2016;Romero-Wolf & Ave 2018;Heinze et al. 2019;Das et al. 2019;Wittkowski & Kampert 2019). However, UHE source neutrinos, i.e., neutrinos produced in or around UHECR accelerators, especially in the presence of extreme photon fields, could also be crucial to unravel the sites of production of the highest-energy particles in the Universe (see Batista et al. 2019, and references therein). This study aims to shed more light on these UHE source neutrinos.
In this paper, we propagate particles in an MHD simulation of an AGN jet to keep assumptions to a minimum and find that particles are espresso-accelerated (MC19; MC21). The neutrino spectrum that ensues is presented in §3.
The paper is organized as follows. In §2 we describe our particle acceleration framework, detailing the different interaction routes that lead to losses and neutrino production. In §3 we investigate the effects of losses on the UHECR chemical composition and put constraints on the expected upper bounds of the neutrino spectrum resulting from UHECR interactions. We also discuss the acceleration mechanism responsible for boosting UHE-CRs that may contribute to the IceCube flux, and finally present our conclusions in §4.

MHD simulation of a relativistic jet
To facilitate the comparison with published results, we model the underlying AGN relativistic jet via the same benchmark simulation used in MC19 and MC21; we refer to those papers for all the details and summarize here the essential information. We consider a 3D relativistic MHD simulation of a powerful AGN jet performed with PLUTO (Mignone et al. 2010), which includes adaptive mesh refinement. The jet, with a magnetization radius R jet , is initialized with Lorentz factor Γ 0 = 7 along the z-direction in a box that measures 48R jet in the x-and y-directions and 100R jet in the z-direction in a grid that has 512 × 512 × 1024 cells with four refinement levels. The jet/ambient density contrast is set to ψ = 10 −3 , the jet sonic and Alfvénic Mach numbers to M s =3, and M A = 1.67 respectively. Once the jet has developed, the effective Lorentz factor in the jet spine is Γ eff ∼ 3.2; this value is important to establish how many Γ 2 shots a particle undergoes during acceleration.

Particle propagation
We propagate ∼ 10 5 test particles in a snapshot of the benchmark jet with a broad range of initial gyroradii R and positions. We include the effects of unresolved turbulence by setting the SGS mean free path to be as small as the particle's gyroradius (Bohm diffusion) to maximize the number of particles within the jet spine and boost the efficiency of particle acceleration (see MC21 for more details on the effects of SGS). This prescription maximizes the particles residence time in the jet (and hence the effects of photon fields), but does not affect particle acceleration at the highest energies, which are invariably accelerated via the SGS-independent espresso mechanism (see MC21).
In addition to protons, we consider four different seed ion species, labelled se =[He, C/N/O,Mg/Al/Si, Fe] with effective atomic number Z se = [2, 7, 13, 26] and mass A se = [4,14,27,56], respectively. The energy flux of these seed galactic CRs below the knee is parameterized as follows: We set the normalizations according to the abundance ratios at 10 12 eV observed in Galactic CRs, such that K se /K H ∼ [0.46, 0.30, 0.07,0.14]. The motivation for using Galactic CR fluxes as fiducial cases hinges on the fact that the knee feature (the maximum seed energy) should be quite independent of the galaxy mass (Caprioli 2015). We acknowledge that the actual seed fluxes in AGN hosts may be different due to the different injection and confinement properties of other galaxies but, since the chemical enrichment provided by diffusive shock acceleration of seeds in supernova remnants (Caprioli et al. 2010(Caprioli et al. , 2017 is rather universal and since the final UHECR fluxes that we consider are normalized to the observed one, the assumptions above are quite generic. The spectra of different species are calculated as outlined in Appendix A. Particles are propagated in the MHD simulation even if their gyroradii are smaller than the grid size, which is reasonable as long as a sufficiently small rigiditydependent timestep is used to resolve their gyration; every particle's gyroradius is resolved with at least 10 timesteps.

Photon field prescriptions
On top of the jet structure provided by the MHD simulation, we prescribe external photon fields based on the methods presented in Murase et al. (2014), hereafter MID14. There is ostensible uncertainty in modeling these external components due to the vast AGN diversity, but the systematic approach of MID14 allows us to assess the individual effect of such fields based on the apparent bolometric luminosity of the jet.
We consider five different photon backgrounds of different origin: (i) Non-thermal emission: -It originates from synchrotron and inverse-Compton radiation of relativistic leptons and/or hadronic emission, emerging from the blazar zone, a sub-pc region close to the base of the jet (see MID14 for more details), where the emission is dominated by X/γ-rays and is relativistically beamed with an angle ∼ 1/Γ in the black hole frame. We refer to this broadband emission region as the non-thermal cone of influence (see Figure 1).
(ii) Radiation from the broad atomic line region (BLR): -This is the reprocessed emission from cold gas clumps photoionized by the UV and X-ray produced by the accretion disk. These sub-pc spherical clumps are located closer to the base of the jet (< 1pc away) and have a luminosity ∼ 10% of that of the accretion disk (see MID14).
(iii) IR emission from the dusty torus: -This is IR from reprocessed accretion dusty disk radiation with a torus size that can reach ∼ 1pc. Following MID14, we model it as a spherical grey body with temperature ∼500K.
(iv) Stellar light: -The photons from the host-galaxy stars have been shown to have large energy densities compared to other photon fields at a few hundred pc, which makes them important targets for accelerated particles in powerful AGN jets. In the remainder of the paper, we will consider the starlight emission profile for Centaurus A as our fiducial case (Tanada et al. 2019).
(v) The cosmic microwave background (CMB) radiation: -Besides affecting every particle regardless of its location, accounting for the CMB contribution serves as a benchmark to compare the effect of the prescribed photon fields.
In the remainder of the paper, all photon fields are assumed to be isotropic except for the non-thermal component, which is beamed within a cone of aperture 1/Γ eff , where Γ eff is the effective Lorentz factor of the jet. The BLR, IR and non-thermal contributions have intensities that are inversely proportional to the square of the distance from their emitting regions.

Particle interactions
When propagating our test particles, at every time step we: (i) calculate the probability of interaction with the thermal plasma (assuming it is electron-proton) and photon fields; (ii) keep track of each particle's atomic mass A and charge Z; (iii) monitor secondary particle production including neutrinos and secondary protons. assumed to escape without experiencing further interactions; their place of production and escaping direction are recorded, though. In order to study the effects of UHECR photodisintegration and the resulting neutrino flux, we consider the most relevant attenuation mechanisms, i.e., inelastic proton-proton (pp), photomeson production (pγ), and neutron decay following photodisintegration of heavy nuclei. More specifically, photomeson production inter-
actions of photons and nucleons also result in the production of pions that subsequently decay to create photons and neutrinos. On the other hand, photodisintegration interactions are nuclear processes that cause the photon-absorbing nucleus to change to another chemical specie and release either a proton or a neutron. Table 1 summarizes the interaction routes that lead to neutrinos and photodisintegration of heavy elements. For simplicity, we ignore the Bethe-Heitler process in this work, which is sufficient for our current purpose. This process can be important for setting the maximum energy only for very powerful blazars (MID14). A more detailed explanation of the interaction probabilities at every time step is included below.

Proton-proton (pp) interactions
Accelerated particles can experience pp scattering to create charged pions and hence ν e and ν µ neutrinos. At every time step, depending on the particle's energy E p and position x, there is an interaction probability P (E p , x) ∼ n(x)σ pp (E p )R jet ∆t, where σ pp is the pp interaction cross section (Tanabashi et al. 2018), n is the position-dependent density, and ∆t is the time step. The neutrino spectrum that results from every pp interaction is calculated based on the parametrization by Kelner et al. (2006).
Powerful AGN jets are often inside clusters of galaxies (e.g., Begelman et al. 1984;Best et al. 2007;Fang & Murase 2018), so that the ambient medium density should reflect that of the intra-cluster medium, which is of order of n ICM ∼ 10 −3 cm −3 (Walg et al. 2013). This reference value is used in Figure 1, which also shows that the jet itself is expected to have an even lower density. Generally, pp interactions inside the AGN jets are not expected to contribute much to the overall neutrino spectrum, although they can be relevant for cosmic rays that have escaped from AGN jets and are confined in clusters (Fang & Murase 2018).

Photomeson production (pγ) Interactions
At every time step ∆t, a photopion production probability f pγ is calculated such that f pγ = t −1 pγ ∆t, where t pγ is the photomeson cooling time. A detailed account of t pγ calculations for isotropic photon fields and interactions at a known angle for both protons and nuclei with atomic mass A is provided in Appendix B. For pγ interactions, we assume a multiplicity of 1 for neutrinos (with energies E ν ∼ E p /20) produced in one interaction of a relativistic proton, which is valid when the resonance and direct production are relevant. This assumption is valid when target photon spectra are sufficiently soft (Murase et al. 2008).
In general, the cooling time depends on the particle position and on the AGN luminosity; therefore, we cannot use dimensionless quantities but need to introduce physical scales for the magnetic field strength and for the jet size and luminosity. For instance, the left panel of Figure 2 shows the cooling time for protons of different energies located 750 pc away from the base of the jet along the spine (z-direction in Figure 1), for a jet with radius R jet = 15pc and an apparent bolometric luminosity L bol = 10 48 erg s −1 . In this characteristic example, we show the contribution from isotropic photon fields (BLR, IR, stellar light, and CMB, as in the legend) and the angle-dependent interactions with the beamed nonthermal continuum (cyan and purple lines). The left panel of Figure 2 shows that the non-thermal contribution provides the shortest photomeson production cooling time (even for tail-on interactions, i.e., θ = 0). These cooling curves depend on the particle position (distance from the base of the jet and position with respect to the non-thermal continuum cone of influence) as will be further discussed below.
The right panel of Figure 2 shows the dependence of the cooling time on the photon energy when the proton energy is fixed. Only the non-thermal contribution is shown here for simplicity as it is the dominant photon field in our fiducial cases. We can see that the shortest cooling time-which depends on the proton energyoccurs at the threshold energy¯ th for photomeson production interactions (See Appendix B for more details). Importantly, the contribution of X/γ-ray photons with ≥ 100eV to cooling is only significant for lower energy protons.

Photodisintegration interactions
On the same photon fields, nuclei with atomic mass A can also undergo photodisintegration with probability t −1 Aγ ∆t per time step. Appendix C details our calculations, in particular the use of the giant dipole resonance (GDR) cross section (e.g., Murase et al. 2008;Wang et al. 2008) as a fiducial case (see Figure 3). Although this is the simplest application, it is sufficient to demonstrate the effects of photodisintegration in our numerical work. The GDR approximation is valid for soft photon spectra (Murase et al. 2008), and non-thermal photon fields are dominant in our examples. More detailed studies including quasi-deutron emission and fragmentation are left as future work. We also account for photomeson production interactions of heavy nuclei (as described in Equation B6 in Appendix B), and find that their cooling time is comparable to that of protons. Also for nuclei, the most important photon background is typically the non-thermal component, as shown by a comparison of Figure 3 with Figure 2 and Equation B6. Higher-energy nuclei are more likely to be photodisintegrated, as expected. Note that the thresholds for photodisintegration interactions for He look quite close to those of photomeson interactions because they are plotted as a function of E A , hence a scaling by the atomic mass A is necessary.
While propagating particles, we keep track of the atomic mass and charge of nuclei, considering that after a photodisintegration event, the nucleus loses either one neutron or one proton. A neutron produced as a result of photodisintegration decays within a distance 9.15(E n /10 9 GeV) kpc (Anchordoqui et al. 2007) to produce one neutrino with energy 0.48 MeV in the neutron rest frame (Murase & Beacom 2010) such that E n /E ν ≈ 2 × 10 3 in the lab frame. Since secondary particles cannot achieve energies beyond ∼ 10 10 GeV, we assume that all the secondary neutrons β-decay and produce neutrinos.

Dependence on the distance from the central black hole
Since different photon backgrounds have different spatial extents, the cooling time for both photodisintegration and pγ interactions depend on the magnitude of the distance from the base of the jet, D, such that close to the base of the jet t −1 ∝ D −2 where photon fields other than the CMB are most relevant. For a jet with isotropic-equivalent bolometric luminosity L bol = 10 45 erg s −1 , non-thermal emission is dominant until D ∼ 3kpc, beyond which the CMB becomes important. Considering that the photon fields we assume are generated close to the base of the jet, CR interactions beyond a few kpc would not increase with large jet extents, usually associated with the more luminous FR-II jets.  Figure 2 but for photodisintegration interactions of nuclei with energy EA, including photomeson interactions and interactions based on the GDR photodisintegration total cross section (see Equation C10 and C12). Note that the threshold for He photodisintegration is quite close to that of photomeson production because the plot is a function of EA and a factor of A must be taken into consideration.

From scale-free MHD simulations to realistic environments
In our simulations, CR gyroradii are normalized to the jet radius R jet and magnetic field B 0 ; therefore, setting a physical value to R jet and B 0 is equivalent to associating physical energies to the seed particles of charge Ze. In order to calculate the actual neutrino fluxes from realistic AGNs, we need to fix a reference magnetic field and two physical quantities of the jet: its bolometric luminosity, which controls the photon fields, and its radius.
Our simulations have extensively shown that particles are routinely espresso-accelerated up to the jet Hillas limit (MC19, MC21). In our 3D relativistic MHD simulations, we launch the jet with Γ = 7; yet, the effective Lorentz factor of evolved jets turns out to be Γ eff ∼ a few, too small to promote CR seeds with rigidities of a few PeV to actual UHECRs with a single Γ 2 eff boost. Therefore, in order to achieve realistic UHECR energies, we fix the normalization of our jet radius and magnetization such that CR seeds have rigidities as large as 3 × 10 17 V, two orders of magnitude above the CR knee. This choice allows us to include the attenuation losses discussed for realistic photon fields and to calculate the fluxes of UHECRs and HE neutrinos expected from different types of AGNs.
When contemplating assigning magnetic field, radius, and L bol prescriptions to our jet, one needs to consider different types of radio-loud AGNs, which at minimum can be split into FR-I and FR-II sources. FR-I jets are typically decelerated to nonrelativistic bulk flows within 1 kpc (e.g., Wardle & Aaron 1997;Arshakian & Longair 2004;Mullin & Hardcastle 2009), while FR-II jets, show Γ 10 at scales of tens of kpc and beyond (e.g., Sambruna et al. 2002;Siemiginowska et al. 2002;Tavecchio et al. 2004;Harris & Krawczynski 2006). The FR dichotomy likely reflects a combination of jet power and ambient density (Bromberg & Tchekhovskoy 2016), and our fiducial jet propagating in a homogeneous density profile may resemble an FR-I jet more than an FR-II one, with a small R jet /H jet ratio, where H jet is the extent of the jet. Yet, in this work we consider a broad range of luminosities that should span the higherluminosity parameter space. More precisely, we consider the two following cases.
Case I: FSRQ-Like Powerful Jetted AGNs -As a benchmark for a quite powerful jet of limited (∼ 1kpc) extent, we consider isotropic-equivalent bolometric luminosity L bol = 10 48 erg s −1 with an opening angle θ j ∼ 18 • 1 , R jet = 15pc (H jet ∼ 1kpc), and B 0 = 100µG. Such a large magnetic field is routinely inferred in powerful FR-II jets, but should also pertain to the spines of FR-I jets (e.g., Hardcastle et al. 2004;Hawley et al. 2015). This L bol prescription is reminiscent of blazars, including FSRQs and BL Lac objects, that have L bol that go beyond 10 48 erg s −1 (e.g. Ghisellini et al. 2010). This should enhance the effects of photodisintegration and the production of neutrinos considering that this L bol is deemed quite large (Ajello et al. 2013;Tadhunter 2016;Blandford et al. 2019;Mingo et al. 2019) and particles propagate closer to the base of the jet-where most photons are emitted-compared to the expected size of FR-IIs where H jet can reach hundreds of kpc.
Case II: BL-Lac-like Jetted AGNs -A jet with a more moderate bolometric luminosity L bol = 10 45 erg s −1 with θ j ∼ 18 • , R jet = 1pc, and B 0 = 1.5mG is assumed here. Just as in Case I, the strong magnetic field prescription serves only to study the effect of photodisintegration on UHECRs and the production of astrophysical UHE neutrinos. We choose to set the jet radius to 1pc to further increase the probability of particle interactions as the bulk of the photon field energy is emitted at the base (See §2.4.4 for a discussion on the distance dependence), by setting H jet to the smallest FR-I scales (Hawley et al. 2015).

UHECR injection spectrum
1 This is important to note considering that the true jet power is L bol θ 2 j While the spectrum of UHECRs detected at Earth is measured to be ∝ E −2.6 , the actual spectrum injected by their sources is not well constrained because of the uncertainties in the cosmological distribution of sources and in adiabatic and inelastic losses. While several authors have considered a spectrum of UHECRs injected into intergalactic space, E −q , with q = 2 (e.g., Waxman 1995; Katz et al. 2009;Aloisio et al. 2011)-which may be valid for protons-more recent Auger data favor harder spectra with 1 ≤ q ≤ 1.6 (e.g., Aloisio et al. 2011;Gaisser et al. 2013;Aloisio et al. 2014;Taylor et al. 2015;Aab et al. 2017a;Jiang et al. 2021) to explain the observed heavy chemical composition.
In this paper, we do not account for propagation effects, and we bracket our ignorance of actual UHECR spectrum by considering injection slopes 1 ≤ q ≤ 2. In the espresso framework, the injection spectrum turns out to be flatter than the spectrum of the CR seeds, which should be a power law ∝ E −qse , with q se ∼ 2 − 2.7 (Caprioli 2015), because reacceleration tends to push particles close to the jet's Hillas limit (MC19, MC21). It is worth mentioning that q se could approach values of ∼ 2.7 if we consider relatively short jets like FR-I jets that extend to 10 kpc, where CRs within the CR halo could be reaccelerated (Kimura et al. 2018). A systematic study of espresso acceleration in different kinds of AGN jets is ongoing, but in general we find that q − q se ∼ 0.5 − 1, is consistent with the flatter spectra required to explain Auger data (e.g., Aloisio et al. 2014;Taylor et al. 2015). Similar trends are observed in shear acceleration mechanism, for which it was shown that Auger data can be quantitatively fitted (Kimura et al. 2018).
In general, the spectrum of secondary particles and UHE neutrinos is a function of q; therefore we show results for different values of q se . In the remainder of this paper, we fix the value of the UHECR injection spectrum q such that 1 ≤ q ≤ 2.

Effects of energy losses on UHECR spectra
Let us start the discussion of our main findings by assessing the role of losses on the spectra of reaccelerated UHECRs. As discussed above, the parameters for Case I are chosen in order to maximize the potential losses for UHECRs: a powerful, yet compact, source would in fact force particles to propagate closer to the AGN and hence be exposed to the bulk of its non-thermal emission (see, e.g., Dermer 2007;Murase et al. 2012). The left panel of Figure 4 shows the average atomic mass and spectrum of UHECRs for Case I, with the contribution of different chemical species and for values q =1, 1.6, and 2. It is worth noting that the UHECR spectrum is not significantly affected by photodisintegration, even with a prescription that may magnify its effects because espresso acceleration mainly occurs at scales reaching kiloparsec levels, where the photon density becomes lower than that in the blazar region. Figure 2 indicates that the optical depth at ∼ 10 20 eV is ∼ H jet t −1 Aγ /c ∼ 10 < A, implying that the effective optical depth-taking into account the inelasticity (Murase & Beacom 2010)-is not much larger than unity. Hence, a significant fraction of the heavy nuclei survives losses because seed reacceleration occurs throughout the jet extent, and not just in the blazar region. The light component cuts off at a few times 10 18 eV, and the overall spectrum gets heavier with increasing energy, consistent with Auger observations (Aab et al. 2014a;Aab et al. 2017b;Yushkov 2019). The right panel of Figure 4 breaks down the contribution of each chemical specie for q = 2 (solid lines, color coded) and also shows the total spectrum for q = 1 (dotted line). Unger et al. (2015) suggested that heavy ions may be photodisintegrated if the UHECR confinement time were increased around sources due to the presence of magnetic irregularities. In this picture the secondary nuclei originating in these regions would account for the light composition observed below 10 18 eV. Here we observe a similar phenomenology, in the sense that particle scattering in the cocoon produces a light-element bump of reaccelerated secondary protons around 10 18 eV, provided that the seed spectrum is sufficiently flat. This might correspond to the so-called EeV component that is often invoked to fit the low-energy section of the UHECR spectra (e.g, Gaisser et al. 2013;Aloisio et al. 2014). Note that the position of this bump does not depend on the assumed SGS, in that it corresponds to the Hillas limit for protons; also, increased scattering does not significantly increase the confinement time of particles because particles are more likely to escape sideways from the jet as SGS increases (see MC21).

Neutrinos from interactions of UHECRs inside their sources 3.2.1. UHE neutrino flux expected from a given jetted AGN
The black and red curves in Figure 5 show the expected UHECR and neutrino spectra for two different AGN jets (Case I and II, left and right panel, respectively) and for q = 1, 2 (dashed, solid lines), as examples of spectral slopes. Extrapolations to other q's are straightforward based on the results in Appendix D.
Two main trends arise, as expected: 1) the neutrino flux increases with increasing jet luminosity (compare left and right panels in Figure 5); 2) the highest-energy neutrinos are produced via photomeson production interactions (thin blue curves), while at lower energies ( 10 PeV), neutrinos from the neutron decay of photodisintegrated UHECRs dominate the flux (thin yellow curves). Finally, we point out that all the neutrinos produced here are from pγ interactions; pp collisions are negligible and not visible in the plots. These results do not depend on the level of assumed SGS, since most of the HE neutrinos are produced by powerful AGN jets that can be relativistic even at kiloparsec scales, for which CRs are mainly accelerated via the espresso mechanism (MC21). We also note that the differences in the total CR spectra between the two panels is quite small for q=2 because the contribution of secondary protons is not as pronounced.
The red lines in Figure 5 are calculated based on the propagation of test-particles in our fiducial jet, but it is instructive to also calculate the neutrino spectrum resulting from a much simpler approach that ignores the actual jet structure. Such an analytical one-zone model is fully described in Appendix E; in a nutshell, a neutrino flux can be estimated via a rescaling of the UHECR flux based on an effective optical depth f pγ for pγ interactions, with f pγ given by: where H jet is the extent of the jet, typ ≈ 0.5¯ res /γ is the most probable photon energy, γ is the particle Lorentz factor, L is the differential luminosity of the jet at , ξ is the average energy fraction lost to the pion, and E ν = αE (see Appendix E and Table 1 for more details). We can then introduce f ν (corresponding to red curves in Figure 5), which can be thought of as an effective optical depth for UHECR interactions. f ν globally captures the order of magnitude of the full kinetic approach and can be used to quickly estimate the expected neutrino flux from a given AGN that is active as an UHECR source.

Expected flux of UHE source neutrinos
10 5 10 6 10 7 10 8 10 9 10 10 10 11 10 12 log 10 (E/GeV)  We move now to estimating the overall flux of UHE neutrinos produced by a realistic distribution of the AGN jet luminosity. The energy flux in neutrinos that comes from the convolution over the cosmological distribution of their sources can be expressed as (e.g., Murase et al. 2014;Ahlers & Halzen 2014): where ρ(z) is the number density of sources, L is the AGN jets' luminosity, H(z) = H 0 (1 + z) 3 Ω M + Ω Λ is the Hubble parameter, with Ω M ≈ 0.3 and Ω Λ ≈ 0.7 for standard ΛCDM cosmology, and H 0 = 70km s −1 Mpc −1 is the Hubble constant. The E 2 ν J ν [(1 + z)E ν , L] term is the average neutrino source luminosity.
In this paper, we consider a scenario in which all the UHECRs are produced in environments with intense photon fields, which should get us closer to an upper limit on the possible neutrino flux. Such a flux must be anchored to some expected UHECR luminosity, hence we scale the ρ(z)E 2 ν J ν term with EdQ UHECR /dE(E = 10 19 eV) (see §2.6), since most CR interactions leading to neutrinos occur at 10 19 eV. See Murase & Beacom (2010) for upper limits on the diffuse neutrino flux from the sources of UHECR heavy nuclei.
Both the photon and the UHECR luminosities for an AGN jet should scale with its bulk power (also see MC19). Ghisellini et al. (2009) found that there may be hints that could relate the jet bulk power to accretion disk luminosity; this study was later complemented by findings from Ghisellini et al. (2014) where it was asserted that there is a clear correlation between the accretion luminosity and γ-ray jet luminosity. The UHECR injection rate at any redshift can then be written as a function of the local one such that: (4) where ζ z is a cosmological evolution factor defined in equation F23 of Appendix F and L x the X-ray luminosity of the AGN sources. The prescribed photon fields in our framework are related to L x such that L bol /L x = 10 4.21 where the constant of proportionality is obtained by modeling the γ-ray luminosity function through the observed X-ray luminosity (see, e.g., Ueda et al. 2003;Inoue & Totani 2009;Harding & Abazajian 2012). This is broadly consistent with the more recent results on the blazar luminosity function (Ajello et al. 2012(Ajello et al. , 2014.
Assuming that there are different classes of jetted AGNs, we can then rewrite Equation 3 as: where the weights w i of the different classes are defined as the relative X-ray injection in Mpc −3 of each AGN type: such that dψ dLx (L x , z) is the z-dependent X-ray luminosity function of the AGN sources per luminosity per comoving volume (see MID14 for more details).
In this study, we consider contributions from two types of AGN jets with two different isotropic-equivalent bolometric luminosities: i) FSRQ-like powerful jetted AGNs with L bol = L 48 ∼ 10 48 erg s −1 (Case I), and ii) BL-Laclike jetted AGNs with L bol = L 45 ∼ 10 45 erg s −1 (Case II). We note that the associated L x dψ dLx ∝ ∼ L −2 x , so the relative contribution of each AGN type in this case is roughly the same, i.e., w 45 ∼ w 48 and the total energy injected per unit time in UHECRs is comparable in our case (see discussion below). Note that they could, in principle, accelerate particles to different maximum energies (see the discussion in MC21). We finally obtain: × w 45 f ν (L 45 , E, q)ζ z (10 −4.21 L 45 , z) +w 48 f ν (L 48 , E, q)ζ z (10 −4.21 L 48 , z) (7) Figure 6 shows the resulting UHE source neutrino fluxes (black lines) based on Equation 7 and expected cosmogenic neutrino fluxes (green and blue lines). The green line is the expected cosmogenic neutrino flux from radio galaxies assuming an ion composition consistent with ours (taken from Zhang & Murase 2019, where the shear acceleration scenario following Kimura et al. (2018) is considered). The blue bands, on the other hand, show the expected cosmogenic flux for an AGN model that is more proton-rich based on models that fit Auger's spectral features with different confidence levels . The pink and brown dotted lines show nucleus-survival bounds on the diffuse neutrino flux from sources of UHECR heavy nuclei (Murase & Beacom 2010), which is consisten with the intuition that neutrino production is inefficient in photon-poor environments allowing the survival of nuclei. The yellow data points show the fluxes of the high-energy neutrino events (Kopper & IceCube Collaboration 2017). Shown are also the CR data from Auger (Aab et al. 2020), Telescope Array (Jui 2016), and KASCADE-GRANDE (Apel et al. 2013) for reference.
A few points are worth noticing.
• If injection spectra are sufficiently flat, the flux of source neutrinos may dominate the expected flux of cosmogenic neutrinos (see also Murase et al. 2014;Rodrigues et al. 2021) 2 . Even for a moderately flat spectrum of q = 2, their relative contribution at ∼ 10 17 eV turns out to be comparable.
• Given the dependence of the source neutrino flux on q and α, neutrinos from the β decay of photodisintegration byproducts might become a nonnegligible fraction of the flux observed by IceCube below a few PeV, although their fluxes are likely to be lower than those from the photomeson production (see also Murase & Beacom 2010).
• The neutrino fluxes that we obtain for q=2 are lower than the nucleus-survival bounds on the diffuse neutrino flux stemming from UHECRs with A ≥ 4.

Upper limits of the UHE neutrino flux in the reacceleration sceanrio
The results shown in Figure 6 should be regarded as close to upper limits of the neutrino flux from AGN jets since we assumed that all UHECRs are produced in jetted AGNs with considerable isotropic-equivalent bolometric luminosities (between 10 45 and 10 48 erg s −1 ).
We remark that: • The jet simulations that we consider have a limited extent-consistent with an FR-I jet-with a bolometric luminosity that reaches that of the most powerful blazars. This further maximizes the effects of the photon fields that are generated close to the base of the jet.
• Relaxing any of our assumptions would reduce the contribution of the non-thermal photon background to pγ collisions; the CMB, BLR, and dusty torus contributions alone would provide a source neutrino flux that would be a few orders of magnitude smaller (see left panel of Figure 2), subdominant with respect to both the cosmogenic neutrino flux and the astrophysical neutrino flux in the Ice-Cube band.  Murase & Beacom (2010) with EdQ UHECR /dE(z = 0, q = 2) = 0.6 × 10 43 erg Mpc −3 yr −1 and ξz = 7.2. IceCube neutrino data, along with UHECR data from Auger, KASCADE, and TA are also included for reference.
• Our results are sensitive to the luminosity function as that directly affects UHECR normalization, the effect of photodisintegration, and the neutrino flux prediction. For example, lower-luminosity AGNs (and/or other classes of sources) may give more contributions to the observed UHECR flux, e.g., via the shear acceleration mechanism, because their luminosity density is larger (Ajello et al. 2014). Then the expected neutrino flux would be reduced because the expected neutrino power scales with the bolometric luminosity as shown in Figure 5. Our neutrino flux should still be regarded as an upper limit, since including lowerenergy AGNs should not result in a higher neutrino flux.
The results of Figure 5 and 6 lead to two very general considerations. On one hand, neutrinos produced by the β decay of secondary neutrons should be accounted for when calculating the flux of neutrinos in the PeV range. On the other hand, in order to explain the whole Ice-Cube flux with neutrinos from UHECR sources, photo- the Case I prescription to trace the regions where neutrinos are most likely to be produced. Neutrinos are preferentially produced close to the base of jet within ∼ 800pc, but a non-negligible fraction of neutrinos is produced at larger distances. Lower Panel: Distribution of the final angles of escaping neutrinos where µ f is the cosine of the angle between the particle momentum and the jet axis. We can see that neutrinos escape the jet "quasi-isotropically".
disintegration would need to happen at a much higher rate than what is estimated here, which would be inconsistent with the presence of heavy elements at the highest energies. In other words, the optical depth f ν required to produce a sizable source neutrino flux would necessarily lead to the complete photodisintegration of heavy UHECR nuclei; this result is general and independent of the UHECR source or acceleration mechanism. Overall, better UHECR chemical composition measurements and constraints on the UHECR injection slopes may pose more stringent limits on the expected contribution of neutron-decay neutrinos to the observed Ice-Cube flux.

Sites of production and angular distribution of escaping neutrinos
In our simulations we keep track of where neutrinos are produced and in which direction they escape. The top panel of Figure 7 shows that a considerable frac-tion of neutrinos is produced at relatively large distances (∼ 1 kpc) away from the base of the jet as a result of two competing effects: on one hand, the photon field intensity declines as D −2 , and on the other hand, the non-thermal emission cone affects a greater volume as we move away from the base of the jet. We observe similar trends for all of our prescriptions, so we argue that the photodisintegration of heavy UHECR nuclei and neutrino production mainly occur at intermediate distances, not too close to the blazar region but well before the jet's end.
The bottom panel of Figure 7, instead, shows the distribution of the cosines of the final angles of flight of the neutrinos produced through pγ interactions and neutron decay for Case I. Neutrinos are released quasiisotropically, essentially because UHECRs are efficiently isotropized in the cocoon (see also MC19 and MC21). This is important for neutrino astronomy because any correlation between neutrino directions of arrival and the AGN population is intrinsically connected to how neutrinos are released from their sources (see also the discussion in Caprioli 2018). If neutrinos are strongly beamed along the jet axis, then they would preferentially be associated with AGN jets that point towards us such as BL Lac objects and FSRQs (i.e., blazars). On the other hand, if the emission were less beamed, we might expect all radio-loud AGNs to contribute to the flux, generally producing a more isotropic signal. Our results suggest that we should receive neutrinos from non-blazar AGNs, too, which makes it harder to assess AGNs as UHE neutrino sources on a statistical basis. This prediction is different from the neutrino emission model postulated by Rodrigues et al. (2021). This, however, does not preclude the association of individual events with given sources.

CONCLUSIONS
In this paper, we studied the effects of photodisintegration on UHECRs accelerated in powerful AGN jets and estimate the ensuing flux of high-energy neutrinos. Particles are propagated in a mechanism-agnostic way and we find that UHECRs are generically accelerated via the espresso reacceleration of galactic-like CRs (Caprioli 2015)-independently of the jet morphology (MC19) and on the details of particle transport (MC21)-in promising "relativistic" jets, especially for powerful jetted AGN such as FR-II galaxies.
We used a bottom-up approach in which test particles (CR seeds) are propagated in a fiducial 3D MHD simulation of an ultra-relativistic jet, augmented with realistic modeling (à la MID14) of the photon fields responsible for UHECR losses and hence for neutrino production via different channels (see Table 1). We considered different prescriptions for the AGN size and luminosity ( §2.5) to maximize the interaction rate and, hence, give an optimistic estimate on the expected source neutrino flux, under the additional assumption that powerful jetted AGNs are responsible for the total flux of UHECRs on Earth. We note that relaxing any of these hypotheses would generally lead to a lower neutrino flux.
To summarize, the main findings of this paper are as follows: 1. For typical baryon densities and photon fields, pp interactions are negligible with respect to pγ collisions; moreover, for powerful jetted AGN the most relevant photon field is provided by the nonthermal jet emission, which dominates over the IR dusty torus emission, the optical stellar light, and even the CMB.
2. UHECRs are not heavily affected by photodisintegration, even in the most luminous AGNs; the spectrum of the highest-energy particles gets heavier with increasing energy, as reported by Auger (Aab et al. 2014b). This situation is different from that in the more compact blazar region (Murase et al. 2012;Murase et al. 2014;Rodrigues et al. 2021).
3. In general, the expected neutrino flux scales with the AGN luminosity and with the slope of the seed and the reaccelerated CR spectra ( Figure 5); the rather hard spectra (q 1.5) required to explain Auger data (see, e.g., Aloisio et al. 2014;Taylor 2014) would maximize the neutrino yield ≥ 10 17 eV with respect to softer UHECR spectra. 4. For the most optimistic scenarios (flat UHECR injection spectra), UHE neutrinos produced inside AGN jets may dominate, or at least be comparable to, the expected flux of cosmogenic neutrinos produced during UHECR propagation across the universe ( Figure 6).

Even if
AGNs sustained the totality of the UHECR luminosity of the universe, their steady neutrino emission could not account for the entire IceCube astrophysical neutrino flux ( Figure 6). However, a non-negligible fraction of the source neutrinos in the IceCube band may come from the β-decay of photodisintegration byproducts, rather than from photomeson production.
6. In the reacceleration scenario, neutrinos should be released quasi-isotropically which strongly suggests that non-blazar jets that correspond to radio galaxies should contribute to the UHE astrophysical neutrino flux. This prediction is different from models of beamed neutrino emission (e.g., Rodrigues et al. 2021).
Note that CR acceleration and neutrino production may occur in different regions. Large-scale jets eventually become subrelativistic, where UHECRs may be accelerated by the stochastic acceleration mechanism (Kimura et al. 2018;Rieger 2019). However, the photon density is so low that neutrino production there was shown to be inefficient (Zhang & Murase 2019). On the other hand, AGNs with large-scale jets are embedded in galaxy groups and clusters, where the cosmic rays escaping from the jets can be confined in the intracluster medium for a cosmological time, in which high-energy neutrinos are produced predominantly via pp interactions (Fang & Murase 2018). We defer to a future study the estimate of the neutrino flux expected from nearby sources, such as M87 and Centaurus A, which may be detectable by current/future experiments even if such moderately powerful AGNs were not contributing substantially to the observed flux of UHECRs.
We thank the anonymous referee for providing insightful comments that helped us to correct the neutrino flux from neutron decay in Figure 6. Simulations were performed on computational resources provided by the University of Chicago Research Computing Center. We thank B. Theodore Zhang for providing us with tables of the cosmogenic neutrino flux. This work of D.C. was partially supported by NSF through grant No. PHY-2010240, while the work of K.M. was supported by the NSF Grants No. AST-1908689, No. AST-2108466 and No. AST-2108467, and KAKENHI No. 20H01901 and No. 20H05852.

APPENDIX
A. SPECTRA OF NUCLEI We propagate particles of given rigidity R, which for different species with different atomic charge Z corresponds to an energy E = RZ. The normalizations of the different ion species are chosen according to the abundance ratios at 10 12 eV such that K = K se /K H ∼ [1, 0.46, 0.30, 0.14] for He, CNO, MgAlSi and Fe respectively. Hence: With f (R) = f 0 R −q where q is the spectral slope. Then: where E i is the ion energy and R is the rigidity.

B. PHOTOMESON PRODUCTION INTERACTIONS
We introduce t pγ as the cooling time of a proton with energy E p due to photomeson interactions, i.e.: We follow Stecker (1968) and for an isotropic photon field we have: where is the photon energy in the black hole frame, its energy in the proton frame, ξ is the inelasticity, and th ≈ 0.15 GeV is the threshold energy in the proton frame. We integrate over = γ p (1 − β cos θ), where θ is the angle between the particle momentum vectors in the black hole frame and β ≈ 1. For 0 ≤ ≤ 2γ p , we can introduce the effective cross section σ eff = 70 µb  with σ pγ (¯ ) ∼ σ eff H(¯ −¯ th ) and obtain: where d p expresses the spatial dependence of n γ . When considering interactions between the photons and a nucleus of atomic number A and energy E A , we assume σ eff,A = Aσ eff and ξ A = ξ/A. Based on equation B5, the cooling time for nuclei t −1 N γ eventually reads: Non-thermal emission in jets is typically beamed; hence, the angle θ between the proton momentum and the target photon is fixed and we can express t −1 pγ (γ p ) for an individual proton interacting with a beamed photon field as: such that again = γ p (1 − β cos θ).
For isotropic photon spectra, we can simplify equation C9 by posing σ Aγ ∼ σ G ∆ G δ( − G ). Equation C8 becomes: We assume that G −¯ th > 0 where¯ th is the threshold energy (10 MeV for photodisintegration) and G is the energy at which the cross section peaks. Also, Ξ(d p ) = where E is the energy of the primary particle, E ν = αE is the typical energy of the secondary neutrino, and α can be read in the last column of Table 1 for the different processes that yield neutrinos. While computing the expected neutrino flux, we need to account for a factor of 1/2 that comes from the charged to neutral pion ratio, and another factor of 3/4 stemming from pion and muon decay. Then we have: Since UHECRs up the jet's Hillas limit are espresso-accelerated, we can estimate their propagation time t prop to be close to ballistic such that t prop ∼ Hjet c . The effective optical depth f pγ can then be estimated as: where typ is the most likely photon energy, and typ depends on the considered particle energy such that typ ∼ 0.5¯ res /γ where¯ res ∼ 0.34GeV is the resonance energy and γ = E/m A c 2 is the Lorentz factor of a particle with energy E. Note that n γ (E) is energy dependent. For example, L ∝ 2−β (with β 1), we have n γ (E) ∝ E β−1 . Finally we get an estimate for E 2 ν F UHEν (E ν ) such that: In the example in the left panel of Figure 4, L bol = 10 48 erg s −1 , H jet ∼ 1kpc, and γ > 10 7 so m ∼ 15eV in the black hole frame; the values of α, σ eff , and ξ are discussed more in detail in Appendix B. such that z c = 1.9 and L a = 10 44.6 erg s −1 .