A new view on Auger data and cosmogenic neutrinos in light of different nuclear disintegration and air-shower models

We study the implications of Ultra-High Energy Cosmic Ray (UHECR) data from the Pierre Auger Observatory for potential accelerator candidates and cosmogenic neutrino fluxes for different combinations of nuclear disintegration and air-shower models. We exploit the most recent spectral and mass composition data (2017) with a new, computationally very efficient simulation code PriNCe. We extend the systematic framework originally developed by the Pierre Auger Collaboration with the cosmological source evolution as an additional free parameter. In this framework, an ensemble of generalized UHECR accelerators is characterized by a universal spectral index (equal for all injection species), a maximal rigidity, and the normalizations for five nuclear element groups. We find that the 2017 data favor a small but constrained contribution of heavy elements (iron) at the source. We demonstrate that the results moderately depend on the nuclear disintegration (PSB, Peanut, or Talys) model, and more strongly on the air-shower (EPOS-LHC, Sibyll-2.3, or QGSjet-II-04) model. Variations of these models result in different source evolutions and spectral indices, limiting the interpretation in terms of a particular class of cosmic accelerators. Better constrained parameters include the maximal rigidity and the mass composition at the source. Hence, the cosmogenic neutrino flux can be robustly predicted, since it originates from interactions with the cosmic infrared background and peaks at $10^8 \, \mathrm{GeV}$. Depending on the source evolution at high redshifts the flux is likely out of reach of future neutrino observatories in most cases, and a minimal cosmogenic neutrino flux cannot be claimed from data without assuming a cosmological distribution of the sources.


