SimProp: UHECR Propagation Simulator

• SimProp is a one-dimensional Monte Carlo code for UHECR propagation that tracks energy losses, photo-disintegration, and secondary production in extragalactic settings.
• It employs both continuous and stochastic methods to simulate adiabatic losses, pair production, and photopion as well as photonuclear interactions for rapid parameter scans.
• Validated against kinetic solutions and codes like CRPropa, SimProp serves as a valuable tool for exploring cosmogenic neutrinos, gamma rays, and composition studies.

SimProp is a publicly available one-dimensional Monte Carlo code for the extragalactic propagation of ultra-high-energy cosmic rays (UHECRs) in the absence of magnetic deflections. It was introduced to follow protons and nuclei from a source redshift down to $z=0$, accounting for cosmological energy losses, photo-disintegration, photo-pion production, and later the production of cosmogenic neutrinos, photons, and electron-positron pairs. Across its development, SimProp has remained defined by a deliberate trade-off: CRPropa is described as a detailed and extensive simulation code, whereas SimProp aims to achieve acceptable results using a simpler code, with particular emphasis on rapid scans of source spectra, composition, and cosmological evolution (Aloisio et al., 2012, Batista et al., 2015, Aloisio et al., 2017).

1. Origins, scope, and conceptual position

Aloisio et al. introduced SimProp as a mono-dimensional Monte Carlo framework for UHECR propagation in cosmological backgrounds, explicitly neglecting magnetic deflections and following particles only in redshift or cosmic time. In its early formulation, the code was inspired by the kinetic approach of Aloisio et al., and used a reduced nuclear description: either a small set of representative nuclei, chosen so that their photo-disintegration thresholds and cross-section shapes span the range of interest, or a single stable isotope per mass number $A$, depending on version and application (Aloisio et al., 2013, Aloisio et al., 2012).

This design placed SimProp in a specific methodological niche. It was intended to predict energy spectra and mass composition at Earth for arbitrary source distributions and injection spectra with enough speed for large parameter scans, rather than to provide a full three-dimensional trajectory integration or a complete isotope-resolved nuclear network. The main approximations stated repeatedly in the literature are rectilinear propagation, simplified nuclear channel structure, and reduced photohadronic kinematics relative to more detailed transport packages (Batista et al., 2015).

A frequent misconception is to treat SimProp as a generic UHECR transport framework without qualification. In the published descriptions, however, its defining scope is narrower: it is a one-dimensional extragalactic propagation code, and its original public releases were tailored to primaries between $^{1}\mathrm{H}$ and $^{56}\mathrm{Fe}$, with later work extending this range only through explicit code modifications (Matteo, 4 Sep 2025).

2. Propagation formalism and Monte Carlo workflow

The core workflow is event-by-event. A primary cosmic ray, either a proton or a nucleus, is injected at a source redshift drawn from a user-defined distribution, commonly uniform in comoving volume up to some $z_{\max}$, with an energy spectrum

$$\frac{dN}{dE}\propto E^{-\gamma}\exp\!\left[-\frac{E}{Z R_{\rm cut}}\right].$$

The particle is then propagated from its source to $z=0$ in small redshift intervals $\Delta z$, or equivalently in small $\Delta t$ steps in later implementations (Batista et al., 2015).

At each step, SimProp applies continuous energy losses from Hubble expansion and, where relevant, $e^+e^-$ pair production. The adiabatic term is written as

$A$0

or equivalently $A$1, depending on notation. For pair production, the code uses a continuous energy-loss approximation, with the nuclear rate scaled from the proton rate as $A$2 in the formulations that keep the explicit charge dependence (Batista et al., 2015, Aloisio et al., 2016).

Stochastic interactions are sampled through survival probabilities or optical depths. In one formulation, interaction points are drawn from

$A$3

where $A$4 is the total interaction rate. In another equivalent implementation, the code advances the cumulative interaction probability step by step and triggers an interaction once the integrated optical depth exceeds a random draw. Once an interaction occurs, the corresponding channel is selected probabilistically and the outgoing fragments are created with updated $A$5; each secondary is then propagated in turn through the same stack-based or branch-based loop (Batista et al., 2015, Aloisio et al., 2016, Aloisio et al., 2015).

Unstable products are usually treated as decaying instantaneously. This includes neutrons, pions, muons, and, in versions with $A$6-decay enabled, unstable nuclei. Neutrinos are then propagated trivially, with only adiabatic redshift losses, so that a neutrino produced with energy $A$7 at redshift $A$8 arrives with $A$9 (Aloisio et al., 2015).

