CRPropa 3.3: Multi-Messenger Simulation

• CRPropa 3.3 is a Monte Carlo simulation framework that models cosmic rays, gamma rays, neutrinos, and electromagnetic cascades across astrophysical settings.
• It features a modular, candidate-based design integrating source modeling, transport, interactions, and observable predictions using C++ and Python with parallel processing.
• Recent enhancements include explicit time tracking, spatially varying radiation fields, and an expanded nuclear database up to lead, broadening its astrophysical applications.

CRPropa 3.3 is a publicly available Monte Carlo framework for simulating the propagation of high-energy particles in astrophysical environments. It is designed to model cosmic rays, gamma rays, neutrinos, and secondaries in magnetic fields and radiation backgrounds, and the 3.3 release is framed as a step toward a unified multi-messenger simulation platform spanning GeV to ZeV energies. In practice, the code combines source modeling, transport, interactions, cascades, and observable prediction within a single modular environment, extending the CRPropa 3 lineage from ultra-high-energy cosmic-ray propagation to a broader high-energy astroparticle scope (Aerdker et al., 15 Jul 2025).

1. Framework lineage and software architecture

CRPropa 3.3 belongs to the CRPropa 3 family, whose defining architectural change was the move from a monolithic propagator to a modular, candidate-based framework. A simulation is assembled as a user-defined sequence of independent modules acting on a Candidate object; the candidate stores particle state, and each module may modify it until a stopping condition or detection is reached. This architecture was introduced as the foundation of CRPropa 3 and remained the core design principle for later minor releases, including 3.3 (Batista et al., 2016).

The framework is written in C++ and exposed to Python via SWIG, with cross-language polymorphism so that users can implement custom modules in Python and combine them with C++ modules. Shared-memory multiprocessing with OpenMP was introduced to exploit the fact that cosmic-ray propagation is embarrassingly parallel at the level of independent particles. Public development infrastructure is part of the framework’s stated design: a public git repository, issue tracker, peer-reviewed contributions, unit tests, documentation, and examples are explicitly emphasized in later overview papers (Batista et al., 2014, Batista et al., 2019).

A consequence of this design is extensibility beyond the original UHECR use case. RadioPropa, for example, was built as a fork of CRPropa that reused the geometry, boundary handling, data containers, file I/O, and modular “module chain” architecture, while replacing the particle-transport physics with ray-tracing modules for in-matter radio propagation. That reuse was presented as evidence that CRPropa already provided a mature modular backbone, even when a separate fork was justified by nontrivial changes to the central data structure and stepping logic (Winchen, 2018).

2. Propagation regimes and transport formalisms

CRPropa supports multiple propagation regimes because different energy ranges and environments require different numerical descriptions. The classical mode is ballistic single-particle propagation, aimed primarily at extragalactic cosmic rays and ultra-high-energy particles, in which individual particles are followed through arbitrary three-dimensional magnetic fields with interactions, decay, and adiabatic losses included. Earlier CRPropa 3 papers emphasized 1D, 3D, and 4D propagation, with 4D introduced specifically so that cosmological effects could be included consistently when magnetic deflections make the effective path length unknown a priori (Batista et al., 2014, Batista et al., 2016).

The galactic regime required a different formalism. Recent versions complement ballistic tracking with a solver for the differential transport equation, implemented through a stochastic equivalent in which pseudo-particles undergo effective displacements in space and momentum that reproduce the transport-equation solution in the many-particle limit. The modified Parker transport equation is explicitly given in the CRPropa literature as