Introduction
The two largest detectors ever built, the Pierre Auger Observatory (Aab et al. 2015) and the Telescope Array (Abu-Zayyad et al. 2013), investigate the origin and nature of ultrahigh-energy cosmic rays (UHECRs) above 10 18 eV with hybrid detection techniques that combine signals from surface and fluorescence detectors to reconstruct extensive air showers, which are giant particle cascades initiated through interactions of the UHECRs with the atmosphere. There is evidence for an extragalactic origin of the UHECRs (Aab et al. 2017c), and studies of the UHECR arrival directions uncovered interesting patterns, such as a strong dipole anisotropy and a correlation with nearby source directions (Aab et al. 2018). However, an association with a concrete source or class of sources is not yet in reach. The chemical composition is likely to be a mixture of different nuclear masses (Aab et al. 2016a), ranging from protons up to nitrogen or heavier nuclei (Aab et al. 2017b). While the mass-sensitive experimental observables are statistically in agreement between the two experiments, their interpretation in terms of physical mass composition is still subject to discussion (de Souza 2018).
Various astrophysical phenomena, typically associated with the emission of high-energy photons, have been proposed as potential accelerators of UHECRs. Gamma-ray bursts (GRBs), provided that a significant fraction of baryons is accelerated in their jets, can be capable of emitting UHECRs and also producing high-energy neutrinos due to photohadronic interactions of protons or heavier nuclei with the target photons (Waxman & Bahcall 1997). Blazars, a subset of powerful active galactic nuclei with their jets pointing at the observer, are numerous and powerful enough to sustain the UHECR spectrum and have been considered as sources of UHECRs and high-energy neutrinos (Stecker et al. 1991;Murase et al. 2014;Rodrigues et al. 2018). The absence of an associated neutrino signal in the IceCube detector (Aartsen et al. 2017a(Aartsen et al. , 2017b constrains the density of cosmic rays in GRBs and blazars but does not necessarily exclude these classes of sources as UHECR accelerators. Other compact source classes, such as jetted tidal disruption events (TDEs; Farrar & Piran 2014) or low-luminosity GRBs (LL-GRBs; Murase et al. 2006), are potentially luminous or copious enough to power the UHECR and high-energy neutrino sky. Starburst galaxies constitute a sample of sources in which the reacceleration of PeV cosmic rays to ultrahigh energies may occur at the termination shocks of kpc-scale "super winds" (Anchordoqui et al. 1999). A higher abundance of young pulsars (Blasi et al. 2000) as an effect of an enhanced supernova rate might also predestine these galaxies as hosts of UHECR accelerators. The anisotropy observed by the Pierre Auger Observatory indeed indicates a directional correlation with a subset of nearby gamma-ray-bright starburst galaxies (Aab et al. 2018). In all cases, the direct association with highenergy neutrinos would be a smoking-gun signature for the origin of the cosmic rays. If, on the other hand, the neutrino production in the sources is inefficient, a directly related neutrino signal will be absent, and indirect methods will be needed to infer the nature of the cosmic-ray accelerators.
Obtaining information on the distribution of sources (such as their evolution as a function of redshift) is one such indirect method to identify the accelerators and will, therefore, be one of the main targets of our study.
The identification of the UHECR sources is complicated by the transport through the intergalactic medium (IGM), where interactions with the cosmic microwave background (CMB) and cosmic infrared background (CIB) photons alter the spectrum and chemical composition compared to the original emission at the source. By assuming a model for the UHECR spectra emitted from the sources and the extragalactic propagation through the IGM, one can infer the free source model parameters through a fit to the available UHECR data. In several such studies (Hooper & Taylor 2010;Aloisio et al. 2014;Globus et al. 2015b;Taylor et al. 2015;Aab et al. 2017a;Wittkowsi 2017), it has been assumed that the sources are identical and isotropically distributed and that the UHECR emission follows power-law spectra with a rigidity-dependent cutoff. Since these sources are representing generic accelerators, the cosmological evolution of the source density is undefined and requires one or multiple additional free parameters. Typically, one assumes piecewise defined evolution functions of the form z 1 m + ( ) , with m the evolution parameter. Due to accumulation of energy losses over large distances, UHECRs, even without considering magnetic fields, experience a horizon or maximal distance they can travel through the IGM, which is approximately equivalent to a redshift of z 1 , or a few Gpc. Therefore, the UHECR spectrum is almost insensitive to the parameterization of the source evolution beyond redshift z 1 . Interactions of UHECRs leave traces, namely, cosmogenic neutrinos that are produced in photohadronic interactions with the target photons. Since neutrinos travel unimpeded through the IGM, the density of UHECRs for z 1 > has an impact on their flux. As a consequence, the cosmogenic neutrino flux can be used to constrain the cosmological source evolution (Ahlers et al. 2009;Gelmini et al. 2012;Aloisio et al. 2015;Heinze et al. 2016;Romero-Wolf & Ave 2018;Das et al. 2018;Møller et al. 2018;Wittkowski & Kampert 2018;Alves Batista et al. 2019b;van Vliet et al. 2019).
The modeling of the transport comes with a number of uncertainties: photonuclear (photodisintegration) reactions (Alves Batista et al. 2015;Boncioli et al. 2017;Soriano et al. 2018) that change the mass composition of nuclei due to interactions with CMB or CIB photons; the hadronic interactions, which are used in the interpretation of air-shower observables in terms of the mass composition; and the CIB spectrum, which is not well known at high redshifts. The interpretation of the UHECR data is affected by these uncertainties, as demonstrated in Alves Batista et al. (2015) and in the combined fit (CF) of the spectrum and composition data by the Pierre Auger Collaboration (Aab et al. 2017a).
While in the CF, different assumptions for source density evolution have been tested for compatibility, no conclusions have been drawn about possible associations with sources. Hence, most attention was devoted to a flat cosmological evolution (nonevolving source densities; Aab et al. 2017a), which, however, cannot be easily related to known accelerator candidates. As an example, sources can evolve similarly to the star-forming rate (SFR), z 1 3.4 + ( ) , for z 1 < , such as GRBs (Wanderman & Piran 2010). Blazars have a typically more complicated luminosity-dependent evolution function and can evolve more steeply with redshift. Some source classes, such as TDEs, may have negative source evolution. As a consequence, any attempt to seek an astrophysical interpretation within the framework of such a fit requires the source evolution to be a free parameter. However, each new parameter is computationally expensive, which has led to different strategies to deal with this problem; for example, the redshift evolution can be included in a coarser way (Alves Batista et al. 2019b) or in a limited range of values (Romero-Wolf & Ave 2018; see also Das et al. 2018Møller et al. 2018; for similar studies).
In this paper, we revisit the approach of the CF, taking into account the dominant model dependencies, and focus on the degeneracies between the fit parameters given a homogeneous distribution of generic UHECR sources. We study the impact of the model uncertainties on the astrophysical interpretation by performing scans in the three parameters, maximum rigidity R GV max [ ] (corresponding to the maximum energy of acceleration divided by the charge of the particle, E Z max ), spectral index γ, and cosmological evolution index m, using different combinations of nuclear disintegration and air-shower models. The computational requirements are significantly reduced through the new numerical code Propagation including Nuclear Cascade equations (PRINCE) that performs the propagation very efficiently under changing physical conditions. We are, therefore, able to investigate the full 3D source parameter space with a comparable resolution in all parameters for different nuclear disintegration models. With Monte Carlo or slower numerical codes, such a study is not feasible due to excessive requirements of computational resources, and thus our result is novel. As an important result, we obtain the allowed parameter space contours that represent the state of the art of current UHECR observations. Under the assumption of one dominant source population that accelerates cosmic-ray nuclei up to a maximal rigidity, we accurately compute the expected cosmogenic neutrino fluxes and discuss the robustness of the predictions by studying the major model uncertainties.

Models of UHECR Transport and Their Sources
In this section, we describe the main model uncertainties affecting our analysis: the photobackgrounds and cross sections for the interactions during propagation, the hadronic interaction models used to infer UHECR properties from the observed air showers, and the implied assumptions about the distribution and characteristics of UHECR sources.

Extragalactic Propagation
During extragalactic propagation, UHECRs interact with the CMB and CIB via photopair (e + e − ) production and photonuclear processes. Additionally, all relativistic particles lose energy adiabatically due to the expansion of the universe. Photonuclear interactions can be subdivided into two regimes: photodisintegration ( 150 MeV r e < ) and photomeson production (above the pion production threshold, 150 MeV r e > ), where r e is the photon energy in the nuclear rest frame. In the photodisintegration regime, the target photons interact with one or two nucleons and collectively excite the nucleus into a resonant state, which subsequently decays emitting (evaporating) nucleons, heavier fragments, or keV-MeV photons. To model the cascading of secondary nuclei during propagation, numerical codes, such as PRINCE, described in Section 3, or Monte Carlo packages, require as input inelastic interaction cross sections and inclusive cross sections (or multiplicities) of secondary particles. Such cross sections can be obtained either empirically from data, as in the Puget-Stecker-Bredekamp (PSB;Puget et al. 1976) parameterization, or by tabulating the output of more realistic nuclear models. In this study, we use TALYS (Koning et al. 2007), a comprehensive pre-equilibrium and Hauser-Feshbach theory-based code, and PEANUT (Fassò et al. 1997(Fassò et al. , 2005, an event generator of the FLUKA package Böhlen et al. 2014) with an intranuclear cascade model at energies 200 MeV r e > and a similar set of statistical models below that (see Boncioli et al. 2017 for a discussion of these models and their uncertainties).
Qualitatively, the distributions of secondaries are similar for the two statistical models, while quantitatively, the results may vary depending on the availability of data for each individual isotope and the degree of parameter optimization for each of these isotopes. We observe that in the default configuration, PEANUT is better optimized to the available data. Unofficial tables for TALYS are available that can improve the description for some isotopes (Alves Batista et al. 2015). Compared to the PSB parameterization, where only one isotope for each mass number is used, PEANUT and TALYS demonstrate a faster disintegration into lighter elements, including the presence of heavier fragments (D, T, 3 He, 4 He). Therefore, the interpretation of the UHECR data in terms of composition at the source is expected to vary with respect to the use of different disintegration models.
Pion production off nuclei in all current propagation codes is handled in a "superposition" approach; i.e., the nucleons are treated as quasi-free. with the number of protons Z and the number of neutrons N.
The dominant pion production process is the Δ-resonance production in the s-channel, p pn The pion takes about 20% of the primary's energy and results in significant energy losses for the projectile. In absence of other processes, the cutoff in the UHECR spectrum at E 4 10 GeV 10 » · could be attributed to this energy loss, as predicted in Greisen (1966) and Zatsepin & Kuzmin (1966) and referred to as the Greisen-Zatsepin-Kuzmin (GZK) cutoff. In the case of nuclei, the Δ-resonance threshold is shifted by A to higher energies. Instead, most interactions take place at the energies of the giant dipole resonance around 20 MeV r e~, leading to a cutoff in the spectrum of UHECR nuclei at energies similar to the GZK cutoff.
As cosmogenic neutrinos are only produced in the photomeson regime, the differences between free nucleons and nuclei are striking. The photodisintegration threshold prevents nuclei from reaching energies A 10 GeV 10 > · , where photomeson production sets in on CMB target photons. Instead, pions and cosmogenic neutrinos are produced by nuclei at energies below the cutoff 10 GeV 9 on the less abundant CIB target photons. There are two consequences: the neutrino flux is expected to peak at lower energies, 10 GeV 8 , and be significantly lower compared to the protons-on-CMB case. The impact of CIB variations on UHECR propagation has been studied in Alves Batista et al. (2015) and Aab et al. (2017a).
While the effect on UHECR spectra is small, it becomes sizable for cosmogenic neutrino fluxes (see, e.g., Aloisio et al. 2015).
Extragalactic and galactic magnetic fields play an important role at the ankle, which is the change of the spectral index at 5 10 GeV 9 · (Fenu 2017), and below. The curvature of the UHECR trajectories effectively elongates the distance to the sources. At sufficiently low rigidities ( 10 18  V), the particles are increasingly trapped in the neighborhood of their accelerator. The quantitative impact has been studied, for example, in Mollerach & Roulet (2013). It results in a hardening of the individual spectra of nuclei at lower energies at Earth and thus can soften the spectral index required at the source. In this work, we neglect the effect of the magnetic fields, assuming a purely ballistic treatment of UHECR transport, as, for example, in Allard et al.

Air-shower Model
When cosmic-ray nuclei enter the atmosphere, the inelastic interactions with air molecules create hadronically (mesons and baryons) and electromagnetically (e ± and photons) interacting particles with smaller energies. This cascading proceeds until most of the initial energy is dissipated as light and long-lived particles (see, e.g., Matthews 2005 for an instructive model). The observation of the light and the secondary particles from these so-called extensive air showers allows the reconstruction of several properties of the original particle, such as the energy, direction, and, to some extent, mass composition (see Kampert & Unger 2012 for a review). At the Pierre Auger Observatory and the Telescope Array, the energy is measured calorimetrically through the integration of the total fluorescence light yield. The direction is inferred through stereoscopy in combination with timing-based measurements at the ground. The nuclear mass of the UHECR is the most challenging property, since it can only be derived indirectly by comparing a large number of observations with model-dependent simulations. Hence, the measurement of the composition is a statistical argument.
The sensitive variable for the mass composition is X E max ( ), the depth at which the energy dissipation of a single air shower is maximal. The X max fluctuates, since the first interaction statistically occurs at different altitudes and because secondary particles can be produced with a multitude of multiplicity and energy configurations. The simplest description that captures the observed distributions is the combination of the mean X max á ñ and the dispersion or variance X max s ( ). The expected values are shown in Figure 1, together with expectations for individual nuclei obtained with different interaction models.
Our simulations of the UHECR transport produce individual spectra for each nuclear mass at the top of the atmosphere for which we compute A ln á ñ and A ln 2 s at each energy of the numerical grid. We exactly follow the procedure from Abreu et al. (2013;Section2) to convert the average of the logarithmic mass and its dispersion (bin-wise in energy) to the experimental observables X max á ñ and X max s ( ) using where X max p á ñ is the mean depth at maximum of the proton showers and f E parameterizes the dependence on the air-shower model and energy, and linearly depends on the dispersion of the masses. All parameters are dependent on the logarithm of the cosmic-ray energy. The values of the parameters are obtained from air-shower simulations that do not take detector effects into account. Instead, this is taken into account by comparing with observables that are already corrected for detector effects. In contrast to the original paper (Abreu et al. 2013), we use an updated set of parameters for the post-LHC interaction models. 2 Essentially, the first moment X max á ñ has a linear dependence on A ln where some nonlinear effects are absorbed in f E . When fitting the data, the different model expectations for X max p á ñ impose shifts of the A ln á ñ that are the result of the propagation simulation and its initial conditions. The second term of the dispersion X 2 max s ( ) becomes small if only a single mass is present or if spectra of similar/neighboring masses are superimposed. It is large in cases where a few masses with large distance in A ln dominate the sum of the spectra. The simultaneous description of both the mean and the variance of X max is indispensable for any serious interpretation of the composition results, since the variables are supplementary and sensitive to different features of the UHECR flux.
For the present study, the differences in the conversion between mass and X max observations are the most relevant feature of Figure 1. For example, at a fixed X max á ñ, the A ln á ñ inferred with SIBYLL 2.3 is heavier compared to the other models. At the same time, the shower-to-shower fluctuations sh 2 s á ñ in Equation (3) are high, implying strong constraints for the mass dispersion term A ln 2 s . While one can simply say "SIBYLL 2.3 is heavier" than EPOS-LHC, the pulls on the fit induced by the properties of the models are highly nontrivial and discussed in a more "applied" way in Section 5.2. Note that some models, like QGSJETII-04, fail to produce a consistent relation between mass and X max variables (Aab et al. 2014;Bellido 2017).

Source Model
Several source candidates-in particular, compact jetted sources, such as GRBs (Globus et al. 2015a;Biehl et al. 2018a;Zhang et al. 2018;Boncioli et al. 2019) or TDEs (Zhang et al. 2017;Biehl et al. 2018b;Guépin et al. 2018)-can describe the UHECR spectrum and composition. Another category of viable UHECR sources are starburst galaxies (Anchordoqui et al. 1999;Anchordoqui 2019) that may also contain populations of powerful accelerators (Fang et al. 2013). The majority of models assume Fermi acceleration as the dominant acceleration process, yielding a power law with spectral indices close to 2 g = at the acceleration site. Hence, charged particles are magnetically confined at the site of acceleration, leading to an additional modification of the spectrum due to the escape mechanism. For example, diffusive or direct escape hardens the in-source flux by up to one power (Baerwald et al. 2013), while advective escape may act as a low-pass filter and suppress the high-energy emission in the presence of a sizable cooling process (Murase et al. 2014). More sophisticated simulations suggest even harder, bell-shaped escape spectra (Ohira et al. 2010;Globus et al. 2015a). Other acceleration mechanisms have been proposed that result in almost monochromatic particle spectra (Lyubarsky & Kirk 2001;Kirk & Giacinti 2017). Therefore, spectra of escaping charged particles that are significantly harder than E 2 are not unexpected for a single source. However, we note that in the current approach, we consider an entire ensemble of sources, and it seems unlikely that all sources will behave in the same way, i.e., reach the same maximal rigidity and have the same mass composition. Therefore, too hard or even peaked ( 0 g < ) spectra may be difficult to reconcile with current knowledge.
In the interest of comparability, we parameterize our generic source population exactly the same way as in the CF (Aab et al. 2017a;Alves Batista et al. 2019b), in which the nuclear species A (here 1 H, 4 He, 14 N, 28 Si, and 56 Fe) share a common spectral index γ and a maximal rigidity The A  are free normalization constants representing the number of particles ejected from the sources  Riehn et al. 2015;and QGSJETII-04, Ostapchenko 2011). The spread between the models (shaded areas) can be regarded as an interpretation uncertainty for the mass composition. per unit of time, comoving volume, and energy. The functional form of the cutoff is arbitrary, and we adopt the definition of the CF: In the CF, the fractions of injection elements f A are defined at a fixed energy point (10 9 GeV) relative to a total normalization. This definition is easily obtained from our A physically more meaningful definition of the mass fractions that does not depend on the arbitrary choice 10 9 GeV in Equation (4) is the integral fraction of the energy density where we choose E 10 GeV min 9 = as the lower boundary. We will mostly refer to I 9 A , providing the f A for comparability with the CF.
In Equation (4), the parameterization for the source evolution with redshift is given by the function For variable m, the function approximates all known continuous source density functions within the UHECR horizon z 1  . However, for the prediction of other messengers, it needs to be extrapolated to higher redshifts. In connection with the cosmogenic neutrino estimates (see Section 6), we will adopt more complex source distributions that include breaks. This flexible parameterization catches many features of theoretical source spectra. However, one has to keep in mind that the assumption of a rigidity-dependent escape is relatively strong and applies only to a subset of sources in which the maximal energy is limited by the size of the source rather than by cooling processes (Biehl et al. 2018a;Rodrigues et al. 2018). Another impacting assumption is that of a single dominant source population. The complexity can be increased by accounting for an additional proton component with higher rigidity (van Vliet et al. 2019) or even by detailed modeling of individual nearby sources (Eichmann et al. 2018). This, however, also vastly increases the degrees of freedom (dof), making a global fit of all free parameters unfeasible given the current statistics of the UHECR data.

Simulation Methods
In this section, we describe the methods of our global fit: the method used for the calculation of UHECR propagation through the IGM and the global fit of the propagated spectra to the observed data.

Propagation of UHECRs with PRINCE
To study the model dependencies in photonuclear cascades, we developed a new original computer code called PRINCE to efficiently solve the cosmic-ray transport problem. Instead of the Monte Carlo methods used in public codes, such as CRPROPA (Alves Batista et al. 2016) or SIMPROP (Aloisio et al. 2017), PRINCE numerically solves a system of coupled partial differential equations (PDEs) for the comoving density )of isotropically emitting and homogeneously distributed cosmic-ray sources. The terms (in order of occurrence) represent adiabatic cooling, pair production, photonuclear interactions (interaction and decays; reinjection) and injection from sources. The system of PDEs in E and z is solved using a sixth-order finite difference operator for the E derivatives and backward differentiation functions (BDFs), essentially an iterative implicit solver, for the redshift dependence. 3 The latter is required, since Equation (8) becomes stiff in z for nuclear systems (more details on the code and numerical methods are given in the Appendix).
Equation (8) is only valid under the assumption of a homogeneous source distribution with a separation much smaller than the diffusion length. For this case, the diffusion in extragalactic magnetic fields can be neglected. This reduces the calculation exclusively to the ballistic regime, in which the propagation becomes a 1D problem (time or redshift). This approximation in particular makes sense if one is interested in the highest energies above the ankle, where the impact of diffusion is small.
While similar codes have been previously developed, as, for example, in Allard et al.  (2015), our code stands out due to its very high computational speed and numerical precision. Even without significant architectural optimizations, PRINCE performs the computation of nuclear and neutrino spectra within 30 s on a single core, integrating an arbitrary injection spectrum that can contain elements with A 56  from a redshift of z=1. While Monte Carlo techniques for UHECR propagation become efficient due to the possibility of reweighting of precomputed events, our code shines when interest is devoted to model uncertainties, since we can essentially change any parameter and recompute within these 30 s, taking into account the impact on all relevant interaction rates. This includes arbitrary variations of the target photon densities without relying on simplified redshift scaling assumptions, as often employed in Monte Carlo methods or common numerical approaches. A detailed description of the numerical methods in PRINCE can be found in the Appendix.

Simulation and Fitting Procedure
This section aims to summarize the relevant setup of the simulations. We choose the five representative injection elements-hydrogen ( 1 H), helium ( 4 He), nitrogen ( 14 N), silicon ( 28 Si), and iron ( 56 Fe)-in accordance with the CF. We verified that choosing different injection elements of the same mass groups yields qualitatively similar results. The generic source model has eight free parameters: R max , γ, m, and free normalizations A  corresponding to the five injection elements.
We allow for a shift E d in energy within the systematic uncertainty given by Auger ( 14%  ; Fenu 2017). The transport equation (Equation (8)) is linear in the normalization factor A  but not in the other source parameters (γ, R max , and m), triggering us to employ a two-stage approach for the fit.
In the first stage, we discretize the parameter space for γ, R max , and m with the ranges and granularity given in Table 1.
For each point of this 3D source parameter grid, we separately compute the spectra at Earth for the five injection elements ( 1.5 10 6 · individual simulations for one choice of the photonuclear interaction model). Since the propagated spectra are linear in the A  values, the all-particle spectrum is calculated as a linear superposition of the results obtained for single-element injection.
In the second stage, we fit the nuclear fractions A  and energy shift E d to the spectrum and the first two moments of X max for each triplet in R m , , max g ( ) using the MINUIT package 4 (James & Roos 1975). The translation from individual mass spectra at the top of the atmosphere to X max á ñ and X max s ( ) is performed with the parameterization from Abreu et al. (2013), using updated parameter sets for SIBYLL 2.3 and EPOS-LHC.
To find the 2 c values for the UHECR fits within the entire 3D parameter space, the simulations are performed starting from redshift z 1 max = . Once the 3σ confidence intervals are localized, we run additional simulations starting from z=3 to compute cosmogenic neutrino fluxes, verifying that the previously derived contours are unaffected by higher redshifts. Both stages have to be repeated for each propagation model, while a change of the air-shower model only requires the repetition of the second stage. In all cases, the CIB model is fixed to Gilmore et al. (2012).
The following 2 c definition is used as the goodness-of-fit estimator: where 2  c refers to each of the three observables  , namely, the combined spectrum, X max á ñ, and X max s ( ). The total 2 c is obtained by summing. A nuisance parameter E d captures the uncertainty in the energy scale, and we assume its distribution to be flat within ±14%. The fit takes into account all data points above E 6 10 GeV min 9 = · . The global best-fit min 2 c is found by minimizing over all points of the 3D parameter space.
We then use 2 2 min 2 c c c D = to draw contours around the best-fit point by projecting to planes of two parameters by minimizing over all other parameters of the scan. While this frequentist approach is sufficient to draw contours and discuss the correlations among source parameters, there are more physical model parameters originating from the combination of discrete model choices, such as that for the photon background, disintegration, and hadronic interaction model. We did not attempt to parameterize these model uncertainties by continuous nuisance parameters, as these are impossible to define in a physically meaningful and unbiased sense. We therefore choose discrete model combinations and discuss their qualitative differences in the fit contours.

Impact of the Updated 2017 Data Set on the 2D Fit
We start the discussion of our results from the state of the CF and study the interesting impact of the updated 2017 data set (Bellido 2017;Fenu 2017) by reproducing a similar procedure to the one in Aab et al. (2017a) with our new code, PRINCE. The source evolution parameter is fixed to m=0 (flat evolution); the nuclear disintegration, CIB, and air-shower model are fixed to PSB (James & Roos 1975), Gilmore et al. (2012), and EPOS-LHC(Pierog et al. 2015), respectively. The minimization runs over the spectral index γ, R max , and the nuclear fractions A  . The energy scale is fixed and not allowed to float within its systematic uncertainty.
The energy range of the CF starts at 5 10 GeV 9 · . We noticed that with the new data set, the 2 c is significantly affected by the small discontinuity next to the X max á ñ point at 5.5 10 GeV 9 · ; i.e., this point alone adds a 35 2 c » to the best fit with a total 102 2 c » . We therefore treat this data point as an outlier and start our fit range at 6 10 GeV 9 · , which does not otherwise qualitatively impact the fit.
The contours are shown in Figure 2, and the best-fit values are summarized in Table 2. For the 2015 data set, we find the same qualitative result as the CF: a flat extended minimum with 1 g < and R 1 10 8 10 GV 9 max 9 < < · · and a second local minimum at 2 g » and R 4 10 GV max 10 » · . The differences in the exact locations of the minima with respect to the CF can be explained by the different propagation code used, as already pointed out in Aab et al. (2017a). Additional small shifts originate from the use of the experimental observables. While we fit the first two X max moments for the composition, the CF uses the full X max distribution. This has the strongest impact on the second minimum at 2 g = , which becomes less significant in our approach. In addition, we directly fit the combined unfolded spectrum and do not use a forward-folding procedure in the fit.
When switching to the 2017 data set, the best-fit parameters do not qualitatively change (see Table 2). However, the 2 c becomes worse due to the higher statistics. The allowed contours become narrower with a stronger preference for positive spectral indices. The second local minimum disappears. The reasons are the reduced statistical errors and a narrower width of the X max distribution at the highest energies of the 2017 data set, leaving less room for the combination of a high R max with somewhat softer spectral indices.
The largest qualitative difference concerns the injected iron fraction. While the 2015 data set did not require iron at the source, the new data suggest a small-but nonzero-integral iron fraction, I 2% Fe 9 » . This is also visible in the comparison of the best-fit spectra in Figure 3: for the 2017 data set (right panels), there is a contribution of heavy elements at the cutoff, which is absent in the fit to the 2017 data set (left panels). This is due to the higher statistics of the three highest-energy data points in the spectrum, which lead to a hardening. Due to the low rigidity found in the fit, reaching these energies requires a high charge number and therefore a significant iron fraction. However, this relies on the assumption of the rigidity dependence of the maximal energy and the fixed energy scale and hence cannot be rigorously interpreted as evidence for a nonzero iron fraction. Note, however, that it will still be visible if we later let the energy scale float. An indication of an iron contribution might also be visible in the composition data above 10 19.4 eV (Unger 2018).

3D Fit
We now include the source evolution m as an additional free parameter and allow the energy scale E d to float within the systematic uncertainties by following the procedure described in Section 3. First, we discuss our "baseline" case, defined by the combination of TALYS as the disintegration model and SIBYLL2.3 as the air-shower model (Section 5.1), before extending to other model combinations (Section 5.2). The impact of the model choices on the injected composition is discussed in Section 5.3.

Baseline Case Characteristics
Our "baseline" case is defined (a posteriori) by the combination of TALYS as the disintegration and SIBYLL2.3 as the air-shower model, motivated by its lowest 2 c out of realistic disintegration model choices. The other model combinations are discussed in Section 5.2.
The parameter space is shown in Figure 4, and the best-fit values are shown in Table 3. We note that the dof 2 c is close to 1, whereas it was close to 3 in the earlier 2D fit with a fixed energy scale and different disintegration and air-shower models; this means that we now actually have a good fit, due to the free source evolution and floating energy scale. The contour in the R max gplane is similar to that in the flat evolution case. Although the 1 g » corresponding to Fermi acceleration with diffusive escape is within the 95% contour, the preferred spectral indices result in flat or almost monochromatic spectra, 1 g < . In contrast to the previous 2D case, a floating E d allows for somewhat softer spectral indices. The R m maxplane exhibits a low-rigidity cutoff for every choice of the source evolution within the 95% CL. This is required by the composition data, in particular the X max s ( ), that suggest a clear separation among the mass spectra. This result can be interpreted as a signature of the preference of the data for the maximum-rigidity scenario with respect to the photodisintegration one. The discrimination among these scenarios is one of the science goals of AugerPrime (Aab et al. 2016b), and what we found constitutes a stronger result with respect to the 2D fit.
The m gparameter plane exhibits a clear anticorrelation, as already noticed, for example, in Unger et al. (2015) and Taylor et al. (2015). Positive source evolution (m 0 > ) results in a pileup from more distant sources, effectively softening the spectrum at Earth. This pileup is compensated by harder spectra at the source. Contrariwise, a high density of local   Note.For all quantities, the 1 σ uncertainties (for 1 dof) are given. , which is excluded from the fit. The expected composition is calculated assuming the EPOS-LHC shower model and comparing to the first two moments of X max distributions. is used to determine the contours, which are given for 1σ, 2σ, and 3σ (for 2 dof). In each 2D panel, the third parameter is treated as a nuisance parameter and minimized over to project the 3D parameter space. sources (m 0 < ) allows for spectral indices compatible with Fermi acceleration. The result clearly favors positive evolution, covering star-forming objects, GRBs, and blazars. The very hard spectra found in this case are consistent with what was found, for example, in Taylor et al. (2015). The 3s contours leave room for negatively evolving sources such as TDEs (Biehl et al. 2018b).
The spectrum and composition corresponding to the best fit of our baseline model are reported in Figure 5, while the corresponding injection spectra at the source (including the respective errors) are illustrated in Figure 6. The pileup effect from higher redshifts is clearly visible: while the injection spectrum is very hard ( 0.8 g = -), the propagated spectra are softer and have a stronger overlap. The best fit for the proton component is zero, and the proton component in the propagated spectrum comes only from propagation. However, the shaded regions in Figure 6 indicate the uncertainty in the normalization, which still allows for a significant proton fraction, as this component is barely contained in the fit range.

Model Dependence of the UHECR Fit
We expand the discussion of the previous sections and study the influence of the propagation and air-shower models by repeating the fit for permutations of the disintegration models PSB, TALYS, and PEANUT and the air-shower models EPOS-LHC, SIBYLL 2.3, and and QGSJETII-04. The results are shown in Figure 7 for the projection to the m gplane, and the corresponding best-fit parameters are reported in Table 4 (Appendix).
Consistent with what was found in the CF, we cannot find reasonable fits for QGSJETII-04 due to the model's broad X max distributions, in combination with a small X max á ñ, opposite to what is observed in data (Bellido 2017). In all other combinations, we find satisfactory best fits with dof 2 c ≈ 1.4-2.0. Clearly, the shower model has a stronger impact on the fit contours than the disintegration model, as can be seen by comparing the columns in Figure 7. Interestingly, for the PSB model in combination with SIBYLL 2.3, negative source evolution is excluded at 3s. This is an effect of the less efficient disintegration, as will be explained in the next section.
The anticorrelation between m and γ is found for all combinations of the disintegration and shower models (excluding QGSJETII-04). However, when exchanging SIBYLL 2.3 with EPOS-LHC, the 3σ contour in Figure 7 is shifted toward more local sources and/or more monochromatic spectra. The reason for this is that EPOS-LHC, compared to SIBYLL 2.3, predicts less shower-to-shower fluctuation, decreasing the X xmax s ( ), while at the same time, its X max á ñ predicts a lighter composition of the measurements. In combination, this allows for less overlap of individual mass spectra. Therefore, local sources are favored for this model,  Note.For all quantities, the 1σ uncertainties (for 1 dof) are given. = · , m=4.2). The shaded regions indicate the 1s uncertainties to the normalization of each injection corresponding to the fit (for R m , , max g fixed). While the bestfit proton fraction is zero, there can be a significant proton contribution within the uncertainty. reducing the impact of photodisintegration, which would increase the mass overlap. At the same time, the maximal rigidity R max is more constrained for EPOS-LHC than for SIBYLL 2.3, again decreasing the impact of photodisintegration (this is not directly evident from Figure 7).
The dof min 2 c is slightly worse when using EPOS-LHC ( 2.0 » ) compared to SIBYLL2.3 ( 1.4 » ), mainly because the fit to the X max á ñ is worse. It is, however, not strong enough to discriminate between these models, as the difference can be somewhat alleviated by allowing for shifts in X max within the systematic uncertainties. We did not include a proper treatment of these systematics.
Our results also show the limitations of what can be inferred from UHECR data alone. While the assumption of a generic rigidity-dependent source candidate describes the data sufficiently well, a strong degeneracy in the parameter space remains. Extending the range of the fit to lower energies could break this degeneracy but would require assumptions about the extragalactic magnetic field and the transition to a (possibly) Galactic component below the ankle, which means that it would add more dof to the model.
With new data from future experiments, the situation is expected to improve. For example, with better information on the UHECR composition from the AugerPrime upgrade, the parameter space will likely be more constrained. A significant improvement of the photodisintegration and air-shower models would be needed as well; otherwise, the ambiguity of the interpretation among different models will remain, as indicated by our results.

