<EMAIL_ADDRESS>
# Modeling Hadronic Interactions in Ultra-High-Energy Cosmic Rays within
Astrophysical Environments: A Parametric Approach

Antonio Condorelli  [  Dipartimento di Fisica “Ettore Pancini”, Università
degli studi di Napoli “Federico II”  INFN - Sezione di Napoli, Complesso Univ.
Monte S. Angelo, I-80126, Napoli, Italy
Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  Sergio
Petrera  Gran Sasso Science Institute, Via Francesco Crispi 7, 67100,
L’Aquila, Italy
INFN/Laboratori Nazionali del Gran Sasso, via G. Acitelli 22, 67100, Assergi
(AQ), Italy

###### Abstract

Interactions of ultra-high energy cosmic-rays (UHECRs) accelerated in
astrophysical environments have been shown to shape the energy production rate
of nuclei escaping from the confinement zone. To address the influence of
hadronic interactions, Hadronic Interaction Models (HIMs) come into play. In
this context, we present a parameterization capable of capturing the outcomes
of two distinct HIMs, namely EPOS-LHC and Sibyll2.3d, in terms of secondary
fluxes, including escaping nuclei, neutrinos, photons, and electrons. Our
parametrization is systematically evaluated against the source codes, both at
fixed energy and mass, as well as in a physical case scenario. The comparison
demonstrates that our parameterization aligns well with the source codes,
establishing its reliability as a viable alternative for analytical or fast
Monte Carlo approaches dedicated to the study of UHECR propagation within
source environments. This suggests the potential for utilizing our
parameterization as a practical substitute in studies focused on the intricate
dynamics of ultra-high energy cosmic rays.

Suggested keywords

††preprint: APS/123-QED

##  I Introduction

The nature of the ultra-high-energy cosmic-rays (UHECRs) still remains
elusive. Yet there are candidate astrophysical sources that appear to possess
the necessary requirements to produce such energetic particles [1, 2, 3, 4,
5]. In most cases they are extended sources where, apart the actual
accelerators, the surrounding environment is filled with magnetic and
radiation fields, as well as gaseous matter.

Cosmic-ray interactions in plausible source environments have proven to be a
critical input for understanding the energy spectrum and mass-composition
provided by Pierre Auger Observatory [6, 7] and Telescope Array data [8, 9] .
To account for the progressive increase in energy of the mean logarithmic mass
number $\langle\ln A\rangle$, it is assumed that each individual nuclear
component has the same magnetic rigidity at the sources. In addition, by
assuming sources with an emission spectrum as $E^{-\gamma}$, the small mixture
of elements detected in data can be replicated if the spectral index
$\gamma\leq 1$ [10, 11, 12]. Such hard spectral indices are surprising at
first glance, as typical expectations from first-order Fermi shock
acceleration predict $\gamma$ to be $\gtrsim 2$. It has been suggested that
hard values for $\gamma$ may represent the effect of in-source interactions in
changing the ejected spectra if the index of protons ejected from the sources
is softer [13, 14, 15, 16, 17, 18]. It was also shown that such a requirement
is consistent with the data in the energy range across the ankle feature [19,
20].

Focusing on the interactions occurring in the environment surrounding the
sources imposes studying the hadronic interactions in these environments.
While negligible for extra-galactic propagation, hadronic interactions could
play an important role in a source environment. In fact, in some extended
sources (e.g. in the nucleus of Starbust Galaxies) where the high density of
target ($\simeq 10^{2}\ cm^{-3}$) and high value of magnetic field ($B\simeq
200\ \mu\rm G$) could significantly enhance their chance to happen.

Hadronic interactions are a fundamental ingredient to study in-source
propagation, because they significantly shape the shower, both in terms of
disintegration of accelerated nuclei and in the secondary production of
neutrinos and electromagnetic particles. As shown in [17], they can contribute
much more than photo-hadronic interactions in the expected neutrino flux from
the source environment.

Cosmic ray propagation codes largely used in UHECR community, such as SimProp
[21] and CRPropa [22], will be soon updated with modules simulating hadronic
interactions [23, 24], to be suited to study in all aspects source
environments.

Proton-proton and proton-nucleus interactions are normally taken into account
by using Monte Carlo (MC) codes based upon different Hadronic Interaction
Models (HIMs). There are a few models available to simulate the hadronic
interactions used in cosmic-ray physics: in particular, EPOS-LHC [25],
QGSJetII-04 [26] and Sibyll 2.3d [27] have been recently re-tuned after LHC
data and are the most commonly used.

The diversity of particles generated by HIM codes strongly hinders
simplification into conventional analytical or parametric models. Yet, for the
study of the production from cosmic rays of stable secondaries such as
neutrinos, photons, and electrons, some simplistic frameworks have been
employed to establish connections between the energies of these particles and
their parent protons. One such approximation revolves around the relationship
between the energy of a proton and that of its resultant pions, resulting in
assuming a fixed value for it ($\delta$-approximation, see e.g. [28]). Further
simplification often assumes fixed energy fractions for the photons, electrons
and neutrinos stemming from the pion and muon decays. However, it is crucial
to recognize the inherent oversimplification in these estimates, which provide
only rough insights alongside the assumption of a primary spectrum.

Other works aim for a more realistic approach, taking into consideration the
distribution of produced pions derived from HIM simulation codes and
consequently also of their secondary products, i.e., neutrinos and photons
[29, 28, 30]. Nevertheless, these parameterizations are limited to proton
primaries and do not explore the highest-energies regime, focusing in the
energy range $[10\ \rm GeV,10\ \rm PeV]$.

In our study, we demonstrate that in the context of the propagation of UHECRs
within an astrophysical source environment, its prevailing conditions permit
the parameterization not only of secondary stable leptons and photons but also
of the other interaction products responsible for the nuclear cascade therein.

Providing a simple and efficient way to simulate hadronic interactions at the
highest energies is thus a timely topic. For this purpose we have developed a
Monte-Carlo simulation code based on parametric functions describing the key
features crucial for in-source propagation. In the following, starting from
HIM simulations, we parameterize the most important quantities, namely the
cross section, the secondary production and the nuclear fragmentation. This
allows us to take into account the effects of the hadronic interactions in a
broad energy range by using only a few parametric functions. This procedure is
currently implemented for EPOS-LHC and Sibyll2.3d. In the future, further
models will be possibly considered.

The paper is organized as follows: we introduce the HIMs we aim to study and
detail the most important quantities for our study in Section II; the
parameterization of the most important quantities in an astrophysical scenario
is described in Section III; we compare the two models to the corresponding
source codes, in order to validate the parametric models and estimate their
accuracies and performances, in Section IV; we discuss the outcomes of our
simulation in a sample in-source case, showing the impact of our assumptions,
in Section V . We finally draw our conclusions in Section VI.

##  II HIMs

The advent of the LHC has brought about remarkable strides in particle physics
and our comprehension of high-energy interactions. One of the key areas of
focus has been the enhancement of Monte-Carlo generators, which are crucial
for simulating particle interactions at extreme energies. The primary emphasis
in this development has been on hadron production. For HIMs used in cosmic-ray
physics, this has led to substantial advancements in the modeling of hadrons,
their interactions and decays, with the utimate goal of understanding the
behavior of forward hadrons within the extensive cascades of secondary
particles. Even minor adjustments in the baryon and resonance production or
the introduction of collective hadronization can have far-reaching
consequences on observables, such as the generation of muons [31].

The propagation of UHECRs in a medium is determined by the following
parameters: the resident time $\tau_{\mathrm{res}}$, i.e. the duration UHECRs
spend within the medium; the interaction time $\tau_{\mathrm{int}}$, i.e. the
time interval before an interaction occurs; the decay time
$\tau_{\mathrm{dec}}$ which is relevant for unstable particles. These time
scales apply to both the primary cosmic rays and any potential secondary
particles generated.

The resident time is influenced by factors such as the size of the medium and
of the presence of magnetic fields. The interaction time depends upon the
cross-sections of the involved processes and the matter density. Decay time
relies on the particle’s mean lifetime and Lorentz factor.

In the atmospheric propagation of UHECRs, the crucial parameters are the
hadronic interaction times and decay times. The decay time, spanning a wide
range of possible values, holds particular significance, although the
interaction time also varies substantially among different particles.
Conversely, the resident time is predominantly dictated by the air column
density, typically measured in hundreds of g cm-2. As a result, in the
evolution of atmospheric cascades, forward hadrons play the primary role, with
their interaction and decay times constantly competing, while only very
forward nuclear fragments and nucleons actively participate.

The propagation regime within an astrophysical source displays significant
variability as far as hadronic interactions are concerned. Depending on the
source characteristics, it can range from a mere transport (or transparent)
regime, where the resident time is considerably shorter than the interaction
time ($\tau_{\mathrm{res}}\ll\tau_{\mathrm{int}}$), to a full cascade regime,
where the resident time greatly exceeds the interaction time
($\tau_{\mathrm{res}}\gg\tau_{\mathrm{int}}$). Even in scenarios where the
cascade regime dominates, the crucial interaction timescales primarily concern
nuclei and nucleons. In fact unstable hadrons in all sources have decay times
substantially shorter than their interaction times
($\tau_{\mathrm{dec}}\ll\tau_{\mathrm{int}}$), thereby impeding their
interaction with matter.

From all that it becomes evident that the role of HIMs in the propagation of
UHECRs within astrophysical sources differs significantly from its role in
extensive air showers. In source environments, a distinct division emerges
between the contributions of nuclei and nucleons, primarily driving the
cascade process, and that of unstable hadrons, which instead generate fluxes
of neutrinos, photons, and electrons.

In order to model the dynamics of hadronic interactions and their effects in
astrophysical sources, three features have been chosen to characterize
hadronic interactions: the inelastic cross sections, the rapidity
distributions of secondary particles, and the fragmentation of nuclei.

  * -

p-p and nucleus-p cross sections: these cross sections are crucial ingredients
to compute the interaction time. The timescale for the hadronic interaction
reads:

$\tau_{\mathrm{int}}=\large({n_{\rm
ISM}\,\sigma_{\mathrm{inel}}\,c\large)}^{-1},$ (1)

where $n_{\rm ISM}$ is the Interstellar medium (ISM) gas density and
$\sigma_{\mathrm{inel}}$ is the inelastic cross section for proton-proton or
nucleus-proton interactions. In this work we assume that matter is composed
only of protons or hydrogen molecules.

  * -

Secondary production: The multiplicity of secondary particles and their
rapidity spectra in both proton-proton and nucleus-proton interactions refer
to the number and variety of particles produced as a result of these
interactions. They affect primarily the observational signatures of
astrophysical sources, such as gamma-rays and neutrinos, but also the cascade
evolution in the source medium. A good description of secondary particles is
crucial for interpreting these observations.

  * -

Fragmentation of nuclei: After a nucleus in the ISM undergoes inelastic
process, it may fragment into various nuclei and particles. The study of these
fragmentation processes is crucial in astrophysics, because they affect the
development of the cascade and, in a second instance, the spectra of secondary
particles.

For this reason, in the following we will provide a description of the
approach used for the parameterization of these three features for the two
different HIMs.

##  III Parameterizations

The software that emulates hadronic interactions described in this paper is
made publicly available, as PARISH (PARametric simulation of In-Source
Hadronic interactions) [32]. This availability ensures transparent access to
the intricate details concerning the parametric assumptions incorporated
within the code. As a consequence, rather than providing in this paper an
exhaustive account of the parametric functions, we focus on demonstrating the
validation of key distributions. These distributions are compared with those
generated by the original HIM simulation codes to illustrate the level of
accuracy achieved for event features crucial for understanding the propagation
in UHECR sources. This validation is presented in the next two sections.

Figure 1: Inelastic cross-sections for nucleus-proton interactions for EPOS-
LHC (solid, transparent colored lines) and Sibyll2.3d (dashed black lines) for
five nuclei (1H, 4He, 14N, 28Si and 56Fe) in the energy-per-nucleon range
$10^{16}\div 10^{20}$ eV/nucleon.

###  III.1 Inelastic cross sections

Nucleus-proton inelastic cross sections have been derived from EPOS-LHC and
Sibyll2.3d simulations. They are easily fitted assuming logarithmic energy and
mass dependence as shown in Fig. 1. The fitted functions, though having
slightly different parameters, appear almost completely superimposed, in the
energy-per-nucleon range $10^{16}\div 10^{20}$ eV/nucleon.

###  III.2 Secondary particles

HIM simulations produce a wide variety of hadrons with number multiplicities
of the order of hundreds up to thousands, for energies in the UHE regime. As
discussed in the previous section, the typical matter density in astrophysical
sources is low enough to impede re-interactions of these particles, so all of
them proceed through consecutive decays until only stable particles survive
(nucleons, photons, electrons and neutrinos). Nucleons111 Neutrons decay with
a mean path of 9.2 $E$/EeV kpc. Whether they escape or decay depends on energy
and source size. Generally, the chance of a neutron interacting as itself,
rather than transforming into a proton post-decay, is very low.  only can
interact again in the source environment under certain conditions of matter
density and diffusion, otherwise escape. Photons and electrons can induce
further cascades in the radiation field depending on its strength, while
neutrinos always escape.

In these circumstances, pions and nucleons can stand as representative
entities for all the hadrons generated in the interactions. This choice is
supported by the following factors:

  * a)