$\frac{\partial n}{\partial t} + \vec{u}\cdot\nabla{n} = \nabla\cdot \left(\hat{\kappa}\nabla{n}\right) + \frac{1}{p^2}\frac{\partial}{\partial p}\left(p^2 \kappa_{pp} \frac{\partial n}{\partial p} \right) + \frac{1}{3}\left(\nabla \vec{u}\right)\frac{\partial n}{\partial \ln{p} + S\left(\vec{x},p,t\right)$

with diffusion anisotropy often characterized by

$\epsilon = \frac{\kappa_\perp}{\kappa_\parallel}.$

This transport machinery was presented as a major extension of CRPropa beyond the original extragalactic UHECR focus (Batista et al., 2019).

By the 3.2 generation, CRPropa was described as covering energies from ZeV down to GeV for gamma rays and electrons, and TeV for cosmic rays and neutrinos, across source environments as well as Galactic and extragalactic propagation domains. The 3.3 paper retains that broad astrophysical scope but emphasizes structured, time-dependent environments more strongly than earlier versions (Aerdker et al., 2023, Aerdker et al., 15 Jul 2025).

3. Interaction physics and multimessenger capabilities

The physical content of CRPropa was originally centered on ultra-high-energy cosmic rays interacting with the cosmic microwave background and extragalactic background light. The framework includes pair production, photo-meson production, photodisintegration, nuclear decay, and adiabatic energy losses from the expansion of the Universe. These processes determine the UHECR spectrum, composition, and arrival directions, and also generate secondary neutrinos and photons, making CRPropa intrinsically a multimessenger code rather than a cosmic-ray-only propagator (Vliet, 2016, Batista et al., 2014).

Electromagnetic secondaries became progressively more central in the CRPropa 3 line. Early CRPropa 3 descriptions relied on external tools such as DINT and EleCa for photon propagation and electromagnetic cascades, with DINT solving one-dimensional cascade transport equations and EleCa providing a Monte Carlo treatment at high energies. Subsequent releases brought photon propagation into the core framework, adding all major energy-loss processes present in DINT and EleCa together with additional photon production channels from photo-disintegration of nuclei, elastic scattering, and radiative decay, plus magnetic deflections of local electromagnetic secondaries (Batista et al., 2016, Heiter et al., 2017, Batista et al., 2019).

This transition materially broadened the observable space CRPropa could address. The code line now supports gamma rays, photons, electrons/positrons, neutrinos, and secondaries from hadronic and nuclear interactions in a common simulation environment. In the CRPropa 3.1 and 3.2 literature, the full electromagnetic cascade chain includes pair production and double pair production for photons, triplet pair production and inverse Compton scattering for electrons and positrons, and synchrotron radiation in magnetic fields. The result is a framework capable of following the coupled evolution of hadronic and electromagnetic messengers across the TeV to sub-EeV bridge as well as in the UHECR regime (Heiter et al., 2017, Batista et al., 2021).

4. Distinctive additions in CRPropa 3.3

The central technical additions in CRPropa 3.3 are explicit time tracking, time-dependent advection fields, position-dependent radiation backgrounds, and an expanded nuclear database up to lead. The new time variable is implemented as an explicit candidate property, incremented at each propagation step, with corresponding outputs and units added. The current implementation still assumes $v=c$, so the explicit time handling is presently formulated for highly relativistic particles; the 3.3 paper notes that relaxing this could enable treatment of lower-energy particles below roughly GeV energies in the future (Aerdker et al., 15 Jul 2025).

Time dependence is also extended to the advection framework, whose velocity field can now depend on both space and time. This was introduced to support use cases such as diffusive shock acceleration at moving shocks, where the local flow pattern evolves and the usual Green’s-function method for time-dependent solutions cannot be used. The paper explicitly notes that such cases require more pseudo-particles, increasing computational cost but avoiding artifacts from time integration and pseudo-particle recycling (Aerdker et al., 15 Jul 2025).

A second major addition is support for spatially varying photon fields through the PhotonBackground and InteractionRates infrastructure. Rather than treating only homogeneous backgrounds such as the CMB, CRPropa 3.3 can now use fields specified on spatial grids of spectral intensities, notably the interstellar radiation field. From these grids, the code computes inverse mean free paths for electromagnetic processes, selecting the nearest grid node to the particle position, with nearest-node search accelerated using the nanoflann library. The stated motivation is accurate modeling of gamma-ray attenuation and cascades in inhomogeneous environments such as the Milky Way, galaxy clusters, and superbubbles (Aerdker et al., 15 Jul 2025).

The nuclear update is equally significant. CRPropa 3.3 adds cross sections for nuclei up to lead, $Z=82$, compiled previously and generated using TALYS 1.9. The paper states that this increases the number of available nuclear species by more than a factor of five relative to the previous version. Two explicit thinning criteria are used: only nuclei with lifetimes exceeding $100\,\mathrm{s}$ are included, and for each nucleus only disintegration channels contributing at least $0.1\%$ to the energy-weighted average cross section are retained, preserving at least $90\%$ of the unthinned total. This reduces the number of channels from

$311{,}403 \to 16{,}956$

which is less than $6\%$ of the original set (Aerdker et al., 15 Jul 2025).

The nuclear interaction update is expressed through the energy-weighted average cross section

$$\tilde \sigma(\varepsilon) = \frac{2}{\varepsilon^2} \int_0^{\varepsilon} \varepsilon' \sigma(\varepsilon') \, \mathrm{d}\varepsilon'$$