Injected Composition
An interesting and reoccurring question is the range of mass compositions permitted by Auger data. While the composition at observation is fixed (within the uncertainty of the air-shower models and data), it can have significantly different interpretations in terms of the composition ejected from the source. Within the limitations of our model, we illustrate the ranges of the injected fractions I A 9 within the 3s contours of our fit in Figure 8 as a function of the source evolution. The figure shows the baseline case TALYS-SIBYLL 2.3, as well as two additional panels changing the air-shower model to EPOS-LHC and the disintegration model to PSB, respectively.
Comparing the fraction ranges for SIBYLL 2.3 (Figure 8, left) with respect to EPOS-LHC (Figure 8, middle), the most striking difference is in the silicon fraction, which is significantly higher for SIBYLL 2.3, while in turn, the nitrogen fraction is higher for EPOS-LHC. This is mainly due to the heavier A ln á ñ predicted by SIBYLL2.3. A significant proton fraction is only found in the case of EPOS-LHC, owing to the slightly lower rigidity found for that model. In both cases, the nitrogen fraction increases at the cost of the helium fraction with higher source evolution. The higher disintegration for distant sources Table 4 Best-fit Parameters for the 3D Parameter Scan with Free Source Evolution for All Nine Model Combinations, as Described in Section 5.  produces more helium during propagation, therefore requiring less helium injected at the source. For the same source evolution, using SIBYLL2.3 with respect to EPOS-LHC leaves the mass fractions less constrained, as the combination of X max á ñ and X max s ( )predicted by SIBYLL2.3 allows for a stronger superposition of different mass spectra. In both cases, the allowed mass fractions widen when going to negative source evolution. This effect is directly connected to the propagation: for a larger concentration of distant sources, the disintegration increases the spread of masses limiting the initial spread, while a larger concentration of local sources allows for a broader spread of isotopes already at the source. This is an explicit demonstration that the X max s ( ) reflects not only the spread of nuclear masses at the sources but also what happens during their propagation to Earth (Abreu et al. 2013).
The impact of the disintegration model is qualitatively different. As mentioned in Section 5.2, negative source evolution is not contained in the 3s contours for the combination of PSB and SIBYLL2.3. This constrains the fraction ranges in Figure 8 (right panel) to positive source evolution. The most relevant features of the disintegration model are the level of α emission and the number of open reaction channels that control how efficiently a nuclear cascade develops. For instance, the absence of α emission in PSB is compensated by higher He fractions at the source, as noticed in Alves Batista et al. (2015) and Aab et al. (2017a). Due to the less efficient photodisintegration in PSB, the necessary development of the nuclear cascade can be ensured only if the sources are distant enough (positive evolution), leading to a rejection of local sources. This finding strengthens the need for using more refined models for photodisintegration, since it demonstrates that the simple PSB model might bias the predictions for source evolution while overestimating the amount of helium at the source. Figure 8, which describes the integral ejection fractions from the sources, can also be interpreted in terms of the physics of the sources. The helium and proton fractions are especially indicative of the amount of disintegration required within the sources. While the isotopes must escape rather intact from the sources for strong evolution, such as active galactic nuclei (AGNs), weaker source evolution seems to allow for higher helium and maybe even proton fractions, which implies that the nuclei may partially disintegrate in the sources. While this gives a rough estimate, a rigid interpretation requires a more sophisticated source model. For higher-luminosity sources that have a stronger disintegration chain, the rigidity dependence of the maximal energy is typically not a valid assumption; see, e.g., Biehl et al. (2018a) and Rodrigues et al. (2018).
A remarkable result is the nonzero iron fraction that we find throughout all model combinations. This is a result of the increased statistics at the cutoff of the updated Auger 2017 data set, as discussed in Section 4.

