Testing hadronic and photo-hadronic interactions as responsible for UHECR and neutrino ﬂuxes from Starburst Galaxies

.


Introduction
The discovery of a diffuse stream of cosmic particles, that extends in energy up to ∼ 10 20 eV, an energy range that significantly exceeds every Earth-based accelerator, is one of the most intriguing astrophysics discoveries of the past century. We have been able to investigate its spectral behavior and composition in terms of atomic nuclei thanks to decades of observations. Despite all, one of the most intriguing open problems in contemporary astrophysics is the mystery concerning ultra-high energy cosmic rays (UHECRs), whose nature and origin are still unknown. In order to provide an answer to such a question, an analysis of measurements by the Pierre Auger Observatory [1] was recently updated in [2], where the correlation between source catalogs and UHECRs at the highest energies is investigated. In particular, a strong correlation has been found (4.2σ, predicted to increase to 5σ in 2026) between the arrival directions of UHECRs and the coordinates of the starburst galaxies (SBGs) in the catalogue [3], taking into account the 23 brightest nearby objects with a radio flux greater than 0.3 Jy, even if the contribution from SBGs is summed to an isotropic background; these studies sup- * e-mail:condorelli@ijclab.in2p3.fr port the hypothesis that SBGs could be a potentially interesting class of sources for UHECRs.
One of the most important aspects to investigate would be how the starburst environment affects the UHECR interactions if UHECRs were created in the most active regions of SBGs, known as starburst nuclei (SBNi). As a matter of fact, in various studies [see e.g. [4][5][6][7][8] it was shown that the post-processing of UHECRs through photo-disintegration of CR nuclei in the environment surrounding a hypothetical source can qualitatively explain the UHECR spectrum and composition across the so-called "ankle," the feature corresponding to the flattening of the spectrum near 5 × 10 18 eV [10].
The photo-disintegration process in these models, which will be referred to as "source-propagation models," functions as a high-pass filter, allowing the highest energy cosmic ray nuclei to escape unscattered while the lowest energy ones are disintegrated inside the source region, resulting in a pile-up of nucleons with energy scaling as 1/A, where A is the mass of the nucleus injected in the acceleration region. Particles escaping the source environment are then propagated through the intergalactic medium and finally the obtained diffuse fluxes are compared to the experimental data at Earth.
Hadronic interactions with the source environment have recently been taken into account for a generic source, demonstrating that they can potentially contribute to the same escape effect [8].
Diffuse fluxes of gamma rays (up to TeV energy) and HE neutrinos (up to PeV energy) have been observed respectively by Fermi-LAT [see 11] and IceCube [see 12, and references therein]. Starburst galaxies have already been suggested as potential candidates for such fluxes, however a comprehensive modeling of UHECR interactions in such environment has not yet been thoroughly explored, primarily because there is currently no acceleration model that can inject particles at the highest energies. In fact, while cosmic accelerators such as powerful supernova remnants [see e.g. 13, and references therein] may be able to accelerate particles up to PeV, it is unclear whether or not the EeV range of acceleration may occur in a starburst environment or in the surrounding wind bubbles [see e.g. [14][15][16][17]. Because there is currently no comprehensive theory for particle acceleration, it is important and relevant to conduct a phenomenological study of UHECRs in starburst environments in order to comprehend the fundamental characteristics of this particle population.
In this work, we assume that SBNi are capable of powering UHECRs and, for the first time, we use a sourcepropagation model to derive the UHECR and high-energy neutrino fluxes from these sources. We develop an extension of the Monte Carlo code SimProp [see 18] to simulate the behavior of UHECRs and study the multimessenger implications in terms of associated High Energy (HE) neutrino flux, focusing on the role of the hadronic and photohadronic interactions in these environments.

UHECR interactions in Starburst Nuclei
Starburst galaxies are peculiar astrophysical structures characterized by an extreme star forming activity that can be as high as ∼ 10 − 10 3 M yr −1 [see 19]. A corresponding increase in the rate of supernovae (SNe) (R SN ∼ 0.1 ÷ 1 yr −1 ) is frequently seen as the star formation rate (SFR) increases along with the infrared luminosity. SBGs are particularly effective cosmic ray (CR) factories thanks to their enhanced SNe rate, and connections between the non-thermal brightness and the SFR have been observed [see e.g. 20, and references therein].
The SBNi environment exhibits extreme conditions such as a gas density as high as (or higher than) n ISM ∼ 10 2 cm −3 , magnetic field (B) at the level of ∼ 0.1 − 1 mG and infrared photon density (U RAD ) often higher than 10 3 eV cm −3 . Additionally, a strong wind might be launched under favorable conditions due to the superposition of many SNe and the intensive star formation activity with an estimated speed: (v W ) of about ∼ 10 2 − 10 3 km s −1 .
PeV and sub-PeV CR protons may lose a constant fraction of their energy through proton-proton interactions before being able to escape, primarily because of the advection in the wind, according to the high level of turbulence expected for the SBN environment [21]. Recent investigations proposed that CRs can be additionally accelerated up to ∼ 10 2 PeV at shocks in the wind bubbles [14,17],   Figure 1. Photon spectrum of the prototype SBG, as inspired by [21]: thermal dust modified black bodies (red line) and optical star black body (green line). The black points refer to the measurements from [22] for two different SBG: M82 and NGC253.
whereas in [16] it is argued that in the same conditions energies up to ∼ 10 2 EeV could be achieved. The protonproton (pp) and proton-gamma (pγ) interactions on the SBN diluted photon field will still cause some of these HE particles to lose some of their energy, but their energy will effectively allow them to diffuse away from the starburst environment.
In this study, we focus on the multi-messenger implications of such particle population in terms of UHECRs and HE neutrinos, assuming that CR nuclei are accelerated in the SBN environment up to the highest energies.

Interactions and escape from starburst environment
The low energy photon field of SBGs is complex and is composed of the superposition of a number of thermal components of various wavelengths, from the far infrared (FIR) to the optical (OPT) and ultraviolet (UV) [22]. In particular, the most prominent spectral components are 1) a blackbody associated to the starlight peaking at ε opt 1 eV and 2) a second thermal component peaking at ε IR 10 meV resulting from the reprocessing of the UV light by interstellar dust. We assume a stereotypical lowenergy spectral energy distribution (SED) approximated by these two thermal components. In our investigation, it was decided to adopt a prototype SBG, that is, an SBG with the specifications specified in in Tab. 1 of [23], in order to analyze the environment surrounding the SBN and how it effects UHECRs. The photon spectrum for our prototype SBG is shown in Fig. 1 where it is compared with the spectra of two nearby SBGs: M82 and NGC253 [22]. Using a modified version of SimProp, typical timescales for photo-hadronic and hadronic interactions of CR particles in the SBN are presented in the following.
Under the assumption of a monochromatic photon field of photon density n γ , the typical interaction rate between a relativistic atomic nucleus (A) and a low energy   Figure 2. Timescales for photo-interactions in our prototype SBG for five different injected CR nuclear species as indicated in the legenda (see [23] for details).
photon is approximately τ −1 Aγ cσ Aγ n γ , where σ Aγ represents the cross section of the process. The interaction rate reads as follows when a more realistic photon density is taken into account along with the cross section's dependence on energy: where Γ is the interacting nucleus's Lorentz factor. The primed symbols (such ) relate to quantities in the nucleus rest frame, whereas the unmarked symbols refer to quantities in the laboratory frame. The interaction timescales, which are the inverse of Eq. (1), are displayed in Fig. 2 photo time for various nuclear species, taking into account the reactions of electron-positron pair photo-production, nuclear photo-disintegration, and photo-pion production. The photo-disintegration on the OPT component causes the drop at low energy, while the FIR peak dominates the interaction timescale at higher energies. The other reactions indicated above have a smaller influence than photodisintegration.
Despite having little effect on the extra-galactic medium, spallation processes between CR nuclei and their surrounding gas play a remarkable role in the ISM of SBNi due to the typical densities encountered in active starforming regions. The timescale of the spallation process is given as follows:: where n ISM is the ISM gas density in the SBN environment and σ refers here to the spallation cross section. This process has been implemented in SimProp adopting the most recent hadronic model, Sibyll 2.3d [24], an event generator designed for Monte Carlo simulations of atmospheric cascades at ultra-high energies. The cross section for proton-nucleus pA and pp interactions can be calculated using Sibyll 2.3d, which then enables the computation of the typical timescale for the spallation  process. Additionally, Sibyll 2.3d allows for the estimation of hadronic interactions while accounting for nuclear fragmentation and the rapidity of secondary particles produced at each interaction. The computation of the longitudinal momentum distribution, in particular, is critical for determining secondary particle fluxes.
High energy particles usually remain confined in an astrophysical environment for a short amount of time before escaping. In reality, particles might escape the system through diffusion or advection in a wind. The advection timescale can be expressed specifically as t adv = R/v W , where R denotes the size of the source and v W denotes the wind speed. The diffusion timescale reads: t D = R 2 /D, where D is the CR diffusion coefficient computed in the context of quasi-linear theory and assuming a coherence length l c ∼ 1 pc for the magnetic field. Such an assumption for l c is motivated by the typical scale at which the turbulence is expected to be injected in the SBN [see e.g. 21]. In reality, supernovae are thought to be the source of the turbulence, and ∼ 1 pc is in line with the typical size of a young supernova remnant. The expression of the diffusion coefficient is: D cr 2−δ L l δ−1 c /3, where r L = E/qB is the particle Larmor radius and δ is the spectral slope of the turbulence, E is the energy and q is the charge of the particle while B is the strength of the magnetic field. In particular, we assume δ = 5/3 as prescribed for a Kolmogorov turbulence cascade. Following [25], we also take into account the change in the diffusion regime occurring when r L l c . In this energy range the diffusion coefficient is estimated as D = D 0 (r L /l c ) 2 , where D 0 is the value of the diffusion coefficient computed at the energy E 0 such that r L (E 0 ) = l c . We finally estimate the escape timescale, shown in Fig. 3, as the minimum between the advection and the diffusion time, namely The typical timelines for interactions and escape in the source environment for our prototype SBG are summarized in Fig. 3. Depending on the source characteristics and the CR spectrum at the acceleration site, the interplay between timescales affects the shape of the CR fluxes to be released in extra-galactic space as well as the mass composition. We note that the spallation has the shortest timescale and so dominates the transport in the lowest energy range (E 10 18 eV). The photo-interaction with infrared photons is the dominating mechanism at higher energies E 10 18 eV). Take into account that the diffusion is ballistic in this energy range, and these two processes are in competition.
Spallation produces more secondary particles than photo-hadronic processes, which leads to more photons, neutrinos, and lighter nuclear fragments. The evolution of the nuclear cascade inside the source environment is changed as a result of the mixed composition of propagating nuclei. The mass distributions of nuclei escaping our prototype SBG are shown in Fig. 4 when the injection is assumed to consist entirely of silicon nuclei with energies of 10 17 and 10 19 eV. We particularly concentrate on differentiating the transport effects when there is only one interaction mechanism occurring inside the source: the photohadronic-only scenario (blue) is displayed separately from the spallation-only scenario (red). It can be seen that the photo-hadronic scenario does not efficiently produce intermediate mass (C, N, and O) nuclei, but the spallation scenario easily produces all lighter nuclei.