Pions account for more than 70% of the total hadrons produced, while baryon-
antibaryon pairs are produced at approximately a 10% level.

  * b)

Pions actively contribute to the production of electrons and neutrinos through
charged pion-muon decays, as well as photons through neutral pion decays.

  * c)

All baryons (antibaryons) heavier than nucleons (antinucleons) undergo
successive decays until only nucleons (antinucleons) remain.

Therefore we adopt this “lightest hadrons” approximation, i.e. assuming that
the secondaries produced are only pions and nucleons. The two components serve
distinct functions: the former contributes to the spectra of neutrinos,
photons, and electrons, while the latter fuels the cascade process within the
source medium. We derive the parametric distributions for both components
(refer to the details below), utilizing simulated events from EPOS-LHC and
Sibyll2.3d222 For Sibyll2.3d, the original source code was used to generate
events. For EPOS-LHC, we used CRMC (Cosmic Ray Monte Carlo) [33], an interface
giving access to different generators. .

Our simulated data sample consists of nucleus-proton interactions in five
logarithmic energy bins, log10 ($E$/($A$ eV)) = {16, 17, 18, 19, 20} and four
mass numbers $A$ = {1, 4, 14, 56}. In each energy-mass bin, we have generated
5000 interactions for each HIM. These interactions include all hadrons that
are produced instantly post-hadronization. Additionally, another sample of
5000 interactions per energy-mass bin has been generated, including only
stable particles (post-decay), with the exception of neutrons, which are
considered pre-decay. The former sample has been used to parameterize
distributions of pions and nucleons, the latter to validate the ones of
neutrinos, electrons and photons. For the parameterization of nuclear
fragmentation also samples of the same size have been produced for 28Si and
40Ca nuclei, to have a better description of the high-mass region.