Cosmogenic Neutrino Fluxes
The source parameters inferred from the fit to UHECR data also lead to a prediction of the cosmogenic neutrino flux. However, cosmogenic neutrino fluxes are significantly affected by the cosmic-ray densities beyond a redshift of 1, while UHECR fluxes are almost insensitive to such distant source populations. Therefore, it is impossible to estimate any confidence interval using a solely data-driven method. Under the assumption that the fit is sensitive up to a redshift of z 1 max = , we draw in Figure 9 the neutrino ranges corresponding to the 1σ, 2σ, and 3σ contours of the fit with the baseline model combination. Essentially, these flux levels can  be regarded as constrained by the present data. In contrast to the 1σ region, which is limited to positive source evolution, the 3σ region is unconstrained toward negative redshifts (compare with Figure 4). Hence, if the sources are local, the expected cosmogenic fluxes are very low.
In the following, we exclusively focus on the 3σ contours. We study the robustness of our results against the changes of the disintegration and air-shower models in Figure 10. In the left (right) panel of Figure 10, the cosmogenic neutrino flux is shown corresponding to the blue UHECR contours for the models in the top row (left column) of Figure 7. The largest model dependence comes from the allowed range for the source evolution. The neutrino spectrum depends on the energy per nucleon; hence, the composition dependence is weak. The variations between the disintegration models are small, resulting in a relatively robust upper bound. For QGSJETII-04, the flux is small, since positive evolution is disfavored. For PSB, a sizable lower limit to the neutrino flux exists, since negative source evolution (local sources) is not allowed.
As the maximum rigidity is strongly constrained by the UHECR fit, the high-energy peak of the neutrino flux stays relatively robust and located at 10 GeV 8 . This is in agreement with Alves Batista et al. (2019b), where equally low fluxes were predicted. 5 A small but relevant difference resides in the propagation code, since Alves Batista et al. (2019b) assume a simplified redshift scaling of the CIB, whose effects in the neutrino fluxes are explained in Alves Batista et al. (2019a). If we apply the same simplified scaling, the cosmogenic neutrino flux in our calculations increases by 50%. Other minor differences come from other details of the propagation code and the fitting procedure. Differences from other works (Romero-Wolf & Ave 2018; Das et al. 2018;Møller et al. 2018) come from their limiting assumptions about the source evolution, injected composition, or cutoff energy.
The most significant impact on the fluxes comes from the extrapolation to redshifts z 1 > , which is unconstrained by UHECR data. For Figure 11, we adopt two approaches.
Left panel: empirical method using a simple continuation of the z 1 m + ( ) parameterization beyond z=1 up to z 3 max = .
We also test a distribution with a break at z=1 and a flat (m = 0) behavior beyond that. Right panel: discrete evolution functions of candidate source classes, where the parameter m is not free; AGN (Hasinger et al. 2005;Stanev 2008), GRB (Wanderman & Piran 2010), SFR (Yuksel et al. 2008; including starburst galaxies), TDE (Lunardini & Winter 2017), and a flat evolution. In this case, z 5 max = is used, which is above the cutoff for all source evolution used.
The most optimistic z 1 m + ( ) extrapolation results in fluxes that are one order of magnitude below the diffuse neutrino flux. It can be considered the upper limit of what is expected in the case of a single dominant UHECR source population with a rigidity-dependent energy cutoff. A flux at a similar level is found for AGN evolution. In either scenario, future radio-based instruments will not be able to distinguish between source types (right panel) or detect any significant cosmogenic neutrino signal. It is important to understand that the expected neutrino flux is (lower-) bounded only if the source evolution is fixed, motivated by a dominant source class. As long as the sources are not known or constrained, a "minimal cosmogenic neutrino flux" (Ahlers & Halzen 2012) is not meaningful.
The low neutrino fluxes are partly related to our choice of generic source model, which leads to fits with low maximal rigidity. Other scenarios are possible in which a small fraction of the UHECR flux originates from proton accelerators that reach GZK energies (van Vliet et al. 2019). These protons would copiously produce cosmogenic neutrinos off the denser CMB and peak at higher energies, while the majority of UHECRs would have a heavier mass composition, in line with current observations. These findings strongly support one of the science goals of the AugerPrime upgrade (Aab et al. 2016b), in which additional hardware is deployed to determine the proton fraction among the observed UHECRs. This should be regarded as being of utter importance for the decisions regarding the next-generation neutrino detectors. On the other hand, this result leaves room for an unambiguous detection of very high energy neutrinos from the sources directly, and it is unlikely that the cosmogenic flux will constitute a substantial background.