Implementation of source-interactions in SimProp
We created an extension of an existing Monte Carlo code called SimProp [18] in order to estimate the escaping flux from an SBG. . This software has so far been developed and used in relation to the extra-galactic propagation of UHECRs, such as in [26,27], while in this work, it has been modified to represent both the transport inside the source . All important quantities are considered to be constant within the SBN, where they are propagated, for the purpose of this analysis. UHECRs are considered to be unaffected by the SBG environment once they have left this area and have started to propagate through the extragalactic medium. According to the following assumptions, the propagation inside the source is carried out in the context of a leaky-box model: Particles are introduced into the SBN in the following ways: 1) they are injected; 2) all typical timescales are independent of position; 3) they escape if the interaction probability is smaller than the escape probability; otherwise, they lose energy and all of their byproducts are taken into account in the following step of propagation; and 4) they are no longer propagated if their energy is outside of our range of interest.
SimProp simulates UHECR propagation in the extragalactic medium under the assumption of a predetermined spectrum of injected particles. It should be noted that the propagation in the source depends on the source parameters but not on the spectral parameters. Due to this, a single flux (∝ E −1 ) is used for each set of source parameters when performing the in-source propagation. Re-weighting the elements spectra will provide the corresponding ejected spectra when spectral parameters (spectral index and rigidity cut) are modified. The benefit of this method is that it drastically cuts down on the amount of time needed to compute to explore the parameter space. For the final stage, the propagation from the SBN to the Earth, we use the same steps outlined in appendix A of [26].