In our parametric model, we focus on few key observables, namely the
multiplicity distribution of all secondaries and the rapidity distributions of
pions and nucleons. The distributions of neutrinos, electrons, and photons are
not parameterized; rather, they are derived through MC simulation of the decay
processes involving pions and muons.

The approach to parameterization involves identifying a simple parametric form
capable of reproducing the distribution from the source code across all
energy-mass bins. Subsequently, we proceed determining the most
straightforward evolution of the best-fit parameters. In the majority of
cases, employing linear evolutions with respect to both the logarithms of
energy-per-nucleon and mass number provides a satisfactory description across
the entire energy-mass range. When dealing with pion rapidity distributions in
the context of EPOS-LHC, identifying a straightforward evolution of the
parameters proved challenging. In this specific scenario, we resorted to a
bilinear interpolation method to achieve an adequate representation.

Nucleon energy spectra are extracted from baryon and antibaryon rapidity
distributions using HIM source codes, operating under the assumption that
these distributions are uniform across all baryonic particles. Unlike pions,
their energy spectra do not demand an accurate modeling, as the cascade is
solely influenced by the nucleon interaction length, which in turn exhibits a
modest energy dependence. Consequently, we opt for a unified nucleon rapidity
distribution with identical parameters depending only on energy for both HIMs,
overlooking the nuances between the two models.

###  III.3 Nuclear fragmentation

Nuclear fragments, along with nucleons from secondaries (Sec. III.2),
contribute to the cascade within the source environment. Notably, each nucleus
interacts with ISM targets with comparable timescales (see eq.(1)), because
inelastic cross sections vary by at most a factor of 10 at a given energy-per-
nucleon.

Nuclear fragmentation is a slow process with respect to hadronic times,
predominantly originated by the evaporation of a nucleus excited by the
hadronic interaction (residual nucleus) [34, 35]. Upon examining the nuclei
ejected in events from HIMs under consideration, we observed that the majority
of fragment distributions exhibit independence of the primary energy, aligning
with expectations from the evaporation model. Consequently, we confidently
assumed the complete preservation of this energy-independence, and accordingly
parameterized the fragmentation distributions solely as a function of the mass
number, aggregating all energies present in the data sample.

Cosmic-ray collisions are predominantly peripheral, affecting only a small
number of nucleons and imparting relatively modest energy to the residual
nucleus (see e.g. [36]). In addressing the subsequent de-excitation process,
each HIM employs some statistical approach, involving the emission of
nucleons, deuterons or heavier fragments. In the context of UHECR interactions
and in the source reference system, these fragments maintain the Lorentz
factor of the initial projectile nucleus and follow its original direction.

The comparison between the outputs of the two HIMs reveals differences,
hindering a uniform treatment of the parameterization process. In Sibyll2.3d,
the nuclear fragments are distinctly categorized from the hadrons generated in
the interaction. The fragments included in this class encompasses nucleons and
heavier nuclei. Consequently, the mass of the residual nucleus responsible for
generating these fragments can be derived from the code output, along with the
count of wounded nucleons 333 The count of wounded nucleons is directly linked
to the mass of the residual nucleus. The mass number of the residual nucleus
may significantly differ from that of the projectile nucleus. . In contrast,
the output of EPOS-LHC exclusively identifies nuclear fragments with $A>1$ via
their $(A,Z)$ numbers, while nucleons originating from the evaporation remain
untagged. Our parameterization assumes that, with $A$ representing the mass
number of the projectile nucleus and $A_{\mathrm{frag}}$ the mass number of
the nuclear fragment, the evaporation process involves also
$A-A_{\mathrm{frag}}-1$ nucleons.

Regarding the distribution of fragment masses, they have been parameterized
for both HIMs. EPOS-LHC exhibits significant differences compared to
Sibyll2.3d (see Sec. IV). The former HIM generates a single heavy fragment
($A>1$), with rare exceptions when two 4He nuclei emerge from 8Be. All
together, the produced fragments are too light with too many free nucleons, as
pointed out in [37]. Sibyll2.3d displays a broad mass distribution of
fragments, necessitating a more intricate parameterization. This complexity
arises because the determination of mass numbers for the fragments must
consider the constraint imposed by the total mass number of the residual
nucleus.