3. Physical processes, cross sections, and photon backgrounds

The physical content of SimProp is organized around four main processes: adiabatic energy losses, Bethe-Heitler pair production, photonuclear disintegration of nuclei, and photopion production of nucleons and nuclei. The general interaction rate for a process with cross section $^{1}\mathrm{H}$0 in the nucleus rest frame is written as

$^{1}\mathrm{H}$1

with equivalent forms used in different versions for numerical tabulation and sampling (Batista et al., 2015, Aloisio et al., 2016).

Photodisintegration has been modeled in several ways. Early releases relied on the Puget-Stecker-Bredekamp (PSB) model, with one-nucleon and two-nucleon ejection cross sections parameterized as Gaussians in $^{1}\mathrm{H}$2 between threshold and $^{1}\mathrm{H}$3 MeV and a constant contribution from $^{1}\mathrm{H}$4 to $^{1}\mathrm{H}$5 MeV, while channels ejecting $^{1}\mathrm{H}$6-particles or heavier fragments were neglected in the default PSB setup. SimProp v2r3 added a simplified TALYS-derived “nucleon+$^{1}\mathrm{H}$7” scheme, in which effective nucleon-ejection and $^{1}\mathrm{H}$8-ejection cross sections are built from TALYS 1.6 yields and fitted by Gaussian-plus-plateau forms. Versions v2r3 and v2r4 also exposed Gaussian, Breit-Wigner, double-Breit-Wigner, and double-Gaussian options for user-supplied cross-section parameterizations (Batista et al., 2015, Aloisio et al., 2016, Aloisio et al., 2017).

Photopion production has likewise evolved. In v2r0, protons could be treated with a continuous energy-loss approximation. In v2r1 and later, stochastic treatments were introduced, based on total cross sections from SOPHIA for proton and neutron channels. A persistent simplifying assumption in SimProp is that photohadronic interactions are approximated by single-pion production with simplified center-of-mass kinematics and an isospin-based branching prescription, rather than the full multi-pion treatment implemented in CRPropa through SOPHIA (Batista et al., 2015, Batista et al., 2019).

The photon backgrounds are the CMB and a selectable EBL. The CMB is treated as a blackbody with $^{1}\mathrm{H}$9 K and exact FLRW redshift scaling. EBL support broadened substantially over time: early versions offered Stecker-type and Kneiske-type models, and v2r3/v2r4 added multiple implementations including Domínguez et al. 2011 and Gilmore et al. 2012. In SimProp v2r3, the provided $^{56}\mathrm{Fe}$0 grid is interpolated directly, without further scaling approximations, and the chosen EBL enters both photo-disintegration rates and, in later versions, photopion production on the EBL (Aloisio et al., 2013, Batista et al., 2015, Aloisio et al., 2017).

4. Version history, software structure, and outputs

The code evolved through a sequence of releases that progressively expanded its physics while retaining the one-dimensional Monte Carlo architecture.

Version or stage	Principal additions	Representative source
Early SimProp	One-dimensional UHECR propagation; CEL treatment plus stochastic photo-disintegration	(Aloisio et al., 2012)
v2r2	Full stochastic treatment of major interaction channels; cosmogenic neutrino production	(Aloisio et al., 2015)
v2r3	Expanded EBL choices and photodisintegration setups	(Aloisio et al., 2016)
v2r4	Secondary $^{56}\mathrm{Fe}$1 and $^{56}\mathrm{Fe}$2 production; reduced computation time	(Aloisio et al., 2017)

The v2r2 release was explicitly presented as treating all major interactions stochastically and taking into account the production of secondary cosmogenic neutrinos. The v2r3 release generalized the available EBL models and photodisintegration prescriptions. The v2r4 release introduced the recording of secondary cosmogenic particles such as electron-positron pairs and gamma rays produced during propagation, while also reporting a substantial reduction in computation time (Aloisio et al., 2015, Aloisio et al., 2016, Aloisio et al., 2017).

The implementation is in C++, and v2r4 is described as requiring a C++11-capable compiler and using the ROOT framework for I/O, random numbers, and data structures. Its high-level structure consists of an input parser, an event generator, a propagation engine with continuous and discrete process modules, a decay handler, a secondary writer, and an output writer. Depending on version, outputs are written either as ASCII event records or as ROOT trees. The v2r4 ROOT output includes the trees nuc, ev, and summary, with the compact summary tree storing injection parameters together with arrays for nuclei, photons, electrons, and neutrinos reaching Earth or produced during propagation (Aloisio et al., 2017).