Comparison to experimental data
The purpose of this study is to investigate the claim that SBGs are UHECR's sources. This is done by comparing experimental data to the CR flux as it has been affected by interactions in the SBNi and extra-galactic propagation. We adopt a measurement of the energy spectrum in log 10 (E/eV) bins of 0.1 width from 17.8 eV to 20.2 eV, obtained with the data collected over 15 years with the surface detector array of the Pierre Auger Observatory [28]. As for the depth of the shower maxima X max distributions [29], we consider log 10 (E/eV) bins of 0.1 from 17.8 eV to 19.6 eV, plus one additional larger bin containing events with energies above 10 19.6 eV; each X max distribution is binned in intervals of 20 g cm −2 . Following [26], we use the deviance D = −2ln(L/L sat ) as estimator of the agreement of our parametric model to data, where L is our model and L sat is a model that perfectly describes the data. The total deviance consists of two terms, D J and D X max . The former, which deals with the energy spectrum, is a product of Gaussian distributions, whereas the latter, which is used to fit the X max , is a result of multinomial distributions. Gumbel distribution functions with hadronic interaction model-dependent parameters are used to model them. For the current analysis we adopt EPOS-LHC [30] as the hadronic interaction model. It has been demonstrated in [27] that utilizing Sibyll 2.3d as the hadronic interaction model in the atmosphere has a minimal influence on the deviance value but has no significant negative effects on the main features of the propagated spectrum and composition, leading to the same general scenario.