Figure 2: Pion rapidity distributions for protons and 56Fe nuclei of 1017 and
1019 eV/nucleon colliding with target protons. Left (right) panels correspond
to EPOS-LHC (Sibyll2.3d). Black (red) line are generated with the
corresponding HIM source code (parametric code) simulation. Figure 3:
Distributions of total multiplicity for HIM source code (black) and pion
multiplicity for parametric model (red). Protons and 56Fe nuclei of 1017 and
1019 eV/nucleon colliding with target protons are shown.

##  IV Comparisons

In this section, we aim to demonstrate the validity of the parameterizations
outlined earlier by conducting thorough comparisons between the distributions
derived directly from the HIM source Monte Carlo codes and those generated by
the corresponding parametric model. Additionally, we provide a comparative
analysis with another established parameterization, focused on the production
of secondary photons and neutrinos from protons, which is further discussed
towards the conclusion of this section.

###  IV.1 Secondary particles

Our parametric model is built upon the assumption of the ‘lightest hadrons’
approximation. Concentrating solely on longitudinal kinematics, as extensively
admitted in the context of UHECRs, the generation of secondary particles is
derived exclusively from modeling the rapidity distribution of pions and
nucleons. We have conducted a comparison between our distributions and those
derived from HIM source simulations across each projectile energy-mass bin
explored in this study. For brevity, we present the comparison results at two
energy-per-nucleon (1017 and 1019 eV/nucleon) for two projectiles: protons and
56Fe.

Figure 2 illustrates the pion rapidity distribution across four energy-mass
bins. The black histograms represent data generated using the HIM source
codes, while the red histograms depict results derived from the parametric
model. The left set of panels corresponds to EPOS-LHC, whereas the right
panels correspond to Sibyll2.3d. The shapes of pion rapidity distributions are
reproduced to a satisfactory extent in the entire rapidity range and in all
cases. Parametric model tends to overestimate the HIM source rapidity
distributions (by $10\div 20\%$ for EPOS-LHC and $40\div 50\%$ for
Sibyll2.3d), notably within the central region. This overestimation is
necessitated by the partial inability of pion-muon decays to comprehensively
represent the cascade decays of unaccounted particles.

Fig. 3 illustrates the distributions comparing the total multiplicity for HIM
source codes (black) against the pion multiplicity for the parametric model
(red). From these histograms we derive the ratio between the mean multiplicity
of pions and all secondaries. Specifically, for EPOS-LHC, this ratio falls
within the range of $0.74$ to $0.78$, while for Sibyll2.3d, it ranges between
$0.90$ and 1.

Figure 4: Neutrino log${}_{10}(E$/eV) distributions for protons and 56Fe
nuclei of 1017 and 1019 eV/nucleon colliding with target protons. Same color
codes as in Fig. 2. Figure 5: Log${}_{10}(E$/eV) distributions of photons and
electrons, for protons and 56Fe nuclei of 1017 and 1019 eV/nucleon colliding
with target protons. Same color codes as in Fig. 2.

For each interacting nucleus of a given mass number and energy, the actual
number $n_{\pi}$ of pions is randomly generated from the parametric
multiplicity distribution. Subsequently, $n_{\pi}$ rapidity values are
independently and randomly drawn from the corresponding rapidity distribution.
Each rapidity value is converted into energy using the expression
$E_{\pi}\simeq\mu_{\pi}/2\times\exp{(y)}$, where $\mu_{\pi}\simeq$ 0.38 GeV
represents the fixed assumed mean pion transverse mass in the parametric
model. Following the decay of pions, and the subsequent muon decays for
charged pions, photons, neutrinos, and electrons are finally generated and the
energy spectra of these secondary particles directly stem from the preceding
distributions of pions. In the following, we will combine the energy spectra
of electrons and photons, referring to them as ‘e.m. particles’. Despite
originating from distinct interactions, they both progress together within the
source medium through the electromagnetic cascade in the resident radiation
field. The cascade process is not considered in this analysis, as it depends
on source characteristics that fall outside the scope of the models.

Figures 4 and 5 display the logarithmic energy spectra of neutrinos and e.m.
particles produced in interactions of 1017 and 1019 eV/nucleon protons and
56Fe nuclei with target protons. Again the HIM source histograms are depicted
in black and the parametric ones in red; the left panels refer to EPOS-LHC and
the right ones to Sibyll2.3d. The reproduced neutrino spectra demonstrate a
reasonable level of fidelity. In the case of EPOS-LHC, the parametric spectrum
aligns within 20% to the source spectrum, except for energies below about 1
GeV. However, for Sibyll2.3d, the parametric spectrum tends to underestimate
the source spectrum by approximately 20% to 40%, particularly noticeable in
certain energy intervals. Similarly, the energy spectra of e.m. particles are
also reproduced to a reasonable extent, even if with a tendency for the
parametric model to underestimate the HIM source spectrum. This
underestimation amounts to about 40% to 50% for EPOS-LHC and 20% to 40% for
Sibyll2.3d. One can argue that most of these disparities, for both neutrinos
and e.m. particles, stem from hadrons whose decay modes and kinematics are
poorly represented by pion-muon decays.

Figure 6: Distribution of the ratio of all baryons to the number of
secondaries for protons and 56Fe at an energy of 1019 eV/nucleon. Up (down)
panel corresponds to EPOS-LHC (Sibyll2.3d). Black (red) line are generated
with the corresponding HIM source code (parametric code) simulation. Figure 7:
Rapidity distributions of baryonic particles, for 56Fe nuclei of 1019
eV/nucleon colliding with target protons. Left (right) panel corresponds to
EPOS-LHC (Sibyll2.3d). Same colors as for Fig. 6.

In relation to nucleons, their quantity is determined by modeling the ratio of
all baryons to the number of secondaries, and studying its variation with the
logarithmic energy-per-nucleon and mass number, as derived from the HIM source
codes. This ratio is reasonably replicated across all energy-mass bins. In
Figure 6, the distributions of the baryonic fraction is presented for EPOS-LHC
(above) and Sibyll2.3d (below), along with their respective parametric models,
specifically for protons and 56Fe at an energy of 1019 eV/nucleon.

Nucleon rapidity distributions have been studied for all the energy-mass bins.
As an example, Figure 7 shows the rapidity distributions for 56Fe nuclei of
1019 eV/nucleon colliding with target protons. Modeling encompasses both the
central region and the projectile leading peak. Instead, the target peak is
excluded from consideration, given that the associated energy range falls
significantly below the objectives outlined in this paper. It is important to
note that the rapidity distributions include all baryonic particles in the HIM
source code (depicted by the black histogram), while the parametric code
represents the distribution for only protons or neutrons (in red), in
accordance with the ‘lighest hadrons’ approximation. Generally, the unified
nucleon rapidity distribution assumed for both models better reproduces
Sibyll2.3d than EPOS-LHC, which exhibits more complexity. However, the
resulting differences in nucleon energy lead to minimal changes in the
interaction length, as observed earlier.