User control is version-dependent. Some descriptions use an ASCII configuration file; later public versions are configured almost entirely through command-line switches controlling the injection mass number, energy range, redshift interval, EBL model, photodisintegration model, $^{56}\mathrm{Fe}$3-decay treatment, pion-production treatment, and the level of secondary-particle output (Aloisio et al., 2013, Aloisio et al., 2015, Aloisio et al., 2017).

5. Validation, uncertainties, and comparison with CRPropa

SimProp has been benchmarked against both semi-analytical kinetic solutions and other Monte Carlo codes. In the 2012 comparison with direct kinetic solutions for pure Fe injection, agreement was reported within $^{56}\mathrm{Fe}$4 for primaries and secondaries over $^{56}\mathrm{Fe}$5–$^{56}\mathrm{Fe}$6 eV. Comparisons with other Monte Carlo approaches gave broader but still structured agreement: with Allard et al. at the $^{56}\mathrm{Fe}$7–$^{56}\mathrm{Fe}$8 level, and with CRPropa 2 at the $^{56}\mathrm{Fe}$9–$z_{\max}$0 level, with the largest differences attributed to different nuclear-cascade treatments (Aloisio et al., 2012).

Later comparisons isolated residuals attributable to specific modeling choices. Batista et al. showed that, when SimProp v2r3 and CRPropa 3 are run with identical EBL and TALYS cross sections, the two codes differ by $z_{\max}$1 in the all-proton spectrum from $z_{\max}$2 to $z_{\max}$3 eV, corresponding to an effective $z_{\max}$4. For hard nitrogen injection, differences in $z_{\max}$5 and $z_{\max}$6 remain $z_{\max}$7, growing to $z_{\max}$8 above $z_{\max}$9 eV. The stated causes are SimProp’s simplified $$\frac{dN}{dE}\propto E^{-\gamma}\exp\!\left[-\frac{E}{Z R_{\rm cut}}\right].$$0 channel scheme, approximate redshift scaling of the EBL in CRPropa, and different light-nuclei cross sections (Batista et al., 2015).

The same work quantified the leading uncertainties internal to propagation modeling. Switching between Gilmore and Domínguez EBL models can produce up to $$\frac{dN}{dE}\propto E^{-\gamma}\exp\!\left[-\frac{E}{Z R_{\rm cut}}\right].$$1 differences in the all-particle spectrum at $$\frac{dN}{dE}\propto E^{-\gamma}\exp\!\left[-\frac{E}{Z R_{\rm cut}}\right].$$2 eV for hard-iron injection, with $$\frac{dN}{dE}\propto E^{-\gamma}\exp\!\left[-\frac{E}{Z R_{\rm cut}}\right].$$3 up to $$\frac{dN}{dE}\propto E^{-\gamma}\exp\!\left[-\frac{E}{Z R_{\rm cut}}\right].$$4 in the same scenario. Replacing PSB by TALYS nucleon-only photodisintegration changes the spectrum by less than $$\frac{dN}{dE}\propto E^{-\gamma}\exp\!\left[-\frac{E}{Z R_{\rm cut}}\right].$$5 except for intermediate secondaries, whereas including or rescaling $$\frac{dN}{dE}\propto E^{-\gamma}\exp\!\left[-\frac{E}{Z R_{\rm cut}}\right].$$6-channels can induce more than $$\frac{dN}{dE}\propto E^{-\gamma}\exp\!\left[-\frac{E}{Z R_{\rm cut}}\right].$$7 changes in the spectrum for hard-nitrogen scenarios and $$\frac{dN}{dE}\propto E^{-\gamma}\exp\!\left[-\frac{E}{Z R_{\rm cut}}\right].$$8. By contrast, transport approximations such as continuous versus stochastic pion production produce less than $$\frac{dN}{dE}\propto E^{-\gamma}\exp\!\left[-\frac{E}{Z R_{\rm cut}}\right].$$9 shifts in the proton spectrum above $z=0$0 eV and negligible composition changes (Batista et al., 2015).

Secondary messengers accentuate some of these sensitivities. Batista et al. found that cosmogenic neutrinos are more sensitive to the choice of EBL model than UHECRs, with flux ratios typically of order $z=0$1–$z=0$2 in the PeV–EeV range, depending on composition and source evolution. They also reported significant differences between neutrino fluxes predicted by the latest released versions of CRPropa and SimProp, tracing them mainly to the simplified EBL scaling in CRPropa and the single-pion/isospin approximation in SimProp (Batista et al., 2019).