Characterization of the parameter space
In what follows, we provide the set of free parameters and assumptions used in the current study to characterize the source environment, the injection parameters of the accelerated CRs, and the details of extra-galactic propagation, as well as the additional parameters required for the lowenergy region of the measured spectrum and composition.

Source parameters:
The total infrared luminosity of the SBG, L IR ,and the radius of the SBN region, R, are the free parameters related to the source. The former is allowed to range from the value 10 44 erg s −1 up to 10 46 erg s −1 , featured by the most powerful ultra-luminous infrared galaxies (ULIRGs), such as Arp 220 [22]. The latter is assumed to vary from a minimum value of 150 pc, as inferred for NGC253 [see e.g. 21], up to 250 pc as a standard value for the scale-height of thin disks in spiral galaxies like the Milky Way. The size and luminosity both have a significant impact on the transport of UHECRs. Energy losses are specifically impacted by L IR , whereas R has an impact on both interaction and escape time. The Kennicutt-Schmidt scaling [31] relates the target density n ISM to the star formation rate and, in turn, to the IR luminosity. Finally, it is assumed that the magnetic field B in the SBN is a fixed parameter with a representative value of 200 μG for the SBNi. The magnetic field's coherence length is set to 1 pc [see e.g. 21].

Injection parameters:
We assume that CRs are injected as a power-law spectrum of index γ, with maximum rigidity R cut , and that the injected flux is proportional to E −γ . The range of γ is specifically taken into account, with a minimum value of γ = 1 and a maximum value of γ = 2 as anticipated by the diffusive shock acceleration. Hard spectra from a variety of conceivable circumstances have previously been proposed in the literature. According to recent findings on the combined fit performed by the Pierre Auger Collaboration [26], we assume the range 10 18 − 10 19 eV for the rigidity threshold. For the sake of simplicity, we operate under the presumption that a single heavy nuclear mass, A, was injected into the SBN environment. The exploration of heavy nuclei fragmentation and the resulting formation of lighter byproducts is made possible by such an assumption. In agreement with [5] and [26] we fix A to the atomic mass value of Silicon-28.