###  IV.2 Nuclear fragmentation

Figure 8: Distribution of the mass number of nuclei from the nuclear
fragmentation of 14N and 56Fe. Up (down) panel corresponds to EPOS-LHC
(Sibyll2.3d). Black (red) line are generated with the corresponding HIM source
code (parametric code) simulation; parametric histograms are shaded for a
better comparison with the source HIM ones. Figure 9: Energy spectra of the
nuclear fragments from 14N and 56Fe interactions at 1019 eV/nucleon. Up (down)
panels correspond to EPOS-LHC (Sibyll2.3d). Same colors as for Fig. 8.

In Fig. 8, the distribution of nuclear mass numbers is shown for all
fragments. The primary nuclei selected in this figure are 14N and 56Fe; the
upper panels correspond to EPOS-LHC, while the lower panels represent
Sibyll2.3d. The HIM source distributions are shown as black histograms, the
parametric distributions are drawn as red shaded histograms.

A notable observation is the significant disparity in the treatment of nuclear
fragmentation between the two HIM source codes (black histograms). In EPOS-
LHC, highly unstable nuclei within the lower mass range (e.g., 5 and 8) are
absent from the list of fragments. Additionally, an evident odd-even mass
number effect is observed. These distinct nuclear characteristics are
presumably handled within the source code itself. Conversely, Sibyll2.3d does
not consider nuclear binding energy and stability, leading to the generation
of all nuclei without discrimination, thereby leaving the treatment of
unstable nuclei to the user. In our parameterization, we retained these
individual features unchanged.

When comparing the outcomes of HIM simulations with those of the parametric
model, we note the reproduction of the most significant features observed in
the original HIMs. However, for Sibyll2.3d we notice some deficiency in mass
numbers within the high mass range, particularly near the mass number of the
projectile. This deficiency results in a tendency to favor lower mass numbers.
Regarding EPOS-LHC, while there is noticeable reproduction for 14N,
difficulties manifest when attempting to produce lower mass numbers, in the
case of 56Fe. All these disparities are expected to diminish their impact when
considering propagation in astrophysical environments, as light nuclei are
generated alongside heavier fragments to compose the final residual nucleus
mass. Given the comparable values of interaction mean paths among all nuclear
fragments, the main aspects of the cascading process are substantially
preserved.

The energy spectra of the nuclear fragments exhibit a reasonable level of
reproduction, as depicted in Fig. 9. This outcome is a direct result of
preserving the Lorentz factor during the interaction. Consequently, the shape
of the energy spectrum is determined by the distribution of fragment mass
numbers.

###  IV.3 Usage and performances

As outlined in the previous section, the parametric simulation code (PARISH)
presented here is specifically tailored for integration into simulation
frameworks within the context of extended astrophysical sources. The primary
function within these frameworks is to calculate the spectra of both nuclei
and secondary particles that escape from the source. These simulation
packages, such as CRPropa and SimProp, incorporate parametric models of photo-
nuclear interactions. Consequently, our Monte Carlo simulation of hadronic
interactions is ideally suited to complement and run concurrently with photo-
nuclear interactions within these frameworks.

Efforts have already been made to integrate the original HIM source codes
directly into propagation codes, either by interfacing with the HIM Fortran
code itself [17], by utilizing interaction tables generated through the CRMC
interface [38] or implementic a frontend interface to HIM’s [24]. However, the
adoption of a relatively simple C++ code, comprising around 1500 lines, as
opposed to the original extensive source codes (approximately 25,000 and
82,000 lines for Sibyll and EPOS, respectively), facilitates straightforward
insertion and management within the propagation codes.

Another crucial consideration is the CPU time required for generating hadronic
interactions. For all the timing tests we utilized, as a reference, the same
Linux machine with Intel(R) Core(TM) i7-7700 CPU @3.60GHz, 7.7GB of total
memory. The system is running Ubuntu 18.04.6 LTS and the installed C++
compiler is g++ (Ubuntu 7.5.0-3ubuntu1 18.04) 7.5.0. The PARISH execution time
per collision consistently remains below 6 ms (from 0.9 to 5 ms for
Sibyll2.3d, $\approx$ 6 ms in all cases for EPOS-LHC). Remarkably, this value
exhibits minimal variance across a broad energy range spanning from $10^{16}$
to $10^{20}$ eV/nucleon, as well as for various projectile nuclei ranging from
proton to iron. Moreover, of significant importance is the observation that
the execution time remains virtually identical for both the generation of
secondary hadrons and secondary stable particles.

The execution times for the same machine using the HIM source codes exhibit an
increase with energy and projectile mass. Specifically, with Sibyll2.3d (EPOS-
LHC), the duration per collision spans from 1 ms (18 ms) for protons at
$10^{16}$ eV/nucleon to 9 ms (360 ms) for iron nuclei at $10^{20}$ eV/nucleon.

However, the most extensive utilization of HIM source codes involves their
management of decays for all generated unstable hadrons. In this scenario,
execution times per collision increase, ranging from 15 ms (22 ms) for protons
at $10^{16}$ eV/nucleon to 41 ms (430 ms) for iron at $10^{20}$ eV/nucleon
with Sibyll2.3d (EPOS-LHC).

The PARISH simulation, which includes particle decays, proves to be
significantly more expedient than the source codes, with speeds ranging from
approximately 10 to 70 times faster.

Figure 10: Differential photon yield per collision of $10^{17}$ eV/nucleon
primaries. The black histograms represent photons generated using the
Sibyll2.3d source code for proton (lower) and iron (upper) collisions. The red
(blue) histogram corresponds to the Sibyll2.3d parameterization proposed in
this study for protons (iron) collisions. The dashed line illustrates the
analytical parametric function for photons provided in [29]. The vertical gray
line marks the lower limit of applicability for this function.

###  IV.4 Comparison with other models for secondary fluxes

In the field of High Energy cosmic ray physics, models concerning hadronic
interactions primarily focus on the production of secondary neutrinos,
electrons, and photons.