Summary and Conclusions
In this work, we have applied a new numerical highperformance propagation code, PRINCE, to the updated spectrum and the composition data published by the Pierre Auger Observatory in 2017. We have included the source evolution m as an additional free parameter. The savings in computation time have been used in favor of a detailed assessment of the main model dependencies, the nuclear disintegration in the propagation, and the hadronic interactions in the air-shower development. For the emission from generic UHECR sources, we have retained the main assumption from the CF of a single dominant accelerator type. Our results, therefore, refer to an "average" or "generic UHECR accelerator" that emits nuclei at most as heavy as iron with a spectral cutoff at a maximal rigidity.
We have demonstrated that the reduced statistical error of the 2017 data set, in particular at the highest-energy data points, favors for the first time a small but constrained iron fraction almost independent of the model variations. This implies a somewhat lower maximal rigidity.
The extension to three dimensions (γ, R max , and m) confirms and strengthens the finding of a low R max independent of the source evolution. We find a clear indication of a correlation between the spectral index and source evolution: rigiditydependent source candidates must be local, m 0 < , with spectral indices compatible with those obtained in models with diffusive shock acceleration or distributed according to the SFR but with very hard, almost monochromatic, spectral indices. Source classes discussed in the literature corresponding to such scenarios are jetted TDEs (Zhang et al. 2017;Biehl et al. 2018b;Guépin et al. 2018) and LL-GRBs (Zhang et al. 2018;Boncioli et al. 2019) or reacceleration scenarios, such as those proposed for termination shocks in starburst and nearby radio galaxies (Anchordoqui 2018;Eichmann et al. 2018;Winchen & Buitink 2018), respectively. While the inclusion of magnetic fields would soften the spectra at the source, the effect is probably not significant enough to draw an entirely different conclusion. It is challenging to reconcile this result with astrophysics, since a large number of alike sources with very similar R max and mass composition is needed to reproduce the observations.
We have assessed the impact of model variations on the contours in the γ-m plane for all combinations of the disintegration models PSB, PEANUT, and TALYS and the airshower models EPOS-LHC, SIBYLL 2.3, and QGSJETII-04. The largest effect comes from changes in the air-shower modeling, which means that a better understanding of hadronic interactions would provide useful constraints. However, the 3σ contours enclose the entire range of m, implying that there is no clear preference for a candidate source type. While the model variations lead to unconstrained distributions of the source, their mass composition is limited, preferring a mixture of nitrogen and helium with an admixture of silicon depending on the level and efficiency of nuclear disintegration during the transport. We have shown that the use of simplified disintegration models prevents the possibility of investigating the whole parameter space, including local sources. Other choices in the number or type of elements do not significantly affect the result.
By using the contours that represent compatibility with UHECR observations, we have studied the cosmogenic neutrino fluxes; compared to a purely theoretical prediction, this can be regarded as a postdiction from UHECR data. Because the allowed range in m is unbounded, no meaningful lower bound can be derived for cosmogenic neutrinos, since local sources cannot be excluded by the fit. On the other hand, we find that the upper bound is relatively robust under model variations. The fluxes are only constrained under fixed assumptions for the cosmic distribution of sources motivated by specific source classes.
In all cases, the expected flux is small and peaks at energies around 10 8 GeV, making the detection by the proposed future radio-based detectors unlikely. On the other hand, this result means that if very high energy neutrinos from sources exist at energies beyond 10 8 GeV, the expected background from diffuse cosmogenic neutrinos is expected to be small. This conclusion applies if UHECRs are produced in one dominant type of accelerator with rigidity-dependent maximal energy cutoffs. If there are multiple types-for instance, including a subset of proton-rich sources-then the fluxes can look Figure 11. Allowed range for the neutrino flux (all flavors) in the 3s region for different source evolution. Left: The purple range corresponds to z 1 max = (same as Figure 9). For the other curves, the source evolution is continued to z 3 max = either by continuing as z 1 m + ( ) (yellow) or with a break to flat evolution at z=1 (green). Right: The ranges are shown for the source evolution fixed to different source classes and for flat evolution. significantly different. Additional clues from high-precision composition measurements are highly valuable, which the AugerPrime upgrade is expected to deliver a few years from now.
In conclusion, the fit is relatively sensitive to the disintegration model and, even more, the air-shower model, which still leads to a strong ambiguity in the interpretation of the data and therefore needs future improvements. The predicted cosmogenic neutrino flux is relatively robust with respect to these models and probably out of the reach of future experiments in all cases. A significant enhancement to the neutrino flux can come from redshifts beyond 1, which cannot be constrained from UHECR data alone.
We thank A. van Vliet for useful feedback on the draft of this paper and T. Piran for inspiring discussions. We also thank the colleagues from the Pierre Auger Collaboration. This work has been supported by the European Research Council (ERC) under the European Unionʼs Horizon 2020 research and innovation program (grant No. 646623).