Extra-galactic propagation:
Two photo-disintegration cross section models are implemented by SimProp: TALYS [32] and PSB [33] as well as two potential EBL models: Gilmore [34] and Dominguez [35]. In this study, we used the TALYS and Gilmore photo-disintegration cross sections and the EBL model, respectively, to compute the interactions in the source environment and in extra-galactic space. Finally, we make the assumption that the redshift distribution of our UHECR sources follows the evolution of the star formation rate (SFR) up to z max = 6.

Low-energy component:
We also add a significant CR flux below the ankle, similar to [5]. This is necessary because the measured mass composition does not match what is expected by the disintegration of nuclei in the source environment, which would result in only light nuclear fragments in the low energy range. Such a spectral component might belong to a distinct category of extra-galactic sources. [see e.g. 17] or the re-acceleration of Galactic CRs at the Galactic wind termination shock are examples of such phenomena. In our work, we take this new component to have a fixed spectral index γ = 4.2, with the nitrogen mass group as its dominant mass group [27], while we allow for a free normalization.

Results
We run a parameter space scan to determine the best configuration, starting from the set of parameters of our prototype SBG (see Tab. 1 of [23]), which is indicative of the most prevalent class of mild starbursts (including nearby sources like M82 and NGC253). The shape of the starformation-rate function, which shows that M82-like galaxies are the most prevalent and plentiful in the local Universe, led us to choose mild starbursts as the central focus of our investigation. The best source parameters found correspond to a luminosity ∼ 5 times higher than our starting prototype; such a luminosity is typical of a more active class of galaxies known as luminous infrared galaxies (LIRGs). Given the properties of the luminosity functions [36], SBGs with these characteristics are somehow less frequent than the prototype model. On the other hand, ULIRGs like Arp220 can be discovered at the highest end of the luminosity function, which is not where LIRGs may be found. As a result, it is possible to infer that galaxies with infrared luminosities above a particular threshold are more likely to have UHECR acceleration sites.
The spectrum and mass composition of UHECRs at Earth are shown in relation to the best fit parameters in Fig. 5. While our calculation produces a complex evolution of the mass composition with energy that closely resembles the ankle feature, the hierarchical order of the par-   Figure 5. Top: All-particle best-fit scenario and partial spectra related to different detected mass groups compared to the measured spectrum [37]. Bottom: Average (left panel) and standard deviation (right panel) of the experimental (black dots as reported in [38]) and expected (lines) X max distribution. tial fluxes follows the first two moments of the X max distributions. This is despite the simplicity of our model, which allows us to qualitatively reproduce the ankle feature with a precision of the order of ∼ 10% (at the energy of the ankle). Particularly, at low energies, light secondary nucleons predominate while at high energies, the composition becomes heavier as a result of the dominant escape rate over the interaction rate. Only below ∼ 10 18 eV the additional low-energy component is relevant; at higher energies, the contribution from SBGs dominates. We emphasize that, rather than using a more appropriate luminosity function, our results were obtained under the hypothesis of a single injected nuclear species and a single prototypical SBG representative for the entire population. The effects of such assumptions are discussed in the section 5.
The number of SBGs required to describe the data should also not be greater than the number of such galaxies detected in the nearby universe. It is possible to estimate the necessary emissivity ε needed to power the UHECRs at redshift z = 0 by comparing the model prediction with data as follows: where J inj is the spectrum injected at the source, before considering the interactions, and E min is an arbitrarily low energy value (here E min = 10 17 eV). As already discussed in [23], we define α as the ratio between the CR luminosity (L CR ) and IR luminosity obtained from the best fit (L IR = 1.2 · 10 45 erg/s). The number density of sources, n SBG , can be estimated as where α is normalized to 0.1 under the assumption of a sub-equipartition of the UHECRs compared to the background photon fields.
In [36] the luminosity density is computed as a function of the redshift using the luminosity functions. By integrating these functions for luminosity above the best fit one , we find: It can be seen that, assuming a fiducial value of α = 10%, the number density of sources derived from the integral of the luminosity function (eq. 5) is an order of magnitude greater than the one produced using our model (eq. 4). This discovery is encouraging since it can be explained by the fact that UHECRs in the starburst environment requires a lower energy budget.