and the interaction rate in a photon background

$\epsilon = \frac{\kappa_\perp}{\kappa_\parallel}.$0

The paper emphasizes that $\epsilon = \frac{\kappa_\perp}{\kappa_\parallel}.$1 is independent of any specific target photon field and depends only on the cross section as a function of photon energy in the nucleus rest frame, which is why the reduced channel set remains applicable across different photon spectra (Aerdker et al., 15 Jul 2025).

CRPropa 3.3 also introduces additional environment models and workflow improvements. Two updated Galactic magnetic-field implementations are added: UF23Field, which implements the eight models by Unger & Farrar (UF24), and the KST24 halo field by Korochkin et al. New EBL models include IRB_Saldana and IRB_Finke22. The framework can now store the full state of all particles in a simulation chain and restart from a stopped point, a feature explicitly motivated by long-running simulations and wall-time limits on large computing facilities (Aerdker et al., 15 Jul 2025).

5. Representative applications

The 3.3 paper highlights two showcase applications that directly exploit the new features. The first is gamma-ray propagation in the Milky Way with spatially dependent interstellar radiation field models. Two ISRF implementations are discussed, F98 and R12, and the example of the first LHAASO catalog source 1LHAASOJ1825-1256u is used to compare propagated spectra including Galactic photon backgrounds with a naive extrapolation that ignores propagation effects. The stated conclusion is that above tens of TeV, gamma rays can be significantly suppressed by pair production on the ISRF, which is precisely the class of effect that homogeneous-background approximations cannot capture (Aerdker et al., 15 Jul 2025).

The second showcase is diffusive shock acceleration at a moving shock, enabled by time-dependent advection. The example uses the OneDimensionalTimeDependentShock advection field with parameters $\epsilon = \frac{\kappa_\perp}{\kappa_\parallel}.$2, $\epsilon = \frac{\kappa_\perp}{\kappa_\parallel}.$3, $\epsilon = \frac{\kappa_\perp}{\kappa_\parallel}.$4, and $\epsilon = \frac{\kappa_\perp}{\kappa_\parallel}.$5, with the shock profile smoothed by a hyperbolic tangent. The flow velocities are required to satisfy the Rankine-Hugoniot conditions in the shock frame, with upstream velocity $\epsilon = \frac{\kappa_\perp}{\kappa_\parallel}.$6, downstream velocity $\epsilon = \frac{\kappa_\perp}{\kappa_\parallel}.$7, and compression ratio $\epsilon = \frac{\kappa_\perp}{\kappa_\parallel}.$8; for a strong shock,

$\epsilon = \frac{\kappa_\perp}{\kappa_\parallel}.$9

The example uses lab-frame values $v=c$0 and $v=c$1, corresponding to a shock-frame compression consistent with $v=c$2 (Aerdker et al., 15 Jul 2025).