Simplistic approaches approximate the behavior of secondary leptons and
photons as delta functions at a certain fraction of the primary energy. These
offer only rudimentary insights alongside the assumption of a CR power-law
differential spectrum, since pion multiplicity and the kinematics of their
decays play an important role. The most widely accepted model is the one
introduced by Kelner et al. [29]. The model provides simple analytical
parametric functions for secondary photons and leptons, adaptable to any
primary proton distribution. Only proton primaries are considered. Secondary
particle spectra are characterized in terms of Feynman-$x$, the ratio of
energy transferred from the incident proton to the secondary particle,
relative to the primary proton energy. The model parameters are fitted using
inclusive secondary meson spectra derived from Sibyll distributions, within
the energy range of primary protons spanning from $10^{11}$ to $10^{17}$ eV,
with the resulting secondary spectra being applicable for $x>10^{-3}$.

Figure 10 illustrates a comparison of photon yields between the
parameterization proposed in this study and that formulated by Kelner et al.
for primaries at $10^{17}$ eV/nucleon. Additionally, the spectra from the
Sibyll2.3d source code are depicted for both proton and iron collisions within
the same figure. The photon distributions for protons from both
parameterizations align well with those generated by Sibyll2.3d and with each
other, particularly for photon energies exceeding $10^{14}$ eV (i.e.,
$x>10^{-3}$). Furthermore, our parameterization demonstrates remarkable
fidelity in reproducing photon spectra for iron nuclei as well.

##  V A physical case:
UHECR propagation in Ultra-luminous infrared galaxies

Figure 11: Escaping fluxes from the Starburst nucleus of Arp220 considering
different injected nuclear species: protons (top-left panel), helium (top-
right panel), nitrogen (bottom-left panel) and silicon (bottom-right panel).
Dashed lines refer to Sibyll2.3d, solid lines to our parametric model. The
escaping fluxes are normalized to the injection spectrum. The mass numbers of
escaping fluxes are grouped as follows: $A=1$ (red), $2\leq A\leq 4$ (grey),
$5\leq A\leq 22$ (green), $23\leq A\leq 28$ (cyan). The escaping fluxes are
normalized to the injection spectrum.

To assess the efficacy of the proposed parameterization, we give a concrete
illustration where both the hadronic interaction model and the
parameterization are employed in an astrophysical scenario. n this particular
example, Sibyll2.3d has been taken into consideration, utilizing SimProp as
propagation code.

When UHECRs propagate within the source environment, a cascade involving
nucleons and nuclear fragments is generated, whose evolution depends on
various source features such as size, radiation, and magnetic fields. Hadronic
interactions occur at progressively lower energies, shaping the escaping
nuclei, which form the observable energy spectrum and composition injected
into the extra-galactic space. Regarding the elemental energy spectra detected
at UHECR observatories, the chosen energy range for parameterization is
considered adequate for producing accurate model predictions.

In contrast, secondary stable particles resulting from the same interactions
possess energies at a per mil level compared to the colliding nucleus. Given
the scientific importance of investigating the spectra of secondary neutrinos
and electromagnetic particles, a key objective in multi-messenger
astrophysics, we have carefully addressed the low-energy extrapolation of our
model to maintain acceptable accuracy and robustness. For this reason, the
multiplicity of secondary particles and the rapidity distribution of these
particles are tuned to best replicate the expected behavior according to the
HIMs down to the 1012 eV energy range.

A perfect environment where UHECR acceleration and high gas density could
coexist is the core of ultra-luminous infrared galaxies (ULIRGs) [39]. They
are very dense environments, with a high rate of star formation and supernova
explosions. They are thought to be sites of cosmic-ray acceleration, and are
predicted to emit $\gamma$-rays in the GeV to TeV range. A high star formation
rate means lots of young, massive, stars, which radiate mostly in the UV and
have short lifespans on the order of tens of millions of years. The UV
emission is absorbed by the interstellar dust and re-emitted in the far-
infrared (FIR).

We focus in particular on Arp220 [40] , the closest ULIRG to Earth, at a
distance of about 75 kpc. It is the product of a merger of two galaxies and
retains two dense nuclei, sites of very high star formation. Due to its high
density and close proximity to Earth, it is considered a good candidate to
search for high energy contribution from its star forming regions. We choose
to model the environment following [41]: a leaky box model is used, where
interactions with photons and protons are taken into account, considering
parameters at the source listed in Table 1. The procedure adopted for this
simulation is similar to [17].

Parameter | Value
---|---
$R\ \rm(pc)$ | 250
$B\ \rm(\mu G)$ | 500
$n_{\rm ISM}\ (\rm cm^{-3})$ | 3500
${U_{\rm eV\ cm^{-3}}^{\rm FIR}\ \bigg{[}\dfrac{kT}{\rm meV}\bigg{]}}$ | 31312 [3.5]
${U^{\rm OPT}_{\rm eV\ cm^{-3}}\ \bigg{[}\dfrac{kT}{\rm meV}\bigg{]}}$ | 1566 [350]
Table 1: Parameters used to compute UHECRs propagation in the Arp220
environment (from [41]). The same notations of [17] are used.

In addition to considering interactions, it is essential to include diffusion
in the environment; in fact, charged particles may remain confined for an
extended period before eventually escaping. The diffusion timescale is given
by $\tau_{\rm D}=R^{2}/D$, where $R\equiv pc/q\approx E/q$ represents the
magnetic rigidity of a UHECR with energy $E$ and charge $q$. The diffusion
coefficient, computed within the framework of quasi-linear theory [42], is
expressed as $D\simeq cr_{L}^{2-\delta}\,l_{\rm c}^{\delta-1}/3$. Here,
$r_{L}=E/qB$ denotes the particle Larmor radius, $l_{\rm c}$ is the coherence
length of the magnetic field, $\delta$ is the slope of the turbulence power
spectrum, and $B$ represents the strength of the turbulent magnetic field. We
adopt $\delta=5/3$ as recommended for a Kolmogorov turbulence cascade.

We assume an accelerated spectra $\Phi(E)$ following the diffusive shock
acceleration with $\Phi(E)\propto E^{-2}$ and no cutoff.

We inject $10^{4}$ particles logarithmically distributed in the energy range
$10^{17}-10^{20}\ \rm eV$, where our parameterizations are valid. Note that
the purpose of this exercise is not to test Arp220 as a possible source of
UHECRs, but rather to estimate the effect of the surrounding environment on
the propagation of UHECRs and thus to show that our parametric model describes
escaping particles as well as HIM source code.

In Fig. 11 we show the transparency, i.e. the ratio between the escaping
fluxes from the Arp220 environment and the injected ones, for different
injected masses, from proton to silicon nuclei; dashed lines refer to
Sibyll2.3d, solid lines to our parametric model. One can notice how the
proposed parameterization can well describe the output of Sibyll2.3d in terms
of secondary nuclei, with significant differences only at very low
transparencies. These differences reflect in very small changes at the escape
of the environment, meaning that the development of the cascade in the medium
is accurately followed. One can also notice how the proton transparency
present a maximum at energy $\simeq 10^{19}\ \rm eV$, which corresponds at the
energy where the interplay between photo-hadronic and hadronic interaction
times takes place.