Constraining the neutrino flux
We track the hadronic and photo-hadronic byproducts such gamma rays and neutrinos in addition to the escaping flux of UHECRs. Due to the strong magnetic fields characteristic of the SBN environment, which is highly opaque to gamma rays with energy 10 TeV, HE electron-positron pairs are predicted to cool via synchrotron rather than starting cascades that could result in spectral features in the TeV range [21]. There is hence no reason to expect a meaningful gammaray counterpart to the presence of UHECRs in the SBN environment. Due to this, we postpone to future works the exploration of the multi-wavelength effects of UHECRs in the SBN environment. Different from gamma rays, neutrinos travel almost undisturbed once they are produced. Both within the SBG environment and during the propagation of UHECRs from their sources to the Earth, we calculate the production of HE neutrinos. We simply take into consideration the adiabatic energy loss effect of the Universe's expansion on the neutrino flux due to their low interaction cross section. After oscillating through cosmological distances, we finally assume an average flavour composition for Earth (1:1:1).
In Fig. 6 we show the diffuse neutrino fluxes associated to the different contributions considered in this work: cosmogenic neutrinos (grey lines), namely those neutrinos produced by the interaction of UHECRs with the cosmic microwave background (CMB) and the extra-galactic background light (EBL), the neutrinos produced by photohadronic interactions of UHECRs in the source (blue lines) as well as the neutrinos produced by hadronic interactions of UHECRs in the source (magenta lines). The neutrino flux resulting from our calculation is compared with the Auger limit, the limit expected by GRAND [40] after three years of operation and the flux observed by Ice-Cube [41,42]. Two different model predictions, solid and dashed lines, are shown and compared in Fig. 6. They Prototype Best fit Figure 6. Single-flavour neutrino fluxes for the prototype case (dashed) and best case (solid), compared to the IceCube neutrino flux (both High Energy Starting Events [41] (points) and single power-law flux from muon-neutrinos [42] (band)), the expected sensitivity of GRAND after three years of operation [40] and to the limits for the cosmogenic neutrinos by the Pierre Auger Collaboration [39]. The shaded area is drawn in order to highlight the difference.
refer respectively to our best fit SBG and the prototype SBG (M82-like assumption). It's interesting to note that in the cases presented, the source-neutrinos produced in pγ interactions are almost equivalent to the expected cosmogenic fluxes, whereas the source-neutrinos produced in pp and pA interactions predominate the neutrino flux, which is of the same order of magnitude as the observed one. This finding suggests that the neutrino flux measured by IceCube above ∼ 100 TeV may be a direct outcome of UHECR accelerators found in SBGs. This multimessenger relationship suggests that we may be able to look at the galaxies that are emitting neutrinos with energies between 1 and 100 PeV in order to learn more about the origins of UHECRs.