Beyond the 3.3 paper, the CRPropa 3 line has been used for multimessenger inference in UHECR source modeling. One widely cited application used the newest CRPropa 3 generation to predict the UHECR spectrum at Earth together with cosmogenic neutrino and photon fluxes produced during propagation through the CMB and EBL. When compared with the Fermi-LAT isotropic diffuse gamma-ray background, especially the $v=c$3–$v=c$4 bin, pure proton source models were found to be strongly constrained unless their source density was strongly peaked at recent times. In that study the gamma-ray channel was explicitly more constraining than IceCube neutrino data for the models considered (Vliet, 2016).

Other application domains were emphasized in the 3.2 and 3.3 overviews: Galactic cosmic-ray density evolution, second-order Fermi acceleration that can generate peaking spectra, diffusive shock acceleration with stochastic differential equations, and gamma-ray propagation in intergalactic and Galactic radiation fields. Taken together, these examples indicate that CRPropa’s typical workflow is no longer confined to source-to-Earth UHECR tracking, but increasingly includes acceleration, transport in structured environments, and secondary-emission prediction within the same framework (Batista et al., 2019, Aerdker et al., 2023).

6. Approximations, validation, and known limitations

CRPropa 3.3 is not approximation-free, and several limitations are explicit in the literature. At the framework level, CRPropa is restricted to the highly relativistic regime, and the new explicit time handling in 3.3 still assumes $v=c$5. This is sufficient for the intended high-energy domain but leaves lower-energy nonrelativistic transport outside the implemented scope. The position-dependent photon-field treatment also uses the nearest grid node rather than a more elaborate interpolation scheme for the local spectral intensity in the implementation described by the 3.3 paper (Batista et al., 2019, Aerdker et al., 15 Jul 2025).

A more substantive limitation appears in precision very-high-energy gamma-ray propagation at moderate and high redshift. A code-comparison study involving CRPropa found that, after correcting electromagnetic interaction-table issues and adjusting magnetic-field scaling, CRPropa can model nearby-source cascades reasonably well: for nearby sources with $v=c$6 and zero IGMF the inter-code discrepancy falls to about $v=c$7, and with nonzero IGMF it remains generally below $v=c$8. However, the same study concluded that CRPropa remains unreliable for precise modeling of spectral, timing, and imaging properties of secondary gamma rays from sources approaching $v=c$9, because its simplified EBL evolution uses a fixed spectral shape with redshift-dependent renormalization. In that regime the disagreement can reach order-of-magnitude levels, especially at the high-energy end of the cascade spectrum (Kalashev et al., 2022).

Secondary-neutrino predictions exhibit a different but related sensitivity. A comparison of SimProp and CRPropa found that when the same EBL and photodisintegration inputs are used, the two codes agree broadly for primary UHECR spectra and composition, but not for neutrinos: CRPropa predicts about $Z=82$0 higher EeV neutrino fluxes than SimProp, and in proton injection scenarios it can exceed SimProp by several times at PeV energies. The paper attributes most of the discrepancy to CRPropa’s approximate EBL redshift treatment and to SimProp’s simplified photohadronic branching, illustrating that apparently modest implementation differences can become important for secondary-messenger predictions (Batista et al., 2019).

The broader uncertainty literature around the CRPropa 3 line also shows that some of the largest variations arise not from numerical bugs but from physical input choices exposed by the framework. In extragalactic nuclei propagation, the EBL model and photodisintegration setup, especially the treatment of $Z=82$1-emission channels, were found to have a considerable impact on the simulated UHECR spectrum and composition. For example, differences between the Gilmore and Domínguez EBL models can exceed $Z=82$2 in hard-injection scenarios, and varying the $Z=82$3-ejection factor in TALYS-based modeling can change the hard-nitrogen spectrum by more than $Z=82$4. This suggests that CRPropa’s value lies as much in controlled model comparison as in nominal forward prediction (Batista et al., 2015).

These caveats do not negate the framework’s central role. Rather, they define the conditions under which CRPropa 3.3 should be interpreted: as a comprehensive modular platform for high-energy astroparticle transport whose predictions inherit both the strengths and the limitations of the interaction models, background-field prescriptions, and transport approximations selected by the user (Aerdker et al., 15 Jul 2025).