Figure 12: Neutrinos (left-hand side) and electromagnetic fluxes (right-hand
side) associated to UHE proton propagation in Arp220. Dashed lines refer to
Sibyll2.3d, solid lines to our parametric model. The escaping fluxes are
normalized to the injection spectrum.

To test the secondary production, we extend our injection range down to
$10^{12}\ \rm eV$, adopting the same number of primaries per energy bin as for
nuclei. The production of neutrinos, photons and electrons in source
environment are instead shown in Fig.12, where we compare the yield,
$y={\frac{dN_{i}}{dE}}/{\frac{dN_{\rm inj}}{dE}}$ i.e. the ratio of the number
of secondary particles produced between $E$ and $E+dE$ to the number of
protons injected in the same energy range.

The expected fluxes of neutrinos, photons and electrons produced in source
environment are influenced by a dual interplay: firstly, they are shaped by
the development of the nuclear cascade, and secondly, by the production of
secondary particles. Also in this case a good agreement is reached, especially
having in mind the uncertainties associated to the modeling of the source
environment. It is crucial to highlight the distinct meanings of the two
fluxes. In the case of neutrinos, the flux represents the actual neutrino
spectrum emerging from the source environment. Conversely, for electromagnetic
particles, determining the actual escaping flux necessitates incorporating the
electromagnetic cascade in the source radiation field, a process not addressed
in the current study.

##  VI Conclusions

In this paper, we introduce a parameterization capable of describing the most
important features of two different HIMs, Sibyll2.3d and EPOS-LHC, when
hadronic interactions occur in an astrophysical source environment. This
parameterization characterizes secondary fluxes, including disintegrated
nuclei, neutrinos, photons, and electrons. Validation procedures demonstrate a
good agreement between our parameterization and the underlying source code.
Additionally, we illustrate this agreement within a physical scenario, showing
how the combined effects of nucleus fragmentation and secondary production
replicate cascades across a broad energy spectrum.

Our software, facilitating efficient emulation of the tested HIMs, is now
publicly accessible, empowering the broader scientific community to enhance
their predictions regarding hadronic interactions within source environments,
without necessitating direct engagement with the HIMs themselves.

This study represents a significant advancement with respect to the previous
literature; Introducing a dedicated parameterization for hadronic interactions
is a novel approach for the UHECR community. We have shown that our proposed
parameterizations offer competitiveness compared to direct HIM usage, with
distinct advantages in terms of usability and computational efficiency.

Furthermore, there is potential for extending our parameterizations to lower
energies, providing a versatile framework applicable to galactic cosmic rays
as well.

The implications of this work extend to multi-messenger astrophysics, where
these interactions play a pivotal role in shaping observed UHECR spectra.
Accurately computing hadronic interactions within source environments is
essential for precisely modeling UHECR energy spectra, compositions, and
arrival directions. This, in turn, offers invaluable insights into the
physical processes within the sources responsible for their acceleration.

###### Acknowledgements.

AC gratefully acknowledges funding from ANR via the grant MultI-messenger
probe of Cosmic Ray Origins (MICRO), ANR-20-CE92-0052. The authors extend
their appreciation to Francesco Salamida for valuable feedback that
contributed to the development of the software. Special thanks to Denise
Boncioli for her unwavering encouragement and support throughout the duration
of this project. Lastly, thanks go to Felix Riehn and Tanguy Pierog for their
assistance in navigating the HIM codes.