These comparisons clarify SimProp’s recommended use cases. The code is described as ideal for rapid parameter scans, sensitivity studies of EBL and cross-section uncertainties, and first-order fits to the UHECR spectrum and $z=0$3. In applications demanding full photohadronic kinematics, detailed secondary yields, or three-dimensional magnetic deflections, a more comprehensive code such as CRPropa 3 is explicitly stated to be preferable (Batista et al., 2015).

6. Extensions, specialized modules, and scientific applications

A major application of SimProp has been the computation of cosmogenic neutrino fluxes. The v2r2 and related 2015 studies assembled all neutrinos produced in $z=0$4, $z=0$5, neutron, and unstable-nucleus decays, recording their production energy and redshift and constructing Earth fluxes after redshift correction. These studies used SimProp to compare source models against IceCube and Pierre Auger Observatory constraints, concluding that the available observations were already close to constraining source composition and cosmological evolution (Aloisio et al., 2015, Aloisio et al., 2015, Aloisio et al., 2015).

The v2r4 release further expanded SimProp into a multi-messenger generator by recording secondary $z=0$6 pairs and gamma rays for subsequent treatment by external cascade codes. Batista et al. then used SimProp together with CRPropa to quantify how uncertainties in the EBL and photodisintegration cross sections propagate into cosmogenic neutrino and photon predictions, finding that overall cosmogenic gamma-ray production rates are relatively independent of propagation details even when neutrino fluxes are not (Aloisio et al., 2017, Batista et al., 2019).

The one-dimensional, no-magnetic-field assumption has also been relaxed in specialized extensions. In 2021, a modified SimProp v2r4 core was coupled to a stochastic differential equation solver that updates the propagation direction $z=0$7 and comoving displacement $z=0$8 in a turbulent extragalactic field. In that framework, the diffusion length is expressed through a fitted Kolmogorov-spectrum form, and the resulting spectral effects are summarized by multiplicative suppression factors $z=0$9 and $\Delta z$0, enabling magnetic-field-modified spectra to be reconstructed from field-free SimProp outputs (González et al., 2021).

Another extension moved SimProp upstream, into source environments. Condorelli et al. embedded a hadronic module for starburst nuclei, treating $\Delta z$1, $\Delta z$2, and nucleus-nucleus spallation via SIBYLL 2.3d together with photo-hadronic interactions on source photon fields, before handing escaping fragments to the standard extragalactic propagation engine. In this setting, spallation reactions could exceed the outcome in neutrinos from photo-hadronic interactions in the source environment and in the extra-galactic space, depending on the source gas density (Condorelli et al., 2022).

SimProp has also served as a base for beyond-standard-propagation studies. One extension introduced Lorentz-invariance-violating neutrino propagation, replacing the default neutrino module with a decay-while-propagating treatment for vacuum pair emission and neutrino splitting, while leaving the cosmic-ray propagation modules unchanged. A related 2026 study interfaced SimProp-generated cosmogenic neutrino fluxes with a flavor-transition calculation in the presence of Lorentz invariance violation, using the resulting ultra-high-energy $\Delta z$3-neutrino fluxes to forecast sensitivities for GRAND and POEMMA (Reyes et al., 2023, Brdar et al., 21 Apr 2026).

The code’s nuclear scope has likewise been revised. A 2025 technical report described minimal modifications to SimProp v2r4 that extend the nuclide list from $\Delta z$4 to $\Delta z$5, enlarge the internal array bounds, replace the approximate pair-production scaling by explicit $\Delta z$6, and fit TALYS 2.0 cross sections for $\Delta z$7. In parallel, Soriano et al. proposed an update of v2r4 light-nucleus photo-disintegration modules using new parametrizations for $\Delta z$8He, $\Delta z$9He, tritium, and deuterium. In their numerical example, the present-day mean free path of ultra-relativistic $\Delta t$0He increases by $\Delta t$1 relative to the default PSB-based SimProp routine, implying a $\Delta t$2 larger survival probability for $\Delta t$3 GeV over $\Delta t$4 Mpc (Matteo, 4 Sep 2025, Soriano et al., 2018).

Taken together, these developments show that SimProp is best understood not as a single immutable transport code, but as a compact Monte Carlo framework whose central identity—fast, one-dimensional UHECR propagation with configurable interaction physics—has supported a broad range of precision studies in UHECR phenomenology, cosmogenic neutrinos, gamma-ray secondaries, magnetic-horizon effects, source-environment processing, and non-standard neutrino propagation (Batista et al., 2015).