Exploring the Components of Cosmogenic UHECR, Neutrino, and Diffuse Gamma Ray from High-Energy Astrophysical Objects

Abstract

The development of multimessenger astrophysics allows us to probe various background particles from the distant early universe. Up to now, much effort has been made researching the emission and radiation of diverse steady or transient astrophysical sources. We review the potential accelerating, escaping, propagating, and radiation process of high-energy particles under specific circumstances for regular astrophysical sources and briefly discuss the underlying contribution from their emissions to the intensity of ultrahigh-energy cosmic ray, TeV-PeV cosmic neutrino, and the diffuse gamma-ray background, aiming to find a possible common origin.

1. Introduction

The origin of ultrahigh-energy cosmic rays (UHECRs; $E\ge {10}^{18}\mathrm{eV}$), with its long history, has remained a mystery for decades. There are two main approaches to solving this problem: (a) investigating the spectral energy distribution (SED) for different nuclear composition of UHECRs, and (b) exploring the corresponding neutral messengers, such as neutrinos, gravitational waves, and very-high-energy (VHE; $E_{\gamma }>100$ GeV) gamma-ray signals. Given the information on SED for various nuclei and protons, the generation, acceleration, and propagation process will be more precisely clarified. The composition of UHECRs seems to be purer, with one dominant mass composition prevailing within a specific energy range, as confirmed by the collaborative efforts of the Pierre Auger Observatory (Auger) [1,2,3,4]. The sharp feature that is hardening the all-particle cosmic-ray spectrum at $E\sim {10}^{18.7}$ eV is known as an “ankle” structure. Observationally, the depth of shower maximum $X_{\max }$ can reflect the primary mass of the nucleus before initialing the extensive atmospheric shower [5]. The increasing of $X_{\max }$ with energy implies the heavy composition gradually dominates at higher energy ($E\sim$ 40 EeV) [4]. Above ∼10^19.7 eV, the spectrum is suppressed, implying the acceleration capability of sources and the energy loss of propagation in the intergalactic space. The telescope array (TA) also obtains a measurement of SED of UHECRs, and the results agree with the Auger experiment regarding the systematic uncertainty but are significantly biased at energies $E>{10}^{19.5}$ eV [6].

It is known that UHECRs have an extragalactic origin because of the suppression at the highest energy predicted by the Greisen–Zatsepin–Kuzmin (GZK) cutoff [7,8,9], and the large-scale anisotropy $d_{\perp }=6.0_{-0.9}^{+1.0}\%$ depicted by equatorial dipole amplitude above 8 EeV [10,11]. The GZK limit is caused by attenuation from significantly enhancing the cross-section for photopion production, $p+{\gamma }_{bg}\to {\Delta }^{+}\to p+{\pi }^0(n+{\pi }^{+})$ due to $\Delta$-resonance. In the context of interacting with the cosmic microwave background (CMB), the ratio of charged and neutral pions production possibility is about 2:1 and charged pions decay into neutrinos and muons while neutral pions decay into gamma rays, meaning that the interaction between UHECRs and extragalactic background photons or the photon field inner the sources could produce neutrinos and PeV gamma rays.

Assuming an extragalactic origin, UHECRs are deflected by the extragalactic magnetic field (EGMF) when traveling in intergalactic space, which makes it difficult to detect the elusive sources of UHECRs. But the neutral signal (e.g., GeV gamma-rays, neutrinos, and gravitational wave) could be detected directly without losing the location information of sources. The acceleration of UHECRs is always associated with the physical environment, such as relativistic shock or magnetic reconnection, although this remains unclear [12,13,14]. However, the environment at the acceleration site remains consistently dense and relativistic [15,16]. Alternatively, for jets of core collapse supernovae (SN) or gamma-ray bursts (GRBs), the nucleus hardly survives at a close launch radius ($r_0\lesssim {10}^9$ cm) or a high initial radiation luminosity ($L_{\mathrm{rad},0}\gtrsim {10}^{48}$ erg/s) [17]. Therefore, the origin of UHECRs is commonly associated with the production of secondary high-energy neutrinos and gamma rays as a result of the hadronic processes occurring within the sources or during intergalactic propagation (e.g., in the case of magnetars [18]).