## References

  * [1] Mauricio Bustamante, Philipp Baerwald, Kohta Murase, and Walter Winter.  Neutrino and cosmic-ray emission from multiple internal shocks in gamma-ray bursts.  Nature Communications, 6(1), April 2015.
  * [2] Luis Alfredo Anchordoqui.  Acceleration of ultrahigh-energy cosmic rays in starburst superwinds.  Phys. Rev. D, 97:063010, Mar 2018.
  * [3] Ke Fang and Kohta Murase.  Linking high-energy cosmic particles by black-hole jets embedded in large-scale structures.  Nature Physics, 14(4):396–398, jan 2018.
  * [4] T. Piran and P. Beniamini.  Ultra high energy cosmic rays from tidal disruption events.  Journal of Cosmology and Astroparticle Physics, 2023(11):049, nov 2023.
  * [5] Enrico Peretti and Markus Ahlers.  Particle acceleration and multi-messenger emission from ultra-fast outflows.  PoS, ICRC2023:361, 2023.
  * [6] P. Abreu, M. Aglietta, J. M. Albury, and al.  The energy spectrum of cosmic rays beyond the turn-down around $10^{17}$ eV as measured with the surface detector of the pierre auger observatory.  The European Physical Journal C, 81(11), nov 2021.
  * [7] Eric William Mayotte et al.  Measurement of the mass composition of ultra-high-energy cosmic rays at the Pierre Auger Observatory.  PoS, ICRC2023:365, 2023.
  * [8] R.U. Abbasi, Y. Abe, T. Abu-Zayyad, and al.  The energy spectrum of cosmic rays measured by the telescope array using 10 years of fluorescence detector data.  Astroparticle Physics, 151:102864, 2023.
  * [9] R. U. Abbasi, M. Abe, T. Abu-Zayyad, and al.  Mass composition of ultrahigh-energy cosmic rays with the telescope array surface detector data.  Phys. Rev. D, 99:022002, Jan 2019.
  * [10] R. Aloisio, V. Berezinsky, and P. Blasi.  Ultra high energy cosmic rays: implications of Auger data for source spectra and chemical composition.  JCAP, 10:020, 2014.
  * [11] Andrew M. Taylor, Markus Ahlers, and Dan Hooper.  Indications of Negative Evolution for the Sources of the Highest Energy Cosmic Rays.  Phys. Rev. D, 92(6):063011, 2015.
  * [12] Alexander Aab et al.  Combined fit of spectrum and composition data as measured by the Pierre Auger Observatory.  JCAP, 04:038, 2017.  [Erratum: JCAP 03, E02 (2018)].
  * [13] Noemie Globus, Denis Allard, and Etienne Parizot.  A complete model of the cosmic ray spectrum and composition across the Galactic to extragalactic transition.  Phys. Rev. D, 92(2):021302, 2015.
  * [14] Michael Unger, Glennys R. Farrar, and Luis A. Anchordoqui.  Origin of the ankle in the ultrahigh energy cosmic ray spectrum, and of the extragalactic protons below it.  Phys. Rev. D, 92(12):123001, 2015.
  * [15] A. D. Supanitsky, A. Cobos, and A. Etchegoyen.  Origin of the light cosmic ray component below the ankle.  Phys. Rev. D, 98(10):103016, 2018.
  * [16] Daniel Biehl, Denise Boncioli, Anatoli Fedynitch, and Walter Winter.  Cosmic-Ray and Neutrino Emission from Gamma-Ray Bursts with a Nuclear Cascade.  Astron. Astrophys., 611:A101, 2018.
  * [17] Antonio Condorelli, Denise Boncioli, Enrico Peretti, and Sergio Petrera.  Testing hadronic and photohadronic interactions as responsible for ultrahigh energy cosmic rays and neutrino fluxes from starburst galaxies.  Phys. Rev. D, 107(8):083009, 2023.
  * [18] Marco Stein Muzio and Glennys R. Farrar.  Constraints on the Hosts of UHECR Accelerators.  Astrophys. J. Lett., 942(2):L39, 2023.
  * [19] Quentin Luce, Sullivan Marafico, Jonathan Biteau, Antonio Condorelli, and Olivier Deligny.  Observational Constraints on Cosmic-Ray Escape from Ultrahigh-energy Accelerators.  Astrophys. J., 936(1):62, 2022.
  * [20] A. Abdul Halim et al.  Constraining the sources of ultra-high-energy cosmic rays across and above the ankle with the spectrum and composition data measured at the Pierre Auger Observatory.  JCAP, 05:024, 2023.
  * [21] R Aloisio, D Boncioli, A.F Grillo, S Petrera, and F Salamida.  Simprop: a simulation code for ultra high energy cosmic ray propagation.  Journal of Cosmology and Astroparticle Physics, 2012(10):007–007, Oct 2012.
  * [22] Eric Armengaud, Günter Sigl, Tristan Beau, and Francesco Miniati.  Crpropa: A numerical tool for the propagation of uhe cosmic rays, $\gamma$-rays and neutrinos.  Astroparticle Physics, 28(4-5):463–471, Dec 2007.
  * [23] Antonio Condorelli.  MultI-messenger probe of Cosmic Ray Origins: MICRO project.  EPJ Web Conf., 283:04011, 2023.
  * [24] Leonel Morejon and Karl-Heinz Kampert.  Implementing hadronic interactions in CRPropa to study bursting sources of UHECRs.  PoS, ICRC2023:285, 2023.
  * [25] T. Pierog, Iu. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner.  Epos lhc: Test of collective hadronization with data measured at the cern large hadron collider.  Physical Review C, 92(3), Sep 2015.
  * [26] S. Ostapchenko.  Monte carlo treatment of hadronic interactions in enhanced pomeron scheme: Qgsjet-ii model.  Physical Review D, 83(1), Jan 2011.
  * [27] Felix Riehn, Ralph Engel, Anatoli Fedynitch, Thomas K. Gaisser, and Todor Stanev.  Hadronic interaction model sibyll 2.3d and extensive air showers.  Physical Review D, 102(6), Sep 2020.
  * [28] Ervin Kafexhiu, Felix Aharonian, Andrew M. Taylor, and Gabriela S. Vila.  Parametrization of gamma-ray production cross sections for p p interactions in a broad proton energy range from the kinematic threshold to PeV energies.  Phys. Rev. D, 90(12):123014, December 2014.
  * [29] S. R. Kelner, Felex A. Aharonian, and V. V. Bugayov.  Energy spectra of gamma-rays, electrons and neutrinos produced at proton-proton interactions in the very high energy regime.  Phys. Rev. D, 74:034018, 2006.  [Erratum: Phys.Rev.D 79, 039901 (2009)].
  * [30] Esteban Roulet and Francesco Vissani.  On the energy of the protons producing the very high-energy astrophysical neutrinos.  Journal of Cosmology and Astroparticle Physics, 2021(03):050, mar 2021.
  * [31] Tanguy Pierog.  Open issues in hadronic interactions for air showers.  EPJ Web Conf., 145:18002, 2017.
  * [32] Antonio Condorelli and Sergio Petrera.  PARametric simulation of In-Source Hadronic interactions (v1.0), doi = https://doi.org/10.5281/zenodo.11057114, 2024.
  * [33] Ralf Ulrich, Tanguy Pierog, and Colin Baus.  The Cosmic Ray Monte Carlo Package, CRMC (v2.0.1), 2021.
  * [34] William A. Friedman.  Rapid massive cluster formation.  Phys. Rev. C, 42:667–673, Aug 1990.
  * [35] William A. Friedman and William G. Lynch.  Yields of medium mass nuclear fragments: Statistical emission.  Phys. Rev. C, 28:950–951, Aug 1983.
  * [36] J. Engel, T. K. Gaisser, T. Stanev, and Paolo Lipari.  Nucleus-nucleus collisions and interpretation of cosmic ray cascades.  Phys. Rev. D, 46:5013–5025, 1992.
  * [37] Tanguy Pierog, Matias Perlin, and Klaus Werner.  Hadronic Interactions Studies at the LHC.  EPJ Web Conf., 283:05001, 2023.
  * [38] Marco Stein Muzio, Glennys R. Farrar, and Michael Unger.  Probing the environments surrounding ultrahigh energy cosmic ray accelerators and their implications for astrophysical neutrinos.  Phys. Rev. D, 105(2):023022, 2022.
  * [39] E. Kilerci Eser, T. Goto, and Y. Doi.  Ultraluminous Infrared Galaxies in the AKARI All-sky Survey.  Astrophys. J. , 797(1):54, December 2014.
  * [40] Martina C. Wiedner, C. D. Wilson, A. Harrison, R. E. Hills, O. P. Lay, and J. E. Carlstrom.  Interferometric observations of the nuclear region of Arp 220 at submillimeter wavelengths.  Astrophys. J., 581:229–240, 2002.
  * [41] Enrico Peretti, Pasquale Blasi, Felix Aharonian, and Giovanni Morlino.  Cosmic ray transport and radiative processes in nuclei of starburst galaxies.  Mon. Not. Roy. Astron. Soc., 487(1):168–180, 2019.
  * [42] Jungpyo Lee, John Wright, Nicola Bertelli, Erwin F. Jaeger, Ernest Valeo, Robert Harvey, and Paul Bonoli.  Quasilinear diffusion coefficients in a finite larmor radius expansion for ion cyclotron heated plasmas.  Physics of Plasmas, 24(5):052502, may 2017.