Discussion and conclusions
In this work, we develop a source-propagation model to quantify the interactions occurring in the SBG environment and in the propagation to the Earth in order to investigate whether SBGs can be the sources of UHECRs. We precisely compute proton-proton and proton-nucleus interactions for the first time in the SBG environment and analyze the effects on the UHECR and neutrino fluxes in addition to photo-hadronic interactions. We compare our model prediction with the energy spectrum and mass composition observed by the Pierre Auger Observatory, working under the assumption that the sources of UHE-CRs are all characterized by some representative properties. We demonstrate how SBGs can qualitatively properly represent the composition and spectrum of the measured UHECR. Additionally, we compute the neutrino flux related to the transport of UHECRs and demonstrate how this enhances the constraining power of our model. This could enable us to consider whether a set of source parameters can adequately describe the UHECR data with-out exceeding the measured neutrino fluxes. In particular, we show that the UHECRs accelerators would contribute significantly to the neutrino flux measured by IceCube at energy 10 2 TeV if SBGs were hosting them. We observe that the expected neutrino flux from sources dominates the cosmogenic one. Additionally, we demonstrate that hadronic interactions might be essential for explaining the observed neutrino flux and deserve further research in source environments , since they could be able to hide the outcome in neutrinos from photo-hadronic interactions in the source environment and in the extra-galactic space. Therefore, in the near future, a complementary search for the origins of UHECRs could be guided by the identification of steady PeV neutrinos emitters.
We find that the data may be precisely described within the investigated parameter space if the UHECR nuclei are injected with a hard spectrum. Here, we emphasize that even if the conventional injection ∼ E −2 , typical of the diffusive shock acceleration, is not preferred, it is not necessarily impossible for such a mechanism to accelerate UHECRs. In fact, it is possible to obtain hard spectra in the context of diffusive shock acceleration under a number of core SBG conditions, such as multiple shocks, converging flows, particle re-acceleration, or transport conditions in the acceleration region that are different from those in the rest of the SBG.
The results that are being given here were obtained under the assumption that just one heavy nuclear species had been injected. The analysis of scenarios in which many masses are injected, each with a distinct relative composition that might also soften the spectral index discovered at the injection, goes beyond the scope of this work and is left for further research. Nonetheless, when the injected masses are heavier than silicon nuclei, the expected scenario fails. We specifically tested the iron nuclei at injection scenario and found that the data description significantly worsens. This is due to two factors: 1) The expected composition at Earth is too heavy compared to the observational data from the Pierre Auger Observatory; 2) The maximum rigidity of the iron nuclei is determined by comparing the expected and measured UHECR spectra; and 3) The maximum energy of the nucleons from the disintegration (1/A of the maximum energy of nuclei) is too low to adequately describe the data at the energy of the ankle.
Different luminosities of the prototype SBG and various particle spectra are just two of the hypotheses that have been investigated in relation to the source characteristics. A clear limitation of this research is that it only considers sources with standard characteristics. On the other hand, this approach can highlight the possibility of some average properties characterizing a class of UHECR sources. Instead of relying on a single prototype, an intriguing upgrade of the current study might make use of a catalog of the luminosity functions of galaxies. In fact, it is reasonable to assume that SBGs with varying luminosities might contribute to the energy spectrum at varying levels and probably better represent the Auger data. The sources of UHECRs could be constrained by many messengers using the HE neutrino flux. The ability of UHECR data to constrain increases as the measured neutrino flow is ex- ceeded. The assumption that all similar sources are distributed throughout the universe has, however, simply been taken into consideration up to a redshift of z = 6 in this work. Instead of the UHECR fluxes, which are predicted to originate close to z = 1, this hypothesis has a greater impact on the predicted neutrino fluxes. Additionally, it is anticipated that using luminosity functions as opposed to a single prototype will reduce and broaden neutrino fluxes.
Future developments of this work will involve the production and propagation of photons within SBGs, as additional multi-wavelength restrictions may be discovered, particularly in hard X-rays, where electrons and positron pairs are predicted to emit through synchrotron on the strong magnetic fields typical of SBGs. The expected photon fluxes can be compared to the actual data as the neutrino fluxes, strengthening the constraining power of our model. Finally, further advancements in sourcepropagation analyses can be expected due to the Pierre Auger Observatory upgrade [43] and its increased precision in determining the mass composition at the highest energies.