Appendix Propagation Code: PRINCE
For our study, we have written an original computer code in order to have a framework in which systematic uncertainties such as cross sections and photon backgrounds can be efficiently varied. This appendix contains details about the numerical methods used to accelerate the computation of the UHECR transport equation.
The two popular public UHECR propagation codes (CRPROPA, Alves Batista et al. 2016;SIMPROP,Aloisio et al. 2017) use a Monte Carlo approach. While these can effectively handle spectral properties by reweighting samples, a rigorous treatment of certain systematics, such as photonuclear cross sections, requires a full, computationally expensive resampling. On the other hand, an iterative numerical solution of the transport equation system requires a constant computational time under the variation of any parameter. The trade-off is that the variation of spectral properties requires a full recomputation as well.
Our code is called PRINCE. The main development goals were as follows.

A time-dependent UHECR transport equation solver
efficient enough to compute a single spectrum within seconds. 2. Fast and easy variation of model input, such as crosssection models and extragalactic photon backgrounds. 3. Accessibility and modularity, such that users can easily modify and extend specific parts of the code through interfaces.
To achieve these goals, PRINCE is written in pure PYTHON using vectorized expressions for the performance-intensive parts, accelerating those using libraries like NUMPY and SCIPY (Jones et al. 2001). This vectorized approach also allows for the code to be implemented for massively parallel accelerators, such as graphics processing units, without much additional effort.
The Boltzmann transport equation for UHECRs is most conveniently solved in terms of the comoving density ). Assuming homogeneous and isotropic sources, the diffusion terms vanish, and the transport equation becomes independent of the spatial coordinate x (propagation theorem; Aloisio & Berezinsky 2004). The coupled differential equation system for the particle species i reads where we introduced the simplified notation Y Y E z , ), which can be transformed between time t and redshift z with the relation dz dt z H z 1 = -+ ( ) ( ). The first two terms describe the continuous energy losses due to adiabatic cooling (HE) and Bethe-Heitler pair production (b e e + -). Here i G is the rate of the photonuclear interactions. The conversion of the particle species j into i is handled by the reinjection terms Q Y j i j  ( ). The decay terms for unstable particles can be treated implicitly, as described below. The last term (J i ) describes the injection from sources. We will discuss the partial and ordinary differential parts separately in the following two sections.