Modern astrophysics has witnessed remarkable advancements in observational achievements. The electromagnetic counterpart of the binary neutron star merger GW 170817 was jointly discovered by Laser Interferometer Gravitational-Wave Observatory (LIGO)/Virgo, the Gamma-ray Burst Monitor (GBM) aboard the Fermi satellite, and the INTEGRAL satellite [19,20,21,22]. Subsequent to the kilonova emission observed a day later, the presence of cosmic heavy r-process nuclei was confirmed [23]. Since the initial direct detection of gravitational waves with GW 150914, advanced LIGO and advanced Virgo have reported over 90 gravitational-wave observations spanning distances of thousands of megaparsecs, with most of these events being binary black-hole mergers [24]. An online cumulative set of gravitational-wave transients is available, which is held by LIGO, Virgo, and KAGRA Collaboration (https://catalog.cardiffgravity.org/, last accessed on 15 November 2024). Contrary to hadronic particles and gamma rays, gravitational waves are generated from the birth and dynamics of compact objects, which demand nonaxisymmetric acceleration for numerous materials to be detected. Here, we focus on the relative process of high-energy particles such as UHECRs, neutrinos, and gamma rays in this review. For this reason, gravitational-wave information is rarely mentioned in the subsequent text.

Located in Antarctica, the IceCube Neutrino Observatory enhances our understanding of extreme cosmic accelerators [25,26]. IceCube can detect neutrinos of all flavors by dealing with the Cherenkov signals from charged current interactions, with a position uncertainty of ∼10–15° for cascades and ∼0.5–1° for tracks. High-energy neutrino astronomy has been significantly propelled by the discovery of the 10 TeV–10 PeV astrophysical neutrino flux [27]. The arrival direction of extraterrestrial neutrinos is almost isotropic, which agrees with the extragalactic origin hypothesis [27]. Multimessenger and multi-wavelength observation towards a blazar known as TXS 0506+056 proceeded during a gamma-ray flare reported by the Fermi Large-Area Telescope (LAT) on 28 September 2017. Although the correlation of the neutrino signal has 3$\sigma$ levels during the flare, the constrained luminosity of a muon-neutrino is similar to the gamma-ray luminosity, indicating that jets in blazar could be extragalactic PeV accelerators [28]. The flux of all flavors of cosmic neutrinos is analyzed by the combined likelihood method, which is sensitive to neutrino energies of 10 TeV–10 PeV, and the diffuse flux comes to ${10}^{-7}$$\mathrm{GeV}/{\mathrm{cm}}^2/\mathrm{s}/\mathrm{sr}$ at 30 TeV [29]. An approximately uniform distribution of three neutrino flavors predicted by neutrino oscillation over the cosmological distance is observed [30,31]. If 10 TeV–10 PeV neutrinos are from hadronic $p\gamma$ interaction, a population of cosmic-ray accelerator hidden 1–100 GeV gamma rays is expected [32]. So far, no obvious evidence has been obtained for the correlation between neutrinos and any known types of cosmic objects [33].

Unlike neutrinos of such extremely small interaction cross-sections as ∼10⁻³⁸ $(E_{\nu }/\mathrm{GeV}){\mathrm{cm}}^2$, the VHE gamma rays undergo attenuation by extragalactic background light (EBL). After subtracting the resolved gamma-ray sources and galactic emission, the all-sky gamma-ray emission seems intriguingly isotropic [34]. The LAT experiment measured the spectrum of the isotropic diffuse gamma-ray background at energies of 100 MeV–820 MeV [35], constraining the search for promising gamma-ray populations. Although gamma rays are the most comprehensive messengers for observing astrophysical sources, the origin of the diffuse gamma-ray background remains a subject of debate. It is known, however, that the extragalactic gamma-ray background is predominantly influenced by blazars [36]. Due to $e^{+}{e}^{-}$ production by EBL attenuation, the flux of the VHE gamma ray is suppressed too acutely to be detectable on a high redshift. Moreover, secondary gamma rays resulting from hadronic interactions can induce an electromagnetic (EM) cascade that can persist down to energies of several GeV despite a significant delay time that can span years, implying a potential observational connection.

With the development of ground-based gamma-ray detectors, four gamma-ray bursts (GRBs) have been found with VHE gamma-ray emissions [37,38,39,40], implying the intriguing speculation that GRBs commonly have VHE emission by synchrotron self-Compton (SSC) radiation [41] and can be potential contributors. In the Milky Way, pulsars and pulsar wind nebulae, supernova remnants, and young massive star clusters could be promising PeVatrons, as photons with energy up to 1.4 PeV have been detected by the Large High Altitude Air Shower Observatory (LHAASO) [42]. Those populations are also mysterious gamma-ray and UHECR sources in extragalactic space. The energy scales of three diffuse cosmic particles (gamma-ray, neutrino, and UHECR) are comparable as ∼10⁴⁴–${10}^{45}$ ${\mathrm{ergMpc}}^{-3}{\mathrm{y}}^{-1}$ [43], constraining the injection index or production mechanism of sources to some extent [32,44]. The isotropic energy flux of the high-energy cosmic background of UHECRs, neutrinos, and gamma rays is shown in Figure 1 from data collected by Fermi-LAT, IceCube, and Auger [1,31,35].

If the hadronic interaction and escaping process are disregarded for simplicity, the maximum capability of particle acceleration is directly related to the radiation mechanisms and geometric properties of sources. In some theoretical models, particles are accelerated by an electric field but confined by a magnetic field. The maximum acceleration space scale, defined by Larmor radius $R_{\mathrm{L}}$, should be smaller than the size of the accelerating region, which is known as the Hillas criterion [45]. Considering different acceleration scenarios (stochastic and inductive; see ref [45] for more details) and radiation energy loss, a Hillas plot can be made to decide the capability of acceleration for various sources. Active galaxies cover a large area on the Hillas plot, with only the most energetic ones, such as radio galaxies, quasars, and BL Lac-type objects, being more likely extragalactic UHECR accelerators [46]. The Hillas plot for 100 EeV protons, considering synchrotron and curvature radiation as primary energy-loss mechanisms, is depicted in Figure 2. The star-formation rate provides clues about the unidentified cosmic accelerator among various potential candidates [47]. For a recent and comprehensive review of observational and theoretical studies on star-formation rate, refer to [48].

Deciphering the origins of the high-energy universal background is a massive and demanding undertaking. This review aims to summarize the current research on possible (or even common) origins of cosmogenic UHECRs, high-energy neutrinos, and diffuse gamma rays. The remainder of this review is organized as follows: Section 2 describes the hadronic and leptonic interactions for particle production processes in sources, as well as the acceleration, escaping, and propagation process of the particles. In Section 3, the potential sources, such as active galaxies, supermassive black holes, and GRBs, are listed and discussed. A rough estimate of the flux contribution based on observations from the Gamma-ray Burst Monitor (GBM) and Large-Area Telescope (LAT) is provided in Section 4, along with a summary.

2. Acceleration and High-Energy Radiation Processes

Cosmic rays are accelerated in the extreme physical environment of astrophysical sources. This section introduces the production processes of neutrinos and gamma rays resulting from the hadronic interactions of primary particles. It encompasses not only hadronic processes but also leptonic cooling mechanisms.

2.1. Acceleration Process

The process of converting energy into relativistic jets or outflows remains incomprehensible and subject to ongoing debate. The Blandford–Znajek mechanism is proposed to explain the conversion of energy from a rotating black hole or the surface of an accretion disk in collapsar scenarios [50,51]. For UHECRs, stochastic and inductive mechanisms are commonly hypothesized to achieve nonthermal energies up to ${10}^{20}$ eV. The most common stochastic mechanisms are Fermi first-order [52] and second-order acceleration [16]. The energy distribution of particles after acceleration follows a power-law distribution $dN/dE\propto E^{-\gamma }$, and $\gamma \sim 2.2$ for relativistic shock by Fermi fist-order acceleration or $\gamma \sim 1,\gamma \sim 2$ for other conditions. For the polar caps of rotating magnetized neutron stars, a nucleus can be accelerated through reconnection, inducing a harder spectrum [53,54]. Second-order acceleration also raises a harder spectrum in highly turbulent magnetic fields for active galactic nuclei (AGN) or starburst galaxies [55]. An exponential cutoff is often used to represent the maximum energy due to the Hillas criterion or radiation loss. Maximum energy is not proportional to charge Z due to significant interactions, and the all-particle spectrum is complexed at the acceleration site [56]. The inductive (one-shot or direct) scenario shows that the particles are accelerated by a large-scale and ordered electric field. The maximum acceleration energy is determined by the configuration of the accelerating field, with the energy-loss process being dominated by curvature or synchrotron radiation [57,58,59].

2.2. Hadronic Process

High-energy cosmic rays (CRs) may undergo interaction with the envelope of progenitors for GRBs or supernovas, depending on the acceleration site. After initial loading, hadronuclear and photohadronic interactions can proceed when CRs interact with the circumstantial matter (i.e., accretion disks or interstellar media) or inner photon fields, which makes it difficult to seek the primary injection components of the nuclei. The hadronuclear process is generally dominated by inelastic $pp$ scatterings, with its cross-section increasing with energy. An approximation of the effective optical depth of $pp$ interactions is characterized as [60]:

$$f_{pp}\approx {\kappa }_{pp}{\sigma }_{pp}ct{n}_N,$$

(1)

where $n_N$ is the proton-number density, t is the period of interacting, the $pp$ inelastic cross-section is estimated as ${\sigma }_{pp}\sim 30$ mb, and the inelasticity is ${\kappa }_{pp}\sim 0.5$. The effective optical depth $f_{pp}$ indicates the energy-loss efficiency of the nuclei for hadronuclear interactions. For CRs penetrating the matter in a rectilinear path, the item $ct$ is approximately equal to the whole length $r_l$ of the interaction.

CRs also undergo attenuation at high energy by the inner photons field or cosmic background photons (i.e., CMB, EBL). Considering a general source, the emission-zone radius is r, and the bulk Lorentz factor is $\Gamma$ with a comoving size $r/\Gamma$. The effective optical depth $f_{p\gamma }$ depicts the energy-loss efficiency for protons interacting with photons. If target photons with comoving energy ${\epsilon }_t^{\prime }\approx {\epsilon }_t/\Gamma$ have a energy distribution of $n_{{\epsilon }_t^{\prime }}\propto {{\epsilon }_t^{\prime }}^{-\alpha }$ and the photomeson production is dominated by the $\Delta$-resonance and direct pion production ($\alpha >1$), then $f_{p\gamma }$ can be given by [32]

$$f_{p\gamma }\approx {\eta }_{p\gamma }\left(\alpha \right){\widehat{\sigma }}_{p\gamma }(r/\Gamma )\left({\epsilon }_t^{\prime }{n}_{{\epsilon }_t^{\prime }}\right){|}_{{\epsilon }_t^{\prime }=0.5m_p{c}^2{\overline{\epsilon }}_{\Delta }/{\epsilon }_p^{\prime }},$$

(2)

where the $p\gamma$ attenuation cross-section is ${\widehat{\sigma }}_{p\gamma }\sim 0.7\times {10}^{-28}$${\mathrm{cm}}^2$ [61,62], ${\eta }_{p\gamma }\left(\alpha \right)\approx 2/(1+\alpha )$, and ${\overline{\epsilon }}_{\Delta }\sim$ 0.3 GeV. Specifically, the effective optical depth $f_{p\gamma }$ is model-dependent. As a conservative approach, the estimation $f_{p\gamma }\gtrsim 0.01$ is reasonable since extragalactic CR flux should not exceed the observed all-particle CR flux at 10 PeV [63].

As a result of hadronuclear and photohadronic interactions, the secondary meson as pion and muon can decay into neutrinos and gamma rays; e.g., ${\pi }^{+}\to {\mu }^{+}{\nu }_{\mu }\to e^{+}{\nu }_{\mu }{\overline{\nu }}_{\mu }{\nu }_e$ and ${\pi }^0\to 2\gamma$. The production rate $Q_{{\pi }^{\pm }}$ of a charged pion is related to CR density $Q_N$ that is involved in the interaction. For nuclear density n, inelasticity $\kappa$, cross-section $\sigma$, and an interaction route l, the efficiency factor for producing pions is $f_{\pi }\simeq 1-exp(-\kappa l\sigma n)$. The energy conversion of CR energy into pion energy reflects the production rate as follows:

$$E_{\pi }^2{Q}_{{\pi }^{\pm }}\left(E_{\pi }\right)\approx f_{\pi }{\displaystyle \frac{K_{\pi }}{1+K_{\pi }}}{\left[E_N^2{Q}_N\left(E_N\right)\right]}_{E_N=E_{\pi }/{\kappa }_{\pi }},$$

(3)

where $K_{\pi }\approx$ 1 for $p\gamma$ and $K_{\pi }\approx$ 2 for $pp$ interactions. The average inelasticity per pion can be estimated as ${\kappa }_{\pi }=\kappa /N_{\pi }\simeq$ 0.2 for pion multiplicity $N_{\pi }$. For secondary neutrinos, the production rate needs to be multiplied by a factor of 3/4 because 1/4 of the energy is divided into the production of $e^{\pm }$ in the former decay chain. If the exotic process, such as inverse Compton or electromagnetic cascade, is not considered, the production rate of gamma rays produced by neutral pions is correlated with $Q_{\nu }$ as follows:

$$E_{\gamma }^2{Q}_{\gamma }\left(E_{\gamma }\right)={\displaystyle \frac{4}{3K_{\pi }}}{\left[E_{\nu }^2{Q}_{\nu }\left(E_{\nu }\right)\right]}_{E_{\nu }=E_{\gamma }/2}.$$

(4)

Equation (4) underestimates true gamma-ray flux since the charged pions can lose energy during the decaying, cooling, or adiabatic processes, and gamma-ray production can be enhanced by ambient electrons and positrons.

2.3. Leptonic Process

Relativistic electrons lose energy and emit photons when moving in a magnetic field through synchrotron radiation. The synchrotron is one of the most common emission mechanisms in astrophysics for leptons. For an electron with comoving energy ${\gamma }_e^{\prime }$ and vertical strength of magnetic field $B_{\perp }^{\prime }$, the typical energy $E_{\gamma }^{syn}$ of gamma ray emitted is:

$$E_{\gamma }^{syn}={\displaystyle \frac{3\Gamma ћ{\gamma }_e^{\prime 2}e{B}_{\perp }^{\prime }}{2m_{e}c}}\sim 200\mathrm{MeV}\Gamma {({\gamma }_e^{\prime }/{10}^8)}^2(B_{\perp }^{\prime }/1\mathrm{G}).$$

(5)

The inducing SED of photons is dependent on the energy distribution of injected electrons. For the Fermi acceleration mechanism, the electrons are accelerated to nonthermal energy as a power-law energy distribution, $dN/dE\propto E^{-\gamma }$, resulting in a nonthermal spectrum of emitted photons.

Another common approach to produce gamma rays is through inverse Compton (IC) scattering. In this process, high-energy electrons scatter off-target photons, resulting in the target photons being upscattered to higher frequencies. In both the Thomson regime and Klein–Nashina regime, the energy of the scattered photon can be approximated by:

$$h{\nu }^{IC}\sim min({\gamma }_e^{\prime 2}h{\nu }_0,{\gamma }_e^{\prime }{m}_e{c}^2),$$

(6)

where $h{\nu }_0$ is the energy of target photons. In a deep Klein–Nashina regime such as ${\nu }_0>{\nu }_{KN}\simeq m_e{c}^2/h{\gamma }_e^{\prime }$, the cross-section is highly suppressed due to Klein–Nashina suppression. The synchrotron self-Compton (SSC) process involves electrons upscattering photons that are emitted from the same electron populations through synchrotron radiation, and it is often proposed to explain the high-energy gamma-ray emission in GRBs external shock [64,65].

2.4. The Source-Dependent Escaping Model

Particles are confined in the acceleration region, and only those with energies close to maximal energy can escape from the source. Depending on the configuration of sources, the escape-limited process for AGN, SNRs, and blazars has been derived and applied in the case of diffusive shock acceleration [66,67,68]. In many cases, the escape model favors a harder spectrum than the index of acceleration $s\sim 2$ for diffusive shock acceleration. Regarding arbitrary acceleration processes, the time-integrated spectrum of escaping CRs from the source can be derived from stationary approximation. If we assume a power-law proton production rate with index s and an exponential cut at $p_m\left(\chi \right)$ as a function of dynamical evolution $\chi$ of the accelerator, the time-integrated spectrum in momentum space can be given by [67]:

$$N_{esc}\left(p\right)\propto {\displaystyle \frac{p^{1-s}K\left(p_m^{-1}\left(p\right)\right)}{p_m^{-1}\left(p\right)[d{p}_m/d\chi ]{|}_{\chi =p_m^{-1}\left(p\right)}}},$$

(7)

where $K\left(\chi \right)$ is the local emission rate of proton production at $\chi$. If both $K\left(\chi \right)$ and $p_m\left(\chi \right)$ follow the power-law forms with index $\beta$ and $-\alpha$, respectively, then the time-integrated spectrum also follows the power law form as $N_{esc}\propto p^{-s_{esc}}$, and the index of escaping particles is $s_{esc}=s+\beta /\alpha$. For the Sedov phase of diffuse shock acceleration, the escaping spectral index is equal to the injection index $s_{esc}=s$, as the kinetic energy of CRs is dominated by relativistic ones [69].

2.5. Maximal Acceleration Energy Related to Radiation and Geometry

The maximal energy for alternating groups of elements can be classified into two cases: a rigidity-dependent effect caused by a Peters cycle in the accelerators or the average energy per nucleon caused by photonuclear spallation processes. The maximal energy of CRs satisfies the Hillas criterion and radiation-loss constraints. For a particle with charge $Ze$ and atomic mass A, it must be confined in the acceleration site as [46]:

$$E_{max,H}\simeq ZeRB(v/c)\sim 9\times {10}^{23}Z(R/1\mathrm{kpc})(B/1\mathrm{G})\mathrm{eV}$$

(8)

The energy-loss power due to radiation for particles also cannot exceed the energy acquisition power. For diffusive acceleration case, the scenario $\overrightarrow{v}\Vert \overrightarrow{E}\Vert \overrightarrow{B}$ rarely occurs for large-scale jets [57]. The energy loss is dominated by the synchrotron process, and the maximal energy is given by [46]:

$$E_{max,r}\sim 3\times {10}^{16}{A}^4{Z}^{-4}{(R/1\mathrm{kpc})}^{-1}{(B/1\mathrm{G})}^{-2}\mathrm{eV}.$$

(9)

For inductive acceleration cases, the energy-loss constraint is determined by the synchrotron and curvature radiation depending on how much the condition $\overrightarrow{v}\Vert \overrightarrow{E}\Vert \overrightarrow{B}$ is satisfied. The maximal energy of the inductive acceleration case is [46]:

$$E_{max,r}\sim 2\times {10}^{20}{A}^2{Z}^{-1.5}{(B/1\mathrm{G})}^{-0.5}\mathrm{eV},$$

(10)

for synchrotron-dominated cases (e.g., jets of active galaxies [57]). When the curvature radiation dominates the energy loss, the maximal energy is:

$$E_{max,r}\sim 1\times {10}^{22}A{Z}^{-0.25}{(R/1\mathrm{kpc})}^{0.5}{(B/1\mathrm{G})}^{0.25}\mathrm{eV}.$$

(11)

Generally, a critical magnetic field strength $B_c$ is determined to decide which criterion dominates between the Hillas criterion and the radiation constraint, which is approximated by [46]:

$$B_c\left(R\right)\simeq 3\times {10}^{-3}{A}^{4/3}{Z}^{-5/3}{(R/1\mathrm{kpc})}^{-2/3}\mathrm{G}.$$

(12)

In the regime $B>B_c$, the radiation loss dominates the acceleration process, and the maximal energy is $E_{max,r}$, while in the regime $B<B_c$, the maximal energy is $E_{max,H}$.

2.6. Propagation of UHECRs and EGMF

As mentioned earlier, relevant hadronic interactions cause CRs to lose energy and produce high-energy neutrinos and gamma rays. For UHECR protons and nuclei, photonuclear interaction with CMB and IRB is the dominating energy-loss process. Models and the cosmological evolution of IR/Opt/UV density have been proposed and agree well in the far-infrared band at low redshift [70,71,72]. For protons, the main energy-loss mechanisms include adiabatic loss (<${10}^{18}$ eV), pair production (∼${10}^{18}$ eV), and pion production (∼$7\times {10}^{19}$ eV). However, unlike the proton case, nuclei mainly undergo photodisintegration and a decrease in the Lorentz factor of the nucleus, which is caused by adiabatic expansion and pair production (See Figure 3). The spectrum after propagation is closely related to various cross-sections, such as giant dipole resonance and baryon resonances for nuclei and protons, as well as evolution with redshift of energy and density for the photon background (see reference [73] for a detailed review).

While propagating in the universe, charged particles are thought to be deflected by EGMF, which is a matter of debate. The origin of EGMF is still a mystery, believed to result from the amplification of so-called seed fields in galaxies or in the intracluster medium [74,75]. To determine the large-scale structure and strength of EGMF, magnetohydrodynamical simulations are involved [76,77]. The average deflection of UHECRs by EGMF can be expressed as [77]:

$$\delta \approx 0.9Z\left({\displaystyle \frac{100\mathrm{EeV}}{E}}\right)\sqrt{l_c/1\mathrm{Mpc}}\sqrt{D/10\mathrm{Mpc}}(B/1\mathrm{nG}),$$

(13)

where $l_{\mathrm{c}}$ and B are the coherence length and strength of EGMF, $E/Z$ is the rigidity of particles, and D is the distance to source.

3. Potential Astrophysical Sources

The detection of UHECRs depends on large-area ground-based experiments due to their extremely low event rate. But the variable models of extensive air shower (e.g., Epos-LHC [78] and Sibyll2.3 [79]) introduce quite an amount of uncertainty to the spectra of different components of nuclei, which, in turn, affects the flux of high-energy neutrinos and diffuse gamma rays originating from hadronic processes. There are two kinds of hypotheses to investigate the possible origin of ultra-high-energy (UHE) particles diffused in the universe. The top-down method involves a hypothesis of source distribution and evolution, as well as a population of primary escaping particles. This method propagates the UHECRs in the cosmic photon background field and calculates the predicted spectra, ignoring the underlying intrasource processes. The bottom-up approach is more explicit and popular, which focuses on typical astrophysical sources. The former case is more likely to have a cosmogenic origin, while an astrophysical origin is more suitable for the latter case. This review primarily focuses on certain explicit astrophysical sources, with cosmogenic research also introduced at the end of this section.

3.1. Origins of UHECRs and Neutrinos

To explain the characteristics of the UHECR spectrum, such as the ankle structure and the so-called GZK cutoff, a “pair production dip” model is proposed. This model interprets the ankle structure as the shortening energy-loss length of pair production and pion production with increasing energy [80]. The dip describes the hardening of the spectra at several EeV energies, which is caused by electron–positron production by extragalactic protons interacting with CMB photons: $p+{\gamma }_{CMB}\to p+e^{+}+e^{-}$. Another natural assumption is that the ankle structure is caused by the transition of galactic CRs and extragalactic CRs, which is known as the transition model.

For UHECRs with energy higher than ${10}^{20}$ eV, the average path is always shorter than about 100 Mpc due to the pion production with CMB photons, defining a GZK radius. Therefore, the different cosmological evolution of sources has little impact on the highest-energy part of the observed spectrum. Radio galaxies are always a promising accelerator for those UHECRs, as their lobes can easily satisfy the Hillas criterion if an amount of energy is divided in the turbulent component of a magnetic field [81]. Considering nearby sources, about 10 powerful FRII radio galaxies reside within a few hundred Mpc (e.g., 3C 98, 3C 442A). The FRII radio galaxies are also a promising origin as their average gamma-ray emission is about ${10}^{46}$ erg${\mathrm{Mpc}}^{-3}{\mathrm{y}}^{-1}$, which is two orders of magnitude higher than the required production rate to power UHECRs [82,83]. The emission characteristics of blazar jets demand that maximal acceleration energy is ∼10¹⁹ eV for protons [84]. Therefore, they can only be the ion accelerator to fulfill the composition of the Auger observation.

If blazars accelerate UHECR protons and ions with a significant fraction, some problems of observation, such as the extra spectral component, will be interpreted [85]. The investigation of correlations between the positions of nearby (0–100 Mpc) AGN and the arrival direction of UHECRs shows that UHECRs are more isotropic compared to AGN, and the latter can hardly be the predominate origin of UHECRs within a small smearing angle of ${\theta }_{\mathrm{s}}<10$°, or a certain subclass of AGN within 60–80 Mpc is a more promising origin [86]. As shown in Figure 4, the arrival direction of CRs from AGN is dependent on the two most nearby, AGN Centaurus A and Messier 87, assuming an equal UHECR luminosity for simplicity. Moreover, some AGN flares caused by tidal disruption events (TDEs) can also accelerate CRs to UHE. TDEs often occur in the vicinity of the cores of galaxies when a star passing close to the central black hole is disrupted by tidal force. The relativistic jet formed during a tidal disruption event could accelerate particles to UHE. TDEs can account for the energy-injection rate for UHECR, especially for events like AT 2018hyz with a power density of ∼${10}^{45}$ erg${\mathrm{Mpc}}^{-3}{\mathrm{y}}^{-1}$ [87]. In the scenario where a white dwarf is disrupted by an intermediate-mass black hole, TDEs could account for the UHECR spectrum with dominant compositions being Ne and Mg, assuming a fixed escaping particle spectral index of $s_{esc}=1$ in a constant circumstellar medium [88].

Regarding GRBs as the origin of UHECRs, the energy-injection rate is an order of magnitude lower than the required flux. If UHECRs originate from the internal shock of GRBs, with a proton hypothesis, the dip model is highly disfavored because the extraordinary baryonic loading $f_{\mathrm{e}}^{-1}>{10}^5$ under the constraint of neutrino bound [89]. As a special population of GRBs, low-luminosity GRBs have lower luminosity, such as $L_{\gamma }<{10}^{49}$ erg/s. They are more transparent to heavy nuclei for interactions with gamma rays within the source. They are convincing contributors of UHECRs within the normal parameters of space compared with high-luminosity GRBs [90,91]. For sources such as star-forming galaxies and misaligned active galactic nuclei, assuming a GRB-like cosmological evolution, the ankle feature can be interpreted as a transition of extragalactic and galactic components [92]. Soft proton spectral indices $\beta \sim 2.4$–2.5 and a strong cosmological evolution are needed to conform to the KASCADE-Grande data, as is shown in Figure 5.

The observation of a heavy-dominated composition at high energy by the Pierre Auger Observatory conflicts with scenarios involving AGN and most GRBs as pure proton sources. However, this composition is compatible with neutron stars such as young pulsars or magnetars. Another plausible explanation is that heavy nuclei could be synthesized within the source environment of GRBs, as well as in transrelativistic supernovae, newborn pulsars, and magnetars. Newborn pulsars, with their energy budget, metal-rich surface, and dense configuration, can be potential UHECR accelerators. For the wind-acceleration coefficient $\eta =0.3$, the galactic magnetic field coherence length $l_{\mathrm{c}}=20$ pc, and the magnetic halo height $H=2$ kpc, fast-spinning pulsars could account for the UHECR spectrum and the composition prediction of Auger observation [54].

Starburst galaxies typically have a star-formation rate that is 1000 times greater than that of normal galaxies, leading to the efficient ejection of gas. The superwinds produced by starburst galaxies create metal-enriched plasma bubbles in the regions of active star formation. Ultra-high-energy cosmic rays can be accelerated in the collisionless plasma shock waves formed by the expansion of these bubbles. The cosmic rays produced can be confined within the dense and strong magnetic regions for millions of years, leading to the production of neutrinos from hadronic cascades resulting from multiple $pp$ interactions. The frequent supernovae in starburst galaxies can give rise to a large number of newborn pulsars that have the potential to act as accelerators of UHECRs. The Auger collaboration has presented a possible correlation between nearby starburst galaxies and UHECRs with energies greater than ${10}^{19.6}$ eV at approximately 4$\sigma$ significance [93]. If these starburst galaxies have a cosmic-ray spectrum that is harder than $E^{-2.2}$, they could contribute a significant fraction of the neutrino flux above 100 TeV.

In the study of the astrophysical origin of high-energy neutrinos, a stacking analysis is often employed to search for unknown sources that correlate with EM observation. However, the statistical significance of these correlations can be diminished by the look-elsewhere effect. Among the sources considered, only two show promise as potential candidates: active galaxies TXS 0506+056 and NGC 1068. For a population of sources with proton spectral indices around $\beta \sim 2.4$, the predicted flux of diffuse high-energy neutrinos can be estimated, as illustrated in Figure 6.

The neutrinos produced by the interactions of UHECRs with matter or photons at their sources carry approximately a 5% fraction of the primary energy. Therefore, the detection of these extremely high-energy neutrinos can help constrain the potential sources of UHECRs. However, despite various theoretical frameworks, the highest-energy neutrinos (>$5\times {10}^{14}$ eV) detected by IceCube suggest a slower cosmological evolution compared to the star-formation rate for astrophysical origins [47]. This finding presents a tension with the conclusions drawn in the study by reference [92].

3.2. Origin of Diffuse Gamma Rays

The study of diffuse gamma rays presents a more complex scenario compared to neutrinos, as gamma rays can be emitted both through hadronic processes and through leptonic nonthermal processes. When considering the leptonic–hadronic process, a variety of astrophysical objects such as AGN, blazars, star-forming galaxies, millisecond pulsars, GRBs, radio galaxies, and even exotic processes like dark matter annihilation, can contribute to the diffuse gamma-ray emission [94,95,96,97,98]. Current research suggests that the extragalactic gamma-ray background, including point sources, is likely dominated by blazars, as shown in the left panel of Figure 7. Flat-spectrum radio quasars, core-dominated radio galaxies, and GRBs contribute only a small fraction to the isotropic gamma-ray background [99,100,101,102], with star-forming galaxies potentially representing a significant portion (approximately 61%) [103]. Recent studies indicate that star-forming galaxies could account for the entire intensity of the isotropic diffuse gamma rays, as illustrated in the right panel of Figure 7.

3.3. Multimessenger Correlation

The investigation into whether high-energy particles share a common origin is being pursued through a correlation study involving IceCube neutrinos, Auger, and TA UHECRs [106]. However, the lack of significance exceeding 3.3$\sigma$ may be attributed to the deflection caused by the EGMF. The ANTARES collaboration has also conducted a search for an excess of neutrinos in the vicinity of UHECR directions and for UHECRs near the directions of the highest-energy neutrinos, but no significant findings have been reported [107]. Despite GRBs and X-ray binaries being considered to be potential sources of high-energy neutrinos, no correlations with statistically significant excess have been identified [108,109]. Two TDEs, AT 2019dsg associated with IceCube-191001A and AT 2019fdr associated with IceCube-200530A, are luminous TDEs with luminosities in the range of ${10}^{44}$–${10}^{45}$ erg/s, suggesting that TDEs could potentially contribute to the all-sky neutrino flux [110,111]. If we consider the diffuse gamma-ray background to be produced by $pp$ interactions or $p\gamma$ interactions from unresolved blazars, the relationship between the flux of diffuse gamma rays and diffuse neutrinos can be established through Equations (3) and (4). In astrophysical scenarios, a large two-photon annihilation optical depth is required, and populations such as choked GRB jets or supermassive black-hole cores may be more promising origins of extragalactic neutrinos [32].

3.4. Cosmogenic Investigation of Diffuse Background

The top-down approach typically assumes cosmic-ray reservoirs, such as galaxies and clusters, with their spatial matter distribution, injection index for compositions, the maximal composition-dependent energy reflecting the capability of source acceleration, and the redshift evolution source related to the star-formation rate. In the case of galaxies and clusters, cosmic rays can be contained within these reservoirs and accelerated by shocks associated with large-scale structure formation or embedded AGN jets. Above 100 TeV, neutrinos align with the predictions of galaxy clusters and groups [112,113]. Using this three-parameter model, the dip-proton model faces significant challenges at a confidence level of more than 95%, independent of source compositions. This is demonstrated by the neutrino excess beyond the upper limit of the IceCube diffuse neutrino flux [114], as illustrated in Figure 8.

4. Rough Estimation

As mentioned in Section 2 and Section 3, under the assumption of a photohadronic origin, the energy flux of diffuse gamma rays is linked to the energy flux of diffuse neutrinos by the relationship $E_{\gamma }=2E_{\nu }$. To provide a comprehensive overview of their interconnection, we examined the energy flux of Fermi-LAT gamma-ray sources. Applying Bayesian statistics with an association probability cut of $P>0.8$ and a detection significance threshold of $\sigma >4$, and excluding galactic sources with $\left|b\right|>20$°, we identified 2761 gamma-ray sources from the Fermi Large-Area Telescope Fourth Source Catalog Data Release 4 (4FGL-DR4) [115]. Among these gamma-ray sources, 1156 were classified as BL Lac-type blazars, 648 as flat-spectrum radio quasars (FSRQs), and 790 as active galaxies of uncertain type. To simplify matters from an experimental perspective, we summed the total energy flux in the 0.1–100 GeV range from LAT power-law fits for each population. For the isotropic diffuse gamma-ray background (IGRB), the extragalactic gamma-ray background (EGB), diffuse neutrinos, and UHECRs, the total energy flux of the IGRB was calculated as:

$${〈E^2{\displaystyle \frac{dN}{dE}}〉}_{\mathrm{IGRB}}={\int }_{0.1\mathrm{GeV}}^{100\mathrm{GeV}}E{\displaystyle \frac{dN}{dE}}dE,$$

(14)

while EGB, neutrinos, and UHECRs are calculated in the same way in their corresponding energy bands.

The summed energy flux of sources from 4FGL-DR4 is depicted in Figure 9, alongside the total energy flux and spectrum of EGB, IGRB, neutrinos, and UHECRs. The majority of the gamma-ray flux from Fermi-LAT sources originates from BL Lac blazars and flat-spectrum radio quasars, contributing a combined flux of $6.8 \times {10}^{-10}$ erg$/{\mathrm{cm}}^2/\mathrm{s}/\mathrm{sr}$ and $7.8 \times {10}^{-10}$ erg$/{\mathrm{cm}}^2/\mathrm{s}/\mathrm{sr}$, respectively, positioning them as primary contributors to the EGB or potential candidates for the IGRB. Among the selected LAT source populations, 81-millisecond pulsars exhibit an energy flux of $7.7 \times {10}^{-11}$ erg$/{\mathrm{cm}}^2/\mathrm{s}/\mathrm{sr}$, which is comparable to the average flux of FSRQs. This suggests that the energy output of millisecond pulsars and FSRQs surpasses that of other populations. Simplifying the consideration by overlooking the details of SED and emission activity, millisecond pulsars and FSRQs emerge as promising candidates for UHECR accelerators in a hadronic scenario.

5. Open Questions and Future Prospects

In this review, we reflect on recent advancements in both observation and theory regarding potential sources of diffuse gamma rays, neutrinos, and UHECRs, as well as their interconnections. Many of these theories rely heavily on models, and concrete observational evidence is often lacking. To identify the most compelling and direct evidence for UHECR accelerators, it is imperative to observe sources with statistical significance. However, due to deflection in EGMF and interactions with CMB, only the lightest composition with the highest energy within the GZK radius can be distinguished from other deflected backgrounds. Furthermore, the composition at the highest energies and its evolution with rigidity is associated with the acceleration mechanism, whether it be the Peters cycle or spallation effects.

The study of UHE neutrinos and gamma rays in astrophysics provides insights into high-energy astrophysical sources and reveals new physics that surpasses the Standard Model at higher energy scales [117]. With the sensitivity of EeV neutrinos reaching ${10}^{-9}$ GeV$/{\mathrm{cm}}^2/\mathrm{s}/\mathrm{sr}$ and considering the evolution of the star-formation rate, the photon fractions of UHE cosmic rays can be indicative of a 10% level predicted by the GZK mechanism. Additionally, in the realm of UHE gamma rays, the interaction length with the CMB extends to tens of megaparsecs as energy increases, rendering the universe more transparent to gamma rays with energies exceeding ${10}^{19}$ eV. These secondary neutral particles play a crucial role in efficiently probing the physics associated with UHE cosmic rays.

UHE multimessenger astronomy necessitates a thorough understanding of extensive air showers. Therefore, large-area ground-based experiments and space experiments capable of detecting fluorescent radiation from extensive air showers are essential. The K-EUSO experiment, led by the JEM-EUSO collaboration and equipped with a Schmidt UV telescope, aims to compare the observations of the Auger and TA experiments in the two hemispheres [118]. On the other hand, the Giant Radio Array for Neutrino Detection (GRAND) is designed to detect cosmic particles with energies surpassing ${10}^8$ GeV [119]. Leveraging radio techniques for air shower detection, GRAND offers unprecedented capabilities for energy reconstruction and $X_{\max }$ measurement. With its extensive exposure, GRAND has the potential to identify approximately 32,000 cosmic-ray events with energies exceeding ${10}^{19.5}$ eV within a five-year timeframe. Some of the challenges pertaining to the energy spectrum and mass composition could be addressed by these future experiments or upgrades. The future detection of cosmogenic gamma rays and EeV neutrinos may pose challenges to certain existing theories in the field.