A.1. Photohadronic Interactions: Ordinary Differential Equation
Our approach to solving the ordinary differential equation (ODE) system that describes the conversion between particle species due to photonuclear interactions follows the method and notation described in Boncioli et al. (2017) and Biehl et al. (2018a). This new approach, however, greatly benefits from rewriting the same equations in terms of matrices.
Photonuclear interactions above a few MeV disintegrate the projectile nucleus, resulting in the production of multiple final state particles. In the system of the ODE, the disintegration happens with the rate E i i i G º G( ), and the (re)injection terms Q Y E , j i j i  ( ) couple the equation systems of different particle species. The general form of the interaction rate on a target photon field is given by an integral over the photon energy ε and pitch angle θ in the comoving frame ) and an additional integral over projectile densities Y E j ( ). The inclusive differential cross section can again be pitch angle-averaged and expressed as a function of y.
In analogy to Equation (12) The decay of unstable particles is governed by a term Y is the lifetime of an unstable particle or nucleus i at rest. The reinjection terms for the decay products have a similar form to Equation (13) but do not depend on the photon field. Hence, the second integral can be omitted: , . 15 The redistribution function dn dE j i i  is, in this case, the inclusive energy distributions of the decay product i in decays of j. To obtain inclusive distributions, all decay channels that contain i are summed with their branching ratio as weight.
Most unstable particles that occur in UHECR propagation have a mean lifetime much smaller than the other relevant timescales. Hence, the decay can be regarded as an instant process at the production vertex.
. For UHECR propagation, we set thresh t to ¥; i.e., all unstable particles decay immediately. A special case arises for secondary nuclei. At high energies (E i  TeV), the impact of the internal nucleon motion can be neglected to a good approximation, resulting in the conservation of the boost of secondary fragments; i.e., the energy per nucleon is conserved. The redistribution function then simplifies to For the discretization (see next section), it is convenient to formulate the equation system in E i A . This makes the treatment of the δ-function in Equation (18) accurate as long as the same grid in E i A is chosen for all nuclear particle species. We use the form of Equation (19) for all nuclear species in the code. However, for the sake of brevity, we will not mention this explicitly in the following and only discuss the more general form of Equation (13).

A.2. Discretization
For the numerical solution of the coupled ODE system (Equation (10)), we introduce a discrete, logarithmic grid in energy, where the grid constant d k can be adjusted independently for the particle and the photon grids to achieve the desired precision. Currently, eight points per energy decade results in a good compromise between precision and computational speed. We use k l m , , as upper indices for the energy grid indices and i j , as lower indices for the particle species. All quantities are represented by their value at the interval centers. In some cases, such as for strongly peaked cross sections, it is necessary to compute precise averages over each interval instead of taking the central value.
On a grid, we rewrite the interaction rate from Equation ( Since each projectile produces only a few secondary particle species, the matrix Φ is sparse, with only 2% » of nonzero elements. The ordering of Y by energy and particle mass results in an upper-triangular shape of Φ and its submatrices, as long as there is no particle acceleration. The calculation of the derivative, a sparse matrix vector dot product, is significantly accelerated by using a sparse matrix storage format from a specialized library. The compact sparse row (CSR) format 6 stores a matrix M as three vectors: a data vector D containing only the nonzero elements, a column index vector C holding the column indices for each element, and a row pointer R pointing to the position of the first element of each row in D and C. The end of each row is given by the next index in R, and an empty row is indicated by a repeated index in R.
For example, the matrix The format is to be read as follows. The first two entries in D and C belong to the first row of M, as R 2 1 = signals that the second row starts with the third entry. With C giving the column position, this means that M 6 00 = and M 1 03 = . A repeated entry in R indicates an empty row, as for R R 3 2 3 = = in the example. The vectors D and C therefore always have a length equal to the number of nonzero elements, while R has a length equal to the number of rows plus one. The compact sparse column (CSC) format is defined analogously. The CSR format is especially effective for multiplication with column vectors.
In our approach, the particle production channels and therefore the nonzero elements of Φ in Equation (25) are fixed. Therefore, the column index vector and row pointer only have to be found once. Instead of recomputing the whole sparsity structure, only the elements of the data vector in the sparse matrix format of Φ have to be replaced in every step, resulting in further computational speed gains.
The computation of elements of Φ can be done in a single matrix expression if  and  are combined into a single crosssection kernel . By ordering  according to the order of the D vector of Φ, the elements of D can be modified in place without additional memory allocations: This arrangement allows for very fast computation of all coefficients of Φ; hence, the handling of the time/redshiftdependent ODE system becomes very efficient. The cross sections can be varied by scaling or replacing elements of the kernels in  between runs without additional initialization overhead.

A.3. Adiabatic Expansion and Pair Production: PDE
The partial differential part of the transport equation comes with two continuous loss terms: with the loss terms b dE dt º . The adiabatic losses due to cosmological expansion are described by more stable for our purpose to use forward-biased differences, e.g., in second order, The code allows us to adjust the order of the finite differences to optimize for the given problem. Currently, we use sixthorder finite differences. While this is probably more than necessary, we find that the impact on performance is small, as the computation time is dominated by the photohadronic part. For applications different from UHECR propagation, we might, however, have to revisit this choice. If the order of the operator does not change, D i kl can be included in the sparse interaction matrix Φ from Equation (23) that is solved as an ODE with the methods described in the next section.

A.4. Differential Equation Solver
Using matrix formulation, we have found an efficient scheme to recalculate the time derivative Y z t ¶ ( ). To solve the problem for Y z ( ), one has to choose an integration scheme in time t (or, for redshift z, by converting with dz = dt z H z 1 -+ ( ) ( )). For a system with light injection, the eigenvalues of the interaction matrix Φ are small enough such that we can use an explicit Euler scheme: For a proton system from redshift z=1 with dz 10 3 = -, the propagation can be solved within a few hundred ms.
For heavier-mass nuclei, the eigenvalues of Φ become very large, and the system becomes stiff, requiring very small time/ redshift steps for a stable explicit integration. In this case, we use an implicit integration scheme based on the SCIPY. INTEGRATE.ODE.BDF solver, which adaptively adjusts the step width and the order. A first-order BDF scheme corresponds to an implicit Euler scheme: .