Cosmic-Ray and Neutrino Emission from Gamma-Ray Bursts with a Nuclear Cascade

We discuss neutrino and cosmic-ray emission from Gamma-Ray Bursts (GRBs) with the injection of nuclei, where we take into account that a nuclear cascade from photo-disintegration can fully develop in the source. One of our main objectives is to test if recent results from the IceCube and the Pierre Auger Observatory can be accommodated with the paradigm that GRBs are the sources of Ultra-High Energy Cosmic Rays (UHECRs). While our key results are obtained using an internal shock model, we discuss how the secondary emission from a GRB shell can be interpreted in terms of other astrophysical models. It is demonstrated that the expected neutrino flux from GRBs weakly depends on the injection composition, which implies that prompt neutrinos from GRBs can efficiently test the GRB-UHECR paradigm even if the UHECRs are nuclei. We show that the UHECR spectrum and composition, as measured by the Pierre Auger Observatory, can be self-consistently reproduced in a combined source-propagation model. In an attempt to describe the energy range including the ankle, we find tension with the IceCube bounds from the GRB stacking analyses. In an alternative scenario, where only the UHECRs beyond the ankle originate from GRBs, the requirement for a joint description of cosmic-ray and neutrino observations favors lower luminosities, which does not correspond to the typical expectation from {\gamma}-ray observations.


Introduction
The origin of ultra-high-energy cosmic rays (UHECRs) at the highest energies 10 18 eV, which are most likely of extragalactic origin, is one of the unsolved mysteries in astroparticle physics. Recent results of the Auger Collaboration favor a mixed chemical composition of cosmic ray particles injected from the sources into the extragalactic medium (Aab et al. 2017). One possible candidate source class are gamma-ray bursts (GRBs; Piran 2004), which are the most energetic electromagnetic outburst class, and which we consider in this study.
GRBs are expected to produce a significant flux of highenergy neutrinos (Waxman & Bahcall 1997) if cosmic ray protons interact with the target photons produced in the prompt gamma-ray emission expected to be created by internal shocks. GRBs, however, are the best-tested source class in stacking searches of the IceCube neutrino telescope using the gammaray counterpart, because these analyses make use of timing and direction to almost completely suppress the atmospheric background. Therefore, neutrino telescopes can efficiently test the GRB-UHECR paradigm (Abbasi et al. 2012). Although earlier prompt neutrino flux predictions (Guetta et al. 2004) have been updated in the meantime Li 2012;He et al. 2012), it is clear that current IceCube data (Aartsen et al. 2017) exert pressure on the allowed parameter space for conventional neutrino emission models and the GRB-UHECR paradigm . Regions in parameter space for which the pion production efficiency is high and neutrons are efficiently produced together with the neutrinos, which can easily escape from the source (Ahlers et al. 2011), are already excluded. If, however, the cosmic ray escape is dominated by directly escaping protons (for which the Larmor radius can reach the size of the shell), fewer neutrinos are expected, and the charged cosmic rays are ejected with a hard spectrum ∝E −1 , whereas neutrons, which are not magnetically confined, escape with a softer spectrum (Baerwald et al. 2013; also discussed in Globus et al. 2015a;Unger et al. 2015) from the source perspective. Such hard spectra are indeed found as escape spectra from the sources in current UHECR fits (Aab et al. 2017).
There are a few important caveats in this picture. First of all, other prompt emission models have been considered (see, e.g., Murase 2008;Zhang & Kumar 2013) and tested (Aartsen et al. 2017); for example, the magnetic reconnection model (Zhang & Yan 2011) is less constrained than the internal shock or photospheric scenarios. Second, all of the mentioned constraints have been obtained in one-zone collision models, that is, the parameters have been assumed to be the same for all collisions. Multi-zone models, however, tend to predict lower neutrino fluxes Globus et al. 2015a;Bustamante et al. 2017). The reason is that the pion production efficiency strongly drops with collision radius, which means that few collisions at inner radii dominate the neutrino flux, whereas cosmic rays and gamma-rays (on average) come from larger collision radii . And third caveat is that all of the IceCube constraints have been derived for proton primaries. To test whether they also apply to a heavier composition in the jets, is one of the motivations of this work.
The behavior of UHECR nuclei in GRB jets has been studied in Anchordoqui et al. (2008); Wang et al. (2008); Murase et al. (2008); Globus et al. (2015a); Joshi et al. (2016) and Boncioli et al. (2017), where in many cases the motivation was to study the necessary conditions for UHECR survival. If, however, the radiation densities are high, a cascade of isotopes lighter than the injected primary emerges due to photodisintegration in the source (Globus et al. 2015a;Boncioli et al. 2017). In this work, we refer to this phenomenon as "the nuclear cascade". In Boncioli et al. (2017) the nuclear cascade in a GRB shell has been self-consistently computed with a disintegration model at a level of sophistication comparable to CRPropa used for the cosmic ray propagation (Kampert et al. 2013), and it has been demonstrated that the feed-down to lower masses can have an impact on the ejected cosmic ray (and also the neutrino) flux. It can also be affected by the photodisintegration modelespecially if a more sophisticated TALYS (Koning et al. 2007) or FLUKA-based (Ferrari et al. 2005;Böhlen et al. 2014) disintegration model is used as compared to the Puget-Stecker-Bredekamp disintegration chain (Puget et al. 1976). The same approach will be used in this study, applying it for the first time to systematic parameter scans over the GRB parameters, taking into account the corresponding neutrino constraints.
In order to test the GRB-UHECR paradigm, a combined source-propagation model is needed: the accelerated nuclei are injected into the radiation zone, where the secondaries are produced, escape from that zone, and are then propagated through the extragalactic space to Earth (see Baerwald et al. 2015;Globus et al. 2015a;Unger et al. 2015, for examples). One of the open questions from the UHECR perspective (apart from the impact of the composition) is where the transition energy between the lower-energy (possibly Galactic) and UHECR contribution occurs. Especially for protons, there are already constraints from cosmogenic neutrino data (see, for example, Heinze et al. 2016) and extragalactic diffuse gamma-ray data (Supanitsky 2016;Berezinsky et al. 2016) if the transition occurs below the ankle; this is the so-called "dip model" (Aloisio et al. 2007). In combined source-propagation models, a transition at the ankle has been considered in Globus et al. (2015aGlobus et al. ( ,b, 2017 for UHECR nuclei. A GRB proton dip model clearly overproduces the prompt neutrinos , while more generic models can effectively describe the transition to a lighter composition below the ankle by disintegrated nucleons (Unger et al. 2015), which, however, does not consider the source-class dependent pion production efficiency in the GRBs.
We test both hypotheses in this study: GRBs as the sources of UHECR nuclei both covering the ankle ("Mixed Composition Dip Model") and above the ankle ("Mixed Composition Ankle Model"). We inject a pure composition of nuclei into the jet in order to study the impact of the injection composition on the UHECR and neutrino fluxes. While it is clear that a mixed injection composition will improve the UHECR fit, using a pure composition allows us to develop a deeper understanding of the behavior of UHECR nuclei in GRBs, and the impact of the injection composition on the predicted prompt and cosmogenic neutrino fluxes. We focus on the one-zone emission model, which is considered in the prompt neutrino analyses, extended to heavier nuclei, and we perform systematic parameter space studies. While most of our results are derived for the internal shock model, we demonstrate how our results can be translated into other emissions scenarios; therefore the key parameters used in this study are emission radius R, gamma-ray luminosity L γ , baryonic loading ξ i , and injection isotope.

Methods
Our simulations are based on the NeuCosmA code, following the implementation for GRBs in Hümmer et al. (2012); Baerwald et al. (2012Baerwald et al. ( , 2013Baerwald et al. ( , 2015 for protons. We focus here on the extensions for nuclei; this description complements the short summary given in Boncioli et al. (2017) for a specific case.

Treatment of the nuclear cascade
Our method is fully deterministic, based upon solving a coupled system of partial differential equations (PDEs), which, for particle species i (such as a nuclear isotope), reads: where b ′ (E) = E ′ t ′−1 loss (with the energy loss rate t ′−1 loss ), and t ′−1 esc is the escape rate. The PDE system is to be solved for the differential particle densities N ′ i [GeV −1 cm −3 ] in the shock rest frame (SRF); primed quantities refer to the SRF. The coupled system arises because of the injection term which allows for injection from an acceleration zone Q ′ i (typically a power law), as well as for injection from other species j with the term Q ′ j→i , such as from photodisintegration or β ± decays. We will discuss the relevant interaction processes below. Fully-stripped ions are identified by A Z, where A is the mass number, Z is the charge (or number of protons), and N = A − Z is the number of neutrons 1 .
A schematic illustration can be found in Fig. 1, where each dashed (gray) box corresponds to one PDE. Several hundred nuclear isotopes will enter later, which means that automatic techniques for the PDE system setup will be used, and pictorial representations are useful to describe the models. In the pictorial representation, the cooling, escape, and injection from the radiation zone Q ′ i act only on one species i, and are therefore "attached" to that, meaning that the PDE system is set up with corresponding terms. The injection term Q ′ j→i acts on species i (and is attached to that PDE), but involves another species' density N ′ j , illustrated by the arrows from a different dashed box corresponding to different species j; it couples the PDEs.
In example A) the computations in Hümmer et al. (2012); Baerwald et al. (2012Baerwald et al. ( , 2013Baerwald et al. ( , 2015, which are effectively performed in the optically thin regime to photohadronic (pγ, nγ) interactions, are extended to optically thick sources. In this case, there are only two PDEs for the nuclear species (protons and neutrons). Here, only protons are injected from an acceleration zone (the shock). Since the PDEs are coupled by the photohadronic interactions, protons will be converted into neutrons, and vice versa. In the optically thin regime, the interactions pγ hardly take place, which means that the neutron species will not be populated. Typically, cooling processes for protons are synchrotron cooling and adiabatic cooling, whereas the escape term is driven by the photohadronic interaction rate. Since in fact a fraction of the interacting protons produces protons again, we treat the photohadronic energy losses discretely, that is, we reinject these protons at lower energy. Neutrons, on the other hand, can typically escape from the region because they are not magnetically confined, or by photohadronic interactions 2 . Although this simple example is only meant for illustration, we note that center-of-mass energy (if applicable). In fact, for photo-meson production, one has to integrate over dn j→i /dx within the photon energy integral, which is "hidden" in the interaction rate here; see Appendix A. The function describes the distribution of secondaries of type i per final state energy interval dx, the probability distribution is normalized p j→i (x)dx = 1, and the integral over that function gives the average number of secondaries M j→i produced per interaction (multiplicity). The specific implementation of this injection function and the computation of the interaction rate depends on the nuclear process considered (for details see Appendix A).
In the main text, we focus on a qualitative description of the interactions.

Beta decays and spontaneous emission
As a starting point, we choose all potentially relevant known isotopes with A < 56 and 56 Fe, which is our heaviest (stable) injection isotope, taken from the database in Mathematica (2016). These isotopes are shown as boxes in Fig. 2 as a function of N and Z, where the considered injection isotopes (most abundant stable elements) are shown in dark blue. Apart from the blue isotopes, which are stable (or have lifetimes longer than one month), all unstable isotopes undergo various decay processes: β decays (red), possibly followed by the delayed (later than about 10 −14 s) spontaneous emission of one or two nucleons or an α particle (light red), or spontaneous emission of one or two nucleons or an α particle (white). Our framework is already simplified by eliminating branchings <5%. In the figure, we only show the leading branchings, as many isotopes have several decay channels that we take into account in the computations. Typical lifetimes of β ± emitters range from fractions of seconds to hours, and have to be Lorentz-boosted to be compared to the dynamical timescale t ′ dyn of the prompt emission in the SRF (which will be of the order of one second). Beta emitters with rest frame lifetimes τ 0 10 −5 s are potentially interesting for the secondary neutrino production in GRBs, which are the red boxes marked with "×" in the figure. These beta emitters decay faster than the dynamical timescale for γ ′ A ≡ E ′ A /m A 10 5 , which corresponds to neutrino energies 10 4 GeV where the neutrino flux from beta decays may dominate. One can, however, see from the figure that these interesting isotopes are relatively far off the main diagonal, which will be hardly populated. Neutrinos from β ± decays therefore play a minor role from within the GRBs (but are included in the computation), while escaping neutrons and unstable isotopes may decay on the way to Earth. The neutrinos from neutron (and unstable nuclei) decays will only show up at very low energies for GRBs (see, e.g., Baerwald et al. 2011Baerwald et al. , 2012, where the component is explicitly shown). We include these neutrinos as part of our prompt neutrino flux. For this component, we let the neutrons and unstable isotopes decay outside the source, since it can be shown that the neutrons and unstable isotopes will rather decay than interact at these low energies where the beta decay component dominates the neutrino flux. However, we do not show this component explicitly since its contribution will be over-estimated for high energies, where the nuclei may interact before they decay.
The other interesting class of emitters are proton-or neutron-rich spontaneous emitters (white boxes in Fig. 2), which typically decay very quickly by the ejection of nucleons (without accompanying neutrinos). These spontaneous emitters can be integrated out (replaced by their daughters) if the lifetime is short enough, that is, they decay at the highest energies with rates faster than that of all the other radiation processes. We estimate that isotopes with τ 0 10 −10 s can be integrated out (white boxes with dots), because then γ ′ A τ 0 ≪ t ′ dyn at the highest energies (where γ ′ A = E ′ A /m A ∼ 10 10 ). Isotopes with unknown lifetimes are also assumed to decay faster than this threshold. Since isotopes far off the main diagonal will be hardly populated, this mainly affects the light unstable isotopes with A 6.

Photodisintegration and photo-meson production
The disintegration of nuclei in interactions with photons can be qualitatively distinguished by the energy scale ǫ r (photon energy in the nucleus' rest frame). For 8 MeV ǫ r 150 MeV, the "giant dipole resonance" (GDR; Goldhaber & Teller 1948) and other processes lead to an electromagnetic excitation of the primary nucleus with the emission of one or more nucleons (or light nuclei). We refer to this regime below the pion production threshold as "photodisintegration". For ǫ r 150 MeV, higher energy processes, such as baryonic resonances, dominate the disintegration, accompanied by meson production; we refer to this regime as "photo-meson production" regime.
Due to the power-law type of the projectile and photon spectra, the nuclear cascade will, to leading order, be determined by photodisintegration above the threshold for GDR, which requires that target photons with the right energies are available as interaction partners. The photodisintegration within the sources has been discussed in Boncioli et al. (2017) from the perspective of the nuclear interaction model, where the impact of different model choices is demonstrated. For this work, our photodisintegration model is based on cross-section information from TALYS 1.8 for nuclei with A ≥ 12 (Koning et al. 2007) and for lighter nuclei on tabulated data, which is distributed with CRPropa2 (Kampert et al. 2013;see Boncioli et al. 2017, for details).
Our photo-meson interaction model is a superposition model based on the SOPHIA Monte Carlo generator (Mucke et al. 2000). Superposition means that the nuclei are treated as a bulk of independent nucleons with E ′ p/n = E ′ A /A. In each interaction, only one nucleon interacts and is ejected from the nucleus. The energy distributions of secondary pions and nucleons are computed with SOPHIA, whereas for the residual nucleus the energy per nucleon is conserved The probability that a neutron within the nucleus interacts is N/A, and Z/A if a proton interacts, respectively. The inelastic cross section is assumed to scale like σ Aγ = Aσ pγ . Although this model is similar to that being used currently in the literature (see, e.g., Anchordoqui et al. 2008;Murase et al. 2008;Kampert et al. 2013), it has clear deficiencies. First of all, the scaling of the cross section may not be correct; for instance, Schlickeiser (Schlickeiser 2002) directly proposes the "Glauber rule" A 2/3 . It is, however, well established that at energies above 150 MeV and several GeV in ǫ r , the cross section scales proportionally to A for nuclei heavier than 4 He (MacCormick et al. 1997). At higher energies the photon interacts in a more hadron-like way, where the A 2/3 behavior is more appropriate. We motivate our current assumption by the choice of the astrophysical source class, where lower interaction energies are more relevant. Second, the probability for the ejection of a proton or a neutron is close to 50−50, which drives the residual nucleus further away from the main diagonal than one would expect from a spectator model (Rachen 1996).
In the main text of this work, we assume that the minimal photon energy is low enough for the high-energy cosmic ray nuclei to find interaction partners above the GDR threshold. This means that the photo-meson production will be of secondary importance. In the case where the photon spectrum is cut off at low energy, which can occur due to synchrotron self-absorption (Wang et al. 2008;Murase et al. 2008), the disintegration at the highest energies will be dominated by the photomeson regime. We show an example demonstrating the quantitative impact in Appendix B (which can be followed after reading Sect. 3), where it is demonstrated that the neutrino flux is hardly affected, whereas the UHECRs can actually extend to higher energies in that case. We will demonstrate in this work that the neutrino production in GRBs, depending on luminosity and collision radius, is typically dominated either by the photo-meson production off the injection nucleus, by the photo-meson production off the secondary nuclei produced in the nuclear cascade, or by protons and neutrons produced in the nuclear cascade. The impact of the photo-meson model for these two cases is quantitatively different: while this photo-meson model produces reliable results for proton and neutron projectiles, we have to interpret the cases where the primary nucleus dominates with care due to the simplifications listed above. A more dedicated study of photomeson production off nuclei will therefore be needed in the future.
The implementation of the photo-meson model follows a logic similar to Hümmer et al. (2010), where the idea was to discretize one of the integrals in the double integration, needed to compute the secondary spectra, into so-called interaction types to gain linear (compared to quadratic) computational complexity. Our novel approach, described in Appendix A, uses a new discretization scheme that extends the previous scheme to cope with interactions of various isotopes in the future. We apply this scheme to SOPHIA in the Appendix, and demonstrate that it can outperform Hümmer et al. (2010) in terms of efficiency and precision.

Other radiation processes
In addition to decays, photodisintegration, and photo-meson production a number of radiation processes are implemented as continuous energy loss (cooling) processes. At ǫ r 8 MeV, the Aγ interactions are determined by the quantum electrodynamics scale. As a consequence, charged isotopes will produce, via the Bethe-Heitler process, e + e − pairs, which we include following Blumenthal (1970). However, for GRBs, this contribution is typically sub-dominant; we therefore do not show it explicitly in our timescale plots.
All charged species are assumed to cool via synchrotron radiation and adiabatic cooling (roughly determined by t ′ ad ≃ t ′ dyn ). Cosmic ray nuclei are assumed to be magnetically confined, but can directly escape within their Larmor radius from the shell boundaries (see discussion of cosmic ray escape below). This escape contribution will always have at least a sub-dominant impact on the source density. We add, however, an escape term for the neutrons with the free-streaming timescale. This term is consistent with our assumptions for neutron escape and, in fact, necessary to reach the steady state within a few times the dynamical timescale. Since neutrons can only decay at the lowest energies, the steady state can be reached there at around the neutron rest frame lifetime in the absence of this term.

Automated isotope selection scheme
One of the key features of our new code is the automatic isotope selection scheme: it picks the isotopes from the 481 isotopes in Fig. 2 that are relevant for the problem under investigation, depending on the requested precision. We first of all load the isotope data corresponding to Fig. 2, and flag all isotopes with lifetimes τ 0 ≤ 10 −10 s to be integrated out. As a second step, we load the photodisintegration and photo-meson production models. If daughter nuclei are encountered that are to be integrated out, all paths are recursively followed to the next "active" isotopes to be treated. In this case, we replace the interaction by an effective one, directly producing the next active isotopes. Consider, for example, the disintegration of a nucleus a with consecutive fast decays a → i 1 → . . . → i N → b, where the isotopes i 1 . . . i N are to be integrated out because they interact quickly. One can show that the secondary injection of b can be written in the form of Eq. (A.7) where effective χ a→b ≡ are to be used. Thus, when one constructs the isotope chart, one can replace the daughter of that process recursively by all possible end states of fast decay chains with effective values of χ and g. The implemented techniques are especially powerful and sophisticated for photomeson production, where the interaction types (x-values, see Appendix A) have to be re-defined automatically. Figure 2 illustrates that the lifetimes of the β ± emitters are longer than the chosen lifetime threshold for integrating out isotopes for GRBs (i.e., there are no dot-marked red boxes), which means that neutrino fluxes from integrated out β ± emitters do not need to be taken into account.
As a next step, we identify the relevant channels and set up the PDE system Eq. (1) automatically. Starting from an injection isotope (such as 56 Fe), we recursively follow all photodisintegration and decay channels (including mixed paths) for different (close enough) center-of-mass energies, and flag the encountered isotopes as relevant for the computation. We impose a threshold on the secondary multiplicity for these recursive paths, as daughters with very low multiplicities will be hardly filled. Our software supports several different selection methods and activation schemes, but the results depend on the chosen threshold values rather than the scheme. For a "two-dimensional" disintegration model (in N and Z), such as TALYS, one should not pick a threshold that is too large, as many isotopes with low multiplicities will be filled off the main diagonal, which means that the nuclear cascade would otherwise cease. Figure 3 illustrates the selected isotopes for the injection of 56 Fe for different threshold multiplicities. The thresholds M > 0.2 and M > 0.05 correspond to the most important isotope selections "Astro, priority 1" and "Astro, priority 2", respectively, in Boncioli et al. (2017; which were, however, computed for all possible injection elements there). In this work, we choose conservatively M > 0.01, which leads to 233 selected isotopes with about 5000 attached photodisintegration channels (from originally 41 000 inclusive channels extracted from TALYS), ten decay channels, and about 1300 photo-meson channels for the GRBs. We have checked that the results hardly change for smaller values of M. In fact, for the code optimization, one would first of all choose a small threshold and then increase it as long as the results are not affected. Of course, for different injection isotopes, Fig. 3 will look different; it is therefore our procedure that defines the isotope chart, which can also be applied to different source classes.

Energetics of the source
The photon spectrum in the source is assumed to follow observations. The spectrum for long-duration GRBs can be typically described by a broken power law where C ′ γ is a normalization factor and ε ′ γ,br ≃ 1 keV (fixed in this work). Our minimal and maximal photon energies are chosen as small and large enough, respectively, such that there are no longer any significant effects on the neutrino flux from the photon spectrum. This especially implies that nuclei will always find interaction partners for disintegration at the GDR (see Appendix B for the effect of the minimal photon energy cutoff).
In order to compute the photon density in the SRF, we define an "isotropic volume" of the interaction region with shell width ∆d ′ and the radius (distance from emitter) of the emission region R. Because of the intermittent nature of GRBs, the total fluence is typically assumed to be coming from ∆T/t v such interaction regions, where ∆T is the duration of the prompt emission. We model the emission from one such region first, as this will allow for simpler interpretations in terms of multi-zone models such as Bustamante et al. (2015Bustamante et al. ( , 2017. The normalization of the photon spectrum in Eq. (5) is obtained from where L γ is the isotropic equivalent luminosity in gamma-rays and Γ is the Lorentz boost factor. The integration limits are taken from the Fermi-GBM range from 8 keV to 40 MeV in the observer's frame. The factor Γ 2 comes from boosting the luminosity from the source (engine) frame to the SRF, and does not come from geometry estimators, that is, Eq. (7) shows how the target photon density scales with Γ and R without any geometry relationship for the radius. It scales with R −2 , which comes directly from the volume of the production region in Eq. (6) and therefore strongly depends on the radius. It is easy to show that the same dependence on R enters the pion production efficiency (fraction of proton energy dumped into pion production in the optically thin case) f pγ ∝ L γ t v /(R 2 ǫ γ,br (1 + z) 2 ) often used for analytical computations (assuming that t v is indicative for the shell width, see Eq. (9)). We note that this result is different from He et al. (2012, Eq. (9)) and Zhang & Kumar (2013, Eq. (6)), where the pion production efficiency scales ∝1/R. Our formula is identical to Eq. (9) in He et al. (2012) if the model-dependent relationship Eq. (8) is applied to one power of R. The next subsection shows the implications of different astrophysical models. If one defines a "magnetic loading" ξ B ≃ 1 as the ratio between energy in the magnetic field and gamma-rays, one can easily derive the magnetic field from the magnetic field energy density being equal to ξ B · u ′ γ in Eq. (7).

Implications of different astrophysical models
A review of the possible model assumptions on neutrino production models can be found in He et al. (2012); Zhang & Kumar (2013) and Bustamante et al. (2017). Here are some examples: -Internal shock model (one-zone). In this scenario, the shells of plasma are assumed to collide at a radius and the variability timescale is assumed to be indicative for the shell width All collisions are assumed to collide at the same radius in the one-zone model. As a consequence, Eq. (7) strongly depends on Γ and t v , and the pion production efficiency becomes f pγ ∝ L γ /(Γ 4 t v ǫ γ,br ) as in Guetta et al. (2004). -Internal shock model (multi-zone). In the multi-zone collision models (e.g., Kobayashi et al. 1997;Daigne & Mochkovitch 1998), shells of plasma are ejected from the central emitter, colliding at varying collision radii centered around a mean value. This case is similar to the one-zone case, but the parameters of each collision m can be very different: Γ m (Lorentz factor of merged shell), R m collision radius, ∆d ′ m (width of merged shell, related to intermittent timescale of the central emitter), and L sh γ,m (radiated energy in gamma-rays). The observed variability timescale (from the lightcurve) roughly matches the intermittent timescale of the central emitter, and thus also the shell width -as for the onezone model, but the shell width is a result of the model, rather than an input. Therefore, the key parameter determining the size of the interaction region is R, which covers a wider range than in the one-zone model, whereas the dependence on the other parameters is milder. -Photospheric models. Below the photosphere, Thomson scattering hinders the escape of gamma-rays. In photospheric A101, page 6 of 31 prompt emission models (e.g., Rees & Meszaros 2005;Giannios 2008), a fraction of energy in gamma-rays is therefore released once the radiation zone becomes optically thin to Thomson scattering, where the photon spectrum is thermal and may be different from our assumed spectrum. The production region is in this case described by the photospheric radius R ph , which is typically smaller than Eq. (8), and the densities are therefore higher. Since the shell width hardly changes in the coasting phase of the GRB, Eq. (9) should still be a good description. -Magnetic reconnection models. Magnetic reconnection models (e.g., Lyutikov et al. 2003;Zhang & Yan 2011) often describe pulses and the fast time variability in the light curves at the same time. There are therefore two variability timescales, a fast one t f v (indicative for the shell width in Eq. (9)) and a slow one t s v (indicative for the collision radius in Eq. (8)). If the slow scale is about a factor of 10−100 larger than the fast scale, R will be about this factor larger and the particle densities correspondingly smaller. The pion production efficiency can be written as Our strategy is to formulate the physics of the nuclear cascade as independently from the model as possible. Our results will therefore be shown as a function of L γ and R for one collision only, such that one can read off the physical implications as a function of the most relevant parameters. Unless stated otherwise, we use z = 2, Γ = 300 and t v = 0.01 s for the other parameters with a milder dependence. In the main text, most results will apply to the internal shock scenarios, that is, Eq. (8) holds. We show, however, in Sect. 4 how our results can be (at least qualitatively) translated to other astrophysical models.

Injection of nuclei and maximal primary energy
Nuclei of species i are assumed to be injected with a cut-off power law with a spectral index k ≃ 2 expected from Fermi shock acceleration from an acceleration zone: We use P = 2 for the cut-off function, unless noted otherwise. We note that there is basically no impact of the cut-off function on the neutrino spectra, whereas a fit to cosmic ray data will be sensitive to it. The maximal energy E ′ i,max is determined automatically by balancing the acceleration rate t ′−1 is the Larmor radius) with the sum of synchrotron loss, adiabatic cooling, photodisintegration, and photo-meson production rates, where we choose η = 1 for the acceleration efficiency. In order to define a "nuclear loading" ξ i of species i, we take into account the optically thick case to Aγ interactions (photodisintegration, photo-meson processes) and normalize the injection luminosity to the γ-ray luminosity. If the photons can escape freely, which is (at least at around the break energy, where τ γγ ≪ 1) typically a valid assumption beyond the photosphere (see, e.g., Bustamante et al. 2017), one finds to determine C ′ i in Eq. (10), with the nuclear loading ξ i , and u ′ γ from Eq. (7). As indicated earlier, we only inject one isotope at a time in this work (pure composition), where we use ξ i = 10 as frequently assumed for protons, unless noted otherwise (such as in the cosmic ray fits), and inject the most abundant stable isotope for each element (Z). With this assumption, we can inject the same luminosity for different compositions, and check how the neutrino flux depends on composition.

Neutrino production and cosmic ray escape
In addition to the nuclei, we add π + , π − , and K + mesons to the system. These are injected from photo-meson production off protons, neutrons, and all nuclear isotopes. Adiabatic cooling, synchrotron cooling, and escape (through decays) are taken explicitly into account. We also include four muon species for left-and right-handed µ + and µ − , since the helicity-dependence of the muon decays is taken into account (Lipari et al. 2007).
For the neutrinos, we have four species at the source (ν e ,ν e , ν µ , ν µ ), which receive injection from the pions and muons according to the usual decay chains, from kaon (ν µ , leading mode only), and from the β ± decays of neutrons and unstable isotopes (ν e , ν e ) from inside and outside the source. Flavor mixing between source and detection is taken into account with the mixing angles θ 12 = 33.48 • , θ 23 = 42.3 • , and θ 13 = 8.50 • , and the CP phase δ CP = 306 • (Gonzalez-Garcia et al. 2014).
The PDE system, after its setup and the pre-computation of interaction rates that do not change as a function of time in this study, is evolved until (after a few times t dyn ) the steady state is reached. For the integration, we re-parameterize the PDE system in terms of E ′2 N ′ and use a Crank-Nicolson solver.
As far as the escape of cosmic rays is concerned, we follow Baerwald et al. (2013), where the "direct" escape from a GRB shell has been discussed in greater detail 3 . It is conservatively assumed that charged particles can only escape from within the Larmor radius of the edges of the shells. This means that a fraction f esc = min(R ′ L (E ′ ), ∆d ′ )/∆d ′ ≤ 1 escapes over the dynamical timescale. We note that f esc = 1 if the Larmor radius reaches (or exceeds) the shell width, which corresponds to the free-streaming limit. It can be shown that f esc is invariant even if the shell expands (Baerwald et al. 2013). The escape fraction can be translated into an effective escape rate . It is interesting to compare that rate to the escape from the whole shell: . For neutrons, the escape rate is given by the free-streaming rate (see Appendix C for a comparison to the model in Baerwald et al. 2013).
As a consequence, protons and nuclei can escape with a rate ∝R ′ L ∝ E ′ and correspondingly hard spectra, which means that they can escape freely when the Larmor radius reaches the shell width. It can be shown that this condition is satisfied at the highest energy for η = 1 if the maximal primary energy is limited by the dynamical timescale. This situation can be typically found in sources where the radiation densities and the maximal primary energy are low. If, on the other hand, the radiation densities are high such that the maximal energy is limited by the photohadronic interactions, the escape at the highest energy will be suppressed.

Propagation of cosmic rays
Once the accelerated particles are able to escape from the source, they propagate through the extragalactic space, encountering photon fields, such as the cosmic microwave background (CMB) and the extragalactic background light (EBL, from infrared to ultraviolet), with which they interact. Similarly to what happens in the source, the energy scale of the process is given by the photon energy in the nucleus rest frame. The dominant processes for nuclei are photodisintegration on CMB at the highest energies and on EBL at intermediate energies, while the photo-meson production is shifted towards A times the threshold for protons. At the lowest energies, nuclei lose energy adiabatically because of the expansion of the Universe. For protons, the intermediate energy range is dominated by energy losses by electron-positron pair production on the CMB. Interactions of cosmic rays with the EBL is more relevant for nuclei than for protons. However, the EBL provides the majority of lowenergy cosmogenic neutrinos and it is, therefore, included in our calculations.
The propagation of cosmic rays ejected from the GRBs is computed using the simulation code SimProp (Aloisio et al. 2012); the isotopes produced during interactions in the source depend on the adopted photodisintegration model. In the present paper we use the SimProp code with the Puget-Stecker-Bredekamp (PSB) model for photodisintegration processes (Puget et al. 1976;Stecker & Salamon 1999), that contains one representative isotope for each nuclear mass. The nuclei ejected from the source are then grouped and the corresponding fluxes are summed, assigning the sum to a mass chosen as representative for each group. The ejected particles are propagated through the extragalactic space: for this study we use the Gilmore extragalactic background light (EBL; Gilmore et al. 2012), that is one of the models available in SimProp. The results presented in the following are dependent on the choices of the physical quantities previously described, such as the EBL and the cross-section models. A detailed description of the effects of using different models for these quantities is given in Alves Batista et al. (2015), comparing the expected observables at Earth with the cosmic ray data. Here, we use the models that correspond to the best fit of the Pierre Auger Observatory data given in Aab et al. (2017). We assume a homogeneous distribution of identical sources up to z = 6. The events are simulated with a uniform distribution of log 10 (E inj /eV) and a uniform distribution of sources in the cosmologically co-moving frame. Each event is then weighted with a source evolution equal to the star formation rate, as given for example in Yuksel et al. (2008). An order four tensor is computed, containing the number of nuclei arriving at Earth with a given mass number and energy, which were ejected with a certain mass number and energy from the sources. The matrix is multiplied with vectors containing the spectrum of the representative mass ejected from the source, in order to obtain the spectrum at Earth. Neutrinos ejected from the sources and those produced during the propagation are treated similarly. The former ones are computed with NeuCosmA following the procedures described in previous chapters. A bi-dimensional propagation matrix is used to apply weights to these fluxes according to the choice of the source evolution. The latter ones are computed with SimProp and their flux at Earth is computed using equal weights as for the spectra of cosmic rays ejected from the sources. The all-particle cosmic ray spectrum at Earth is normalized to the energy spectrum measured by the Pierre Auger Observatory (Valiño et al. 2015). The baryonic loading is computed as in Baerwald et al. (2015).

Comparison with existing computations and approximations
Compared to the computation in Globus et al. (2015a), the method presented here has the advantage that it is fully deterministic and fast (which allows for the parameter space scans we perform), but the disadvantage that certain effects cannot be as easily included as in a Monte Carlo simulation. One example is the re-acceleration of secondary nuclei, which has been included in Globus et al. (2015a). As discussed in Winter et al. (2014), for the re-acceleration of secondary mesons, a prerequiste for the re-acceleration is an efficient transport back to the shock, which can be expected at the highest energies when the Larmor radius is comparable to the distance to the shock front. This means that the secondaries have to be produced at the first place (i.e., disintegration is efficient) but the maximal energy is limited by the size of the region (i.e., disintegration and photo-meson production are inefficient). For nuclei close to the primary, significant effects can therefore be only expected in a small fraction of the parameter space, whereas the secondary proton (and light nuclei) spectra might be slightly enhanced. Our overall features, however, agree with Globus et al. (2015a). Furthermore, we note that the resulting spectra in that reference are a superposition of multiple collisions, which also has an impact on the final spectra.

Nuclear cascade source classes
In order to describe the nuclear cascade in the sources and the characteristics of the ejected nuclear isotopes from a GRB shell (interaction region), we introduce three qualitatively different cases in this section: "Empty Cascade", meaning that the nuclear cascade hardly develops, "Populated Cascade", meaning that the nuclear cascade develops, and "Optically Thick Case", meaning that the source is optically thick to Aγ interactions for nucleons and all nuclei. While there are continuous transitions among these cases, we discuss three prototypes, which exhibit the behavior characteristic for their class. In order to keep the discussion as simple as possible, we inject a pure composition of 56 Fe only, which is the heaviest commonly discussed injection isotope. We will show continuous parameter space scans in the following sections.

Empty cascade
We define the Empty Cascade as the case that is optically thin to Aγ interactions of the injection isotope (and, consequently, all lighter isotopes including nucleons). We have checked that this definition is almost equivalent to asking for the neutrino flux to be dominated by photo-meson production off the injection isotopes. The luminosity of this source class is relatively low (or the production radius is large), since low enough radiation densities in the source are a necessary requirement for that. The rates of the dominant processes (apart from Bethe-Heitler pair production) for primary nuclei and photo-meson production off protons t ′−1 pγ are shown in the upper left panel of Fig. 4 for one example. By comparing t ′−1 Fe γ and t ′−1 dis with t ′−1 dyn ≃ t ′−1 ad ≡ c/∆d ′ (see arrows), one can clearly see that the source is optically thin to Aγ interactions at the highest energy (dominated by the dynamical timescale or adiabatic losses).
The particle spectra in the source at the end of the time evolution are shown in the upper right panel of Fig. 4. The primary (E −2 ) injection spectrum of 56 Fe is indeed hardly modified and extends to the mentioned maximal energy. Secondary . Example for the "Empty Cascade" source class (isotropic luminosity L γ = 10 49 erg s −1 ) for injection of pure 56 Fe. Interaction rates (top left), particle densities in the source (top right), nuclear cascade (bottom left), and ejected cosmic ray fluence per shell (bottom right, without CMB and EBL interactions) as a function of the energy in the observer's frame. The different curves in the right panels correspond to the different isotopes: primary 56 Fe in blue, secondaries produced in nuclear cascade according to legend, where the result from secondary nuclei other than nucleons are summed over. The other GRB parameters are chosen to be R ≃ 10 8.3 km, Γ = 300, ξ Fe = 10, ε ′ γ,br = 1 keV, and z = 2.
nuclei and nucleons are suppressed, and their maximal energies follow basically the conserved Lorentz factor. The spectra are somewhat harder than E −2 because of the high-energy behavior of the Aγ cross section (which is not peaked like a resonance there). The nuclear cascade (integrated energy of isotopes relative to total injection energy) is shown in the lower left panel of Fig. 4. One can clearly see that apart from 56 Fe a few close by isotopes are populated, while most of the nuclear cascade is, relative to the injection luminosity, empty. Nucleons and light nuclei (especially 4 He) are produced as disintegration products, but their occupation is relatively small compared to the injection isotope. Figure 4, lower right panel, illustrates the ejected cosmic ray spectra for the assumptions stated earlier. This type of figure shows the fluence at Earth including adiabatic losses, but no interactions with the CMB and EBL (which deform the spectrum). Here, charged cosmic rays escape via direct escape, while neutrons are not magnetically confined (see Baerwald et al. 2013, for a more detailed discussion). This leads to a characteristically different shape of the ejected spectra: relatively hard ejection spectra for the charged cosmic rays, and softer spectra for the neutrons, which will eventually decay into protons on their way to Earth. There is evidence that such hard ejection spectra are actually the required input for the cosmic ray propagation of UHECR nuclei to describe Auger data (Aab et al. 2017). We will come back to this issue later.
A characteristic of the Empty Cascade is that the cosmic ray ejection spectrum is dominated by the hard spectrum of the injected primaries. Nevertheless, the contribution of the neutrons is, especially at low energies, more substantial than one may think from the upper right panel in the source. The reason is that the neutrons basically escape unmodified, while the nuclei are somewhat depleted even at the highest energies by the escape mechanism (see

Populated cascade
We define the Populated Cascade as the case that is optically thick to Aγ interactions of the injection isotope that will disintegrate and populate the cascade. At the same time, the source environment is optically thin to pγ interactions, so that the proton and neutron fluxes will be hardly affected by these interactions (see arrows in Fig. 5). We have checked that this definition is almost equivalent to asking for the neutrino flux to be dominated by photo-meson production off the secondary isotopes lighter than the primary, but heavier than the nucleons. The luminosity of this source class is intermediate. The source class corresponds to the typical assumption of the optically thin case for protons, extended to nuclei. Figure 5 shows one example for this source class in the same format as Fig. 4. From the upper left panel, we can indeed see that the source is optically thick to Feγ interactions, while it is optically thin to pγ interactions. The maximal primary energy is given by Aγ interactions. The nuclear cascade (lower left panel) is, relative to the injection energy, well populated, where light nucleons are occupied similar to the near-56 Fe isotopes. A101, page 10 of 31 From the densities in the source (upper right panel) one can see a clear depletion of the injection E −2 spectrum beyond the disintegration threshold, while the secondary isotopes are populated with a density comparable to the original primary density. The proton and neutron densities are larger than in the Empty Cascade, but not yet comparable to the injection density.
Nevertheless, the ejected cosmic ray spectrum (lower right panel) is already dominated by the neutrons. The reason is that the maximal energy is dominated by Aγ interactions and therefore the Larmor radius can reach only about 1/30 of the size of the region at the maximal energy (see upper left panel: compare acceleration and dynamical rates at the maximal primary energy). Consequently, the ejected secondary spectra are suppressed by a similar factor at the maximal energy because cosmic rays from the innermost regions of the shell take too long to escape. The effective (all energy) ejection spectrum including neutrons and nuclei will be relatively soft because of the neutron escape component, with a potential dip.
Comparing the Populated and Empty Cascades, we notice that the isotropic luminosity, which determines the photon density in the source, controls the relative height between the neutron and nuclei escape components. Since we use the isotropic luminosity as fit parameter later, the right luminosity to UHECR data will be automatically determined later.

Optically thick case
We define the Optically Thick Case as the one that is optically thick to Aγ interactions of nucleons and nuclei. This can be clearly seen in the upper left panel of Fig. 6, which shows one example with an extremely high luminosity (arrows). It turns out that the neutrino production in this case will be dominated by the protons and neutrons produced in the disintegration chain. It is first of all instructive to compare the nuclear cascade in the lower left panel of Fig. 6 with the corresponding one in Fig. 5. The cascade appears to be less populated off the main diagonal because the intensity of the cascade is actually reduced at intermediate isotopes, and most energy is dumped into nucleons produced in the cascade; the nucleons are populated similarly to the primaries. The maximal primary energy is determined by photodisintegration (see upper left panel).
The primary and secondary isotope densities are strongly suppressed beyond the photodisintegration threshold (see upper right panel), while the protons and neutrons are populated at a level comparable to the injected primary flux. Compared to the Populated Cascade, the nucleons peak at the photo-meson threshold (break in t ′−1 pγ in upper left panel), because they cascade down in energy by multiple interactions.
The ejected charged cosmic rays (lower right panel) dominate at the highest energies but are strongly suppressed because of the low densities in the source and because the Larmor radius at the highest energy is extremely small compared to the size of the region. In fact, it is difficult to reach the UHECR energy range because of the Aγ interactions limiting the maximal primary energy.

Impact of astrophysical parameters and models
In the previous section, we discussed the results for three prototypes from the Empty Cascade, Populated Cascade, and Optically Thick Case regimes. We follow the previous classification using the interaction rates at the maximal energy, and show in Fig. 7 these three regimes as a function of L and R for the injection of pure oxygen (left panel) and iron (right panel). The transition between the regions in terms of the development of the nuclear cascade is continuous, while our classification scheme will produce well-defined different regimes. In Fig. 7 we also show the sub-photospheric region, that is, the region from which photons cannot escape because of Thomson scattering 4 .
Since high luminosities or small collision radii mean high photon densities, it is clear that the optically thick case and even the photosphere are reached in the lower right corners of the panels, whereas the cascade is hardly populated in the upper left corner. Our chosen prototypes (circles in right panel) correspond to extreme examples for the Empty Cascade and Optically Thick Case, whereas the prototype for the Populated Cascade lies well within the corresponding region. The maximal primary energy (gray dashed contours) is typically given by Aγ interactions (lower right) or adiabatic cooling (upper left) 5 . The adiabatic cooling limited case corresponds to a rigidity-dependent maximal energy, and it roughly coincides with the Empty Cascade case.
The transition between the Empty and Populated Cascades depends on the injection composition (compare left and right panels of Fig. 7) because the Aγ interaction rates increase with 4 Assuming that the plasma is electron-(fully stripped) isotopedominated and that the electron density is similar to the proton density (from charge conservation, electron-positron pair production assumed to be sub-dominant), the photospheric radius -where the optical thickness to Thomson scattering is one -is given by R ph ≃ (1/2 ξ i L γ t v σ Th /(4 π κ (1 + z)Γc 2 m p )) 0.5 . Here κ ≃ 0.25 (see e.g. last column, Table 2 in Bustamante et al. 2017, for simulation results) is the re-conversion efficiency from the kinetic into the radiated (nuclei dominated) energy. Note that the factor 1/2 comes from the assumption of heavier nuclei (which carry about as many protons as neutrons), which descreases the photospheric radius by 1/ √ 2. 5 Since the interaction rates t ′−1 Aγ ∝ u ′ γ ∝ L γ /R 2 (cf., Eq. (7)) and t ′−1 acc ∝ B ′ /E ′ ∝ u ′ γ /E ′ ∝ L γ /(E ′ R), one finds the maximal primary energy to be constant along the curves R ∝ L γ in the Aγ dominated region. If, however, adiabatic losses dominate the maximal energy, the dependence on R cancels and the maximal energy hardly depends on R because of the implied internal shock relation Eq. (8) leading to t ′−1 ad = 2Γ/R.  Fig. 7 for 56 Fe, but the internal shock model relationship Eq. (8) among collision radius, t v , and Γ is not applied, and only holds in the region marked "Internal Shock model", which implies that for the prototypes (circles) the same results as in Sect. 3 are obtained. Here the idea is that R is the main control parameter for the target density in the shells, whereas t v determines the shell thickness and Γ the energy shift only, that is, the dependence on these parameters is milder (the pion production efficiency scales ∝t v similar to L γ , but typically does not vary in such a large range; therefore L γ is chosen as parameter here); further details are given in Sect. 2.5. The regions applying to different astrophysical production models ( mass number, which determine this classification. Correspondingly, the transition between adiabatic and interaction dominated maximal energy shifts (see gray dashed contours). Since the transition between Populated Cascade and Optically Thick Case depends on the optical thickness to proton and neutron interactions, it is not affected by the injection composition. Therefore, the lower the injection mass is the smaller the Populated Cascade region will be, because the Aγ interaction rate increasingly approaches the pγ rate. In order to address the astrophysical model-dependence of the nuclear cascade, Fig. 8 shows the same result as Fig. 7 for 56 Fe with t v = 0.01 s fixed (instead of scaling it with the internal shock model relationship Eq. (8) among collision radius, t v , and Γ). We note that the pion production efficiency f pγ ∝ L γ t v /R 2 (see Sect. 2), which means that the dependence on Γ (which shifts the energies) is very mild, and R is the dominant control parameter, which we have chosen on the vertical axis compared to earlier papers Baerwald et al. (2015).
The maximal proton energy behaves similarly to Fig. 7 in the lower right corner, whereas in the adiabatic-loss-dominated upper left corner, the scaling is the same as in the lower right corner. Since t ′−1 acc ∝ B ′ /E ′ ∝ u ′ γ /E ′ ∝ L γ /(E ′ R) is proportional to the (constant) adiabatic cooling rate, one has R ∝ L γ for fixed maximal energy. This observation is interesting because the cosmic ray fit is sensitive to the maximal energy: since maximal energy and Aγ interaction rates scale in the same way in the parameter space shown, we do not expect significant changes of the results if one moves parallel to the lines (for the fixed chosen value of t v ). Consequently, our protoype behavior of the nuclear cascade (compare to the dots in the figure) does not change very much along these lines.
The figure also displays typical regions expected for different prompt emission models, as outlined in the figure caption and discussed in Sect. 2. For example, larger collision radii are expected in the ICMart model compared to the internal shock model because the pulse timescale is indicative for the collision radius. While Eq. (8) holds for the internal shock scenario, the chosen Γ and t v lead to a relatively narrow region for the internal shock model, while the multi-collision model has a large intrinsic spread of the collision radii around this nominal value. Photospheric models typically predict a prompt emission dominated by the photosphere. One can read off from Fig. 8 what happens for a certain choice of L and R, such as for a specific collision in the multi-zone model. Moving along the "diagonals", one can then easily find a prototype that displays the qualitative behavior. For example, the result of our L = 10 49 erg s −1 , R = 10 8.3 km internal shock model prototype will correspond to a L = 10 51 erg s −1 , R = 10 9.3 km ICMart collision in terms of nuclear cascade development, maximal primary energy, and neutrino production efficiency. However, the ICMart collision has higher γ-ray and proton luminosities, which means that in the cosmic ray fit a correspondingly lower baryonic loading will be required.

Prompt neutrino production
Let us now study the neutrino production for the prototypes in Sect. 3. Figure 9 shows the all-flavor neutrino fluence per shell for the injection of pure 56 Fe, split up by the contribution from the injected primary, all secondaries, and the interactions of nucleons produced in the cascade. As indicated earlier, we have checked that the classification of the three regimes (Empty Cascade, Populated Cascade, Optically Thick Case) roughly corresponds to a classification according to the dominant contribution of the neutrino fluence: For the Empty Cascade (see upper left panel), the injected primary dominates the neutrino production; for the Populated Cascade, the contribution from all secondary isotopes produced in the nuclear cascade dominates; and for the Optically Thick Case, the nucleons produced by nuclear disintegration dominate. The contribution from neutron decays is visible as a bump at low energies in all cases, but it never dominates the neutrino fluence. The neutrino fluence grows quadratically with luminosity since photon and baryon density each scale with luminosity.
This separation of the contributions to the neutrino fluence is not only for conceptual reasons, but has more profound implications. As indicated in Sect. 2, the photo-meson production off nucleons is relatively well known and has been studied in detail in the past with the SOPHIA interaction model. However, the photo-meson production off nuclei relies typically on a superposition model, which (in our case) scales the cross section ∝A. Therefore we know that the neutrino fluence prediction off nucleons (red curves in Fig. 9) is more or less robust, whereas the other contributions carry large uncertainties to be quantified in the future. As a consequence, the neutrino prediction A101, page 14 of 31 for the Populated Cascade (where the nucleon contribution is sub-dominant but large) and Optically Thick cases (where the nucleon contribution dominates) is more reliable than in the Empty Cascade case -where, however, the neutrino flux is low anyway.
The impact of the injection isotope (for pure composition injection) is shown in Fig. 10. That is, we inject the same luminosity of different isotopes, and we show the resulting neutrino fluence for pure proton, pure 56 Fe, and pure other (intermediate injection isotope, shaded region) in the figure for the different prototypes. As one of our most important results, the neutrino fluence at the peak hardly depends on the injection composition. This result is a consequence of several factors: Lorentz factor conservation in the disintegration, which splits up the nucleus into lighter nuclei or nucleons, the E −2 injection spectrum, which conserves the energy per decade, the photomeson interaction rate being almost flat beyond the threshold, and strong magnetic field effects on the secondary muons, pions, and kaons, which determine the maximal energy cut-off. One can, for instance, treat the nucleus as a superposition of nucleons and re-write the secondary production in terms of the nucleons (see, e.g., Joshi et al. 2014, for nucleus-nucleus interactions). Then one can see that the secondary production is roughly the same as before if the primary flux × interaction rates roughly scale ∝E −2 and the photo-meson cross section σ Aγ ≃ Aσ pγ . Because of the Lorentz factor conservation in disintegration (there is effectively hardly any energy lost), disintegration does not affect this result.
Regarding the shape of the spectrum, the high-energy cut-off is somewhat higher for lighter injection isotopes because the maximal primary energy does not follow the Peters cycle Peters (1961, rigidity-dependent maximal energy) for the Populated Cascade and Optically Thick cases (as we will discuss in Sect. 7). This means that the maximal energy per nucleon, which /erg/s) Here χ 2 − χ 2 min is shown as a function of R and L γ for the fit to cosmic ray data of the Pierre Auger Observatory (Valiño et al. 2015;Aab et al. 2014) above 10 18 eV; pure silicon injection with k = 1.8 is assumed. The sources are distributed in the range from z = 0 to 6 following the star formation rate (SFR). The blue squares mark the current (90% C.L.) IceCube-excluded region from the GRB stacking analysis from northern and southern sky muon tracks Aartsen et al. (2015Aartsen et al. ( , 2017, while the green squares mark the current (90% C.L.) IceCube-excluded region from the cosmogenic neutrino analysis (Aartsen et al. 2016), applied to ν µ +ν µ . The contours show the nuclear loading (log 10 ξ).
affects the maximal neutrino energy, actually decreases. If this effect is stronger than the cut-off from magnetic field effects on the secondaries, it becomes visible in the neutrino fluence. Furthermore, the neutrino fluence from nuclei exhibits a stronger low energy component, which mostly comes from neutron decays produced in the nuclear cascade.

Description of cosmic ray data
In order to connect with cosmic ray data, we extrapolate now from the previously studied single collision zone to a population of alike GRBs with a fixed duration of ∆T = 10 s, which implies that the emission from each GRB comes from ∆T/t v such collisions. We perform a fit of UHECR observations from the Pierre Auger Observatory, combining the modeling of interactions in the source (computed as described in the previous sections) with the propagation of cosmic rays in the extragalactic space. For the propagation, the SimProp code (Aloisio et al. 2012) has been used with Gilmore EBL (Gilmore et al. 2012) and PSB crosssection model (as defined in Alves Batista et al. 2015). We also consider the extensive air shower in the Earth's atmosphere and EPOS-LHC (Pierog et al. 2015) is assumed for UHECR-air interactions. We find that a good description of the data is obtained by distributing the GRBs as sources of cosmic rays following the star formation rate Yuksel et al. (2008), assuming a pure 28 Si at the injection. The primary spectrum is described in Eq. (10) and we use here k = 1.8 and P = 2. We fix the following parameters: source evolution, spectral index and cutoff shape at injection, and the nuclear species at the injection. In the present work, we keep this procedure as simple as possible in order to show the power of the method, while leaving a more detailed analysis for a future work. The fit is performed above 10 18 eV (Mixed Composition Dip Model) and above 10 19 eV (Mixed Composition Ankle Model) by using the combined spectrum (Valiño et al. 2015) and the shower depth (X max ) distributions (Aab et al. 2014), which contain information about the mass of the nucleus interacting with the atmosphere, with a similar procedure as used in di  and Aab et al. (2017). A scan over (R, L γ ) is performed and for each pair the normalization to the experimental flux is found. For each point of the parameter space the number of expected prompt and cosmogenic neutrino events is calculated following Baerwald et al. (2015). For the prompt neutrino flux, the exposure for muon neutrinos is calculated by summing the exposure relative to the IceCube analyses of Aartsen et al. (2015, 3 yr) and that of Aartsen et al. (2017, 3 yr) for a total of 1014 GRBs that occurred in the Northern Hemisphere, and that of Aartsen et al. (2017, 5 yr) for 664 GRBs that occurred in the Southern Hemisphere, and comparing the total number of bursts of the combined sample with the assumed 667 bursts per year as in Baerwald et al. (2015). For the cosmogenic neutrino flux, the exposure for muon neutrinos is taken from Aartsen et al. (2016). The exclusion regions (90% C.L.) of the prompt and cosmogenic neutrino analyses are calculated assuming both analyses as background free.

Mixed Composition Dip model
We first fit the spectrum and composition above 10 18 eV, including the ankle region. The result of the fit is shown in Fig. 11, where the χ 2 − χ 2 min as a function of R and L γ is given. The blue region indicates the current IceCube exclusion region from Cosmic ray observables obtained with the best-fit parameters in the Mixed Composition Dip Model, corresponding to point A for (log 10 (R/km), log 10 (L γ /(erg s −1 ))) = (8.1,50.5) in Fig. 11. Top: simulated energy spectrum of UHECR (multiplied by E 3 ) compared to data from Valiño et al. (2015). Spectra at Earth are grouped according to the mass number as follows: A = 1 (red), 2 ≤ A ≤ 4 (gray), 5 ≤ A ≤ 22 (green), 23 ≤ A ≤ 28 (cyan), total (brown). Bottom: average and standard deviation of the X max distribution as predicted (assuming EPOS-LHC, Pierog et al. 2015, for UHECR-air interactions) for the model versus pure ( 1 H (red), 4 He (grey), 14 N (green), and 56 Fe (blue)) compositions, compared to data from Porcelli et al. (2015).
the prompt GRB analysis Aartsen et al. (2017), and the green marks refer to the cosmogenic neutrinos limits Aartsen et al. (2016), both at the 90% confidence level (C.L.). The best-fit spectrum at Earth and the composition observables are shown in Fig. 12. The ankle feature is well reproduced without the need of additional components, since the interactions in the source naturally produce a light component at the lowest energies and a heavy one at the highest energies, allowing for the reproduction of the composition trend. However, we emphasize that the current procedure is not optimized to find the best description of the data, for which at least a mixed composition at the injection and a shift in the energy scale would be required. The contours of the baryonic loading in Fig. 11 follow the behavior of the maximal energy, confirming what has been found in Baerwald et al. (2015). The baryonic loading reaches extreme values in the region of the photosphere. The baryonic loading required at the best fit is 3 × 10 4 . The best fit clearly lies in the Populated Cascade region, where the prompt neutrino flux is dominated by photo-meson production off the secondary isotopes in the cascade. This region is excluded by the GRB stacking analysis. In the Optically Thick region, the amount of nucleons emitted at the source is maximal due to the development of the cascade; as a consequence, this region is also excluded by the fit because of the expected cosmogenic neutrino flux. Since a number of parameters are fixed in our calculation, this observation should not be perceived as being too general. Moreover, the goodness of the fit including the ankle is low due to the very small statistical errors at the lowest energies and to the fact that in the current work the systematic uncertainties are not included. In Fig. 13, the cosmic-ray spectra and composition observables are shown for selected points. Point A corresponds to the best fit, while points B and C correspond to cases with the same collision radius as the best-fit case, with a different choice for the luminosity. The B case is at the border between the Empty and Populated Cascade. This can be seen also in the lower energy part of the energy spectrum in Fig. 13, where the amount of propagated nucleons cannot account for the flux. The opposite case is represented by point C, where high photodisintegration in the source results in a large contribution from nucleons at lower energies, which further increases during propagation in the extragalactic space. The composition observables (lower panels of Fig. 13) react on increasing luminosity at the source with a sharper transition from light to heavy masses as energy rises. In Fig. 14  prompt neutrino flux (left panel), but it is still allowed within the cosmogenic neutrino limits (right panel). The prompt flux clearly shows the discussed enhancement as a function of the luminosity of the sources. The cosmogenic flux reaches the IceCube limit only when approaching the Optically Thick Case (C), because of the amount of nucleons injected into the extragalactic space.

Mixed Composition Ankle Model
The fit is performed above 10 19 eV, with a penalty for overshooting the flux at lower energies. The goodness of fit is here enhanced with respect to the Mixed Composition Dip Model due to the absence of the data points at the lowest energies. The best fit is found at low source luminosity and intermediate collision radius (point A in Fig. 15). This result is attributable to the requirement of avoiding the overshooting of the flux at the lowest energies, which can be naturally achieved in the case of the Empty Cascade. Moreover, the position of the best fit is also slightly dependent on the choice of the minimal energy for the fit: while lowering the starting energy for the fit, the flux at the lowest energies has to be enhanced, pushing the preferred solution towards the Populated Cascade Region. The parameters corresponding to the best fit are not excluded by the existent neutrino limits. This case corresponds to the Empty Cascade: the protons in the energy spectrum at Earth (Fig. 16, upper panel) come, therefore, mainly from propagation in the extragalactic space. The obtained GRB parameters require either very low γray luminosities (points A and B), which are significantly lower than expected from γ-ray observations (see Atteia et al. 2017; for an overview), or large collision radii (point C), which may be indicative for magnetic reconnection models. The transition of the mass composition from light to heavy is less sharp than in the best fit for the Mixed Composition Dip Model (compare the lower panels of Figs. 12 and. 16). We also show in Fig. 17 the energy spectra and the composition observables for other selected points in Fig. 15. Point B reproduces the cosmic ray spectrum in a reasonable way, but with a sharper transition from light to heavy masses. The left panel of Fig. 18 shows that enhancing the photodisintegration with a higher source luminosity results in an overshooting of the prompt neutrino limits, while the cosmogenic neutrino flux does not substantially change before reaching the Optically Thick Case. We also show point C, which represents an intermediate luminosity and a high collision radius. In this case, photodisintegration in the source is not efficient and results in a high maximal energy of the primary at the source (compare with Fig. 7), and propagation will lead to very energetic protons at Earth. This can be seen in Fig. 18, where σ(X max ) is almost flat. The same message is given by the corresponding cosmogenic flux shown in the right panel of Fig. 18.  Fig. 11 for the Mixed Composition Dip model. Here we compare to the differential limits obtained from Aartsen et al. (2017, considering northern + southern exposure) and Aartsen et al. (2016). The different points are: A: (log 10 (R/km), log 10 (L γ /(erg s −1 ))) = (8.1, 50.5), B: (8.1, 49.5), C: (8.1, 51.7). Here the differential limits are defined as in Baerwald et al. (2015), such that following the differential limit curve for one decade in energy will yield one event.
/erg/s) Here χ 2 − χ 2 min is shown as a function of R and L γ for the fit to cosmic ray data of the Pierre Auger Observatory (Valiño et al. 2015;Aab et al. 2014) above 10 19 ev, including a penalty for the overshooting of the flux at lower energies. The caption of Fig. 11 gives details.

Impact of injection composition
So far, we have tested several selected injection isotopes: 16 O and 56 Fe to describe the disintegration in a GRB shell as kind of extreme examples for nuclei, and 28 Si for a reasonable fit to Auger data. In this section, we qualitatively discuss the dependence on the injection element (Z) in order to illustrate what changes to expect for different injection isotopes, why 28 Si is a good example to describe Auger data, and why the best-fit parameters are found in the regions determined in the previous section. We will show the results as a function of Z, implying that we inject the most abundant stable isotope.
We anticipate that the UHECR ejection spectrum from the sources (injection spectrum into the intergalactic medium) can be qualitatively characterized by four different estimators (see Fig. 19 for illustration): -Spectral shape. We assume that the ejection spectrum can be roughly described by a power law with a certain spectral index. Since both the neutron and charged cosmic rays contribute to the ejected spectrum, the overall spectrum will be roughly determined by the peaks of these two contributions (see left panel). -Maximal energy. The ejection spectrum will have a characteristic maximal energy E max . This maximal energy is determined by the equilibrium between acceleration rate and the sum of energy loss rates (see Sect. 2.6). -Composition at E max . The composition at the maximal energy will be relevant to describe the observed composition at the highest energy (see right panel for illustration). -Transition energy to heavier composition. If the ankle is to be described by the GRBs, the transition energy E trans has to be (E/eV) in the right ballpark. It is determined from the second derivative of the composition curve (see right panel). We do not discuss the cosmological source evolution as a separate parameter here, although it is well known that there is a certain degeneracy between source evolution and spectral ejection index.
The four estimators are shown as a function of Z for different values of L γ in Fig. 20. The estimated preferred ranges from cosmic ray data are shaded, that is, we expect that the estimator value lies in the indicated band to describe data. In some of the curves, the zig-zag pattern comes from the mass pattern of the most abundant injection isotopes as a function of Z (see Fig. 2; dark blue isotopes). We note that R is fixed in this figure.
From the upper left panel of Fig. 20 one can see that the ejected spectral index prefers 10 49 erg s −1 L γ 10 51 erg s −1 , while the Optically Thick case results in ejected spectra which are too soft (see example in Fig. 6). This result hardly depends on the injection composition Z. The ejected cosmic rays have hard spectral indices for low GRB luminosities, as found in our fit of the Mixed Composition Ankle Model; this is in agreement with what has been found in Aab et al. (2017), where the spectra at injection are comparable to what we find here as ejected from the source.
From the upper right panel, the maximal energy tends to be too low for large luminosities, where the maximal energy is strongly limited by Aγ interactions, while almost all other cases lead to maximal energies in a reasonable range (as long as Z 8). One can also see from this panel, when the often used Peters cycle is justified (as in Aab et al. 2017): in that case, the maximal energy scales with rigidity. For low luminosities, corresponding to the Empty Cascade, we have found that the maximal energy is indeed limited by adiabatic losses, which is identical to this assumption. Therefore, the 10 49 erg s −1 line indeed follows the Peters cycle (see dashed curve for comparison). However, at higher luminosities (as one would typically expect for GRBs), the maximal energy becomes limited by Aγ interactions, and scales much more mildly with Z.
The lower left panel of Fig. 20 shows the ejected composition at the maximal energy, which needs to be at least as heavy as the composition observed by Auger at the highest energy (see, e.g., the lower panels of Fig. 12). While it is clear that the composition at the highest energy is somehow proportional to Z, the ejected ln A depends on the distribution of secondaries in the disintegration chain if the maximal energy does not follow the Peters cycle, which reduces the average mass for higher luminosities. Nevertheless, the dependence on L γ is very mild, while there is an obvious strong dependence on Z.  Fig. 15 for the Mixed Composition Ankle Model. Here A: (log 10 (R/km), log 10 (L γ /(erg s −1 ))) = (8.1, 49), B: (7.8, 50), and C: (9.6, 51). relatively stable in Z, which means that this estimator is the only one sensitive to the injection composition in that range.
Taking the information from the first three estimators, we can now qualitatively explain the fit range 10 49 erg s −1 L γ 10 51 erg s −1 in Fig. 15 for the Mixed Composition Ankle Model. The contours in the two-dimensional space roughly follow the maximal energy contours (similar to Fig. 7), apart from the changes due to the penalty for the overshooting at the lowest energies. The region preferred by the fit is found for intermediate values of the maximal energy, with low isotropic luminosity in order to avoid the overshooting at the lowest energies. The isolated region at high collision radius and intermediate luminosity corresponds to the highest values for the maximal energy of the residual primaries in the source in the Empty Cascade region. Fig. 19. Illustration of quantities used for qualitative analysis using one example. Left panel: ejected cosmic ray spectrum as a function of (observer's frame) energy. Right panel: ejected composition ln A as a function of (observer's frame) energy. Here adiabatic losses are taken into account, but no interactions in CMB and EBL. The left panel illustrates how to determine the spectral index, the right panel the composition at E max and the transition energy E trans . The maximum energy E max is determined by acceleration and loss processes.
In the context of the Mixed Composition Dip Model, the best fit is found corresponding with an isotropic luminosity for which the transition energy is allowed (see bottom right panel in Fig. 20). We discuss this issue in somewhat greater detail in Fig. 21, where the lnA is shown for our three prototypes (Empty Cascade, Populated Cascade, Optically Thick case). For the Empty Cascade, we know that the primary dominates at the highest energy (see Fig. 4), which leads to a relatively smooth transition at relatively low energies, which slightly increases with Z. For the Populated Cascade, the transition hardly depends on the injection composition and only slightly decreases for higher Z. The transition is relatively sudden and becomes smoother for higher L γ , where the cascade increasingly populates the lighter isotopes (see Optically Thick Case). We note that the example shown in Boncioli et al. (2017) has an L γ between the middle and right panels.
Finally, it is difficult to compare the injected Z to theoretical expectations, because these depend on where the jet originates, how charged cosmic rays are injected into the jet, and how they are accelerated. For example, if the jet originates within a Wolf-Rayet star, there may be a significant contribution of He, O, and C. If there is a close connection with supernova explosions, a significant contribution of Si or Fe may be expected. If r-process nucleosynthesis occurs within the jet, even heavier elements may contribute, as for example studied in Metzger et al. (2011). Our approach follows the opposite direction: we trace back the observed UHECR composition into the acceleration zone and identify what is needed there to describe UHECR data. In the future, one would expect that this information will be useful for models of jet formation and injection of UHECR nuclei.

Conclusions
In this study, the question of whether GRBs can be the sources of the UHECRs has been addressed in the multi-messenger context. While the neutrino emission from the prompt phase of GRBs has been used to constrain GRB models in stacking analyses, all these constraints have been derived for proton primaries. However, the UHECR composition measured by the Auger experiment indicates a dominant contribution of nuclei at the highest energies. We have, therefore, studied the UHECR paradigm for GRBs in the presence of nuclei in the sources, by taking into account the nuclear cascade in the sources with a high level of detail, and by combining the UHECR source and propagation physics.
As a first step, we injected nuclei with a pure composition into a GRB shell, and we quantified the nuclear cascade, the emission of cosmic rays, and the production of neutrinos. We identified three qualitatively different regimes as a function of the shell parameters, such as production radius R or gamma-ray luminosity L γ : -Empty Cascade (low L γ or large R). The source is optically thin to photodisintegration of primary (injected) nuclei, where the nuclear cascade does not develop. The ejected cosmic rays are dominated by a hard spectrum of the residual primaries, where (as a function of injection charge) the cosmic rays follow a Peters cycle with a rigidity-dependent maximal energy. This scenario is consistent with the global fit result, recently performed by the Auger Collaboration Aab et al. (2017). The neutrino flux is dominated by photo-meson production off the primary nuclei; it is relatively low because of the low target photon density. -Populated Cascade (medium L γ and R). The source is optically thick to photodisintegration of the primary nuclei, and the nuclear cascade can efficiently develop. The ejected cosmic rays receive a significant contribution from neutrons (disintegration products), which can easily escape because they are electrically neutral. While at the highest energies secondary nuclei produced in the disintegration chain dominate, cosmic rays from neutron decays dominate at lower energies, and the ejected overall spectrum becomes softer. The ejected UHECR composition is clearly mixed even if a pure composition is injected from the acceleration zone. The neutrino production becomes more efficient, and is dominated by photo-meson processes off secondary nuclei produced in the nuclear cascade. -Optically Thick Case (large L γ or low R). The source is optically thick to photohadronic interactions of all nucleons and nuclei. In this case, the nuclear cascade is so efficient that most energy ends up in nucleons (protons and neutrons), which are the end of the disintegration chain. The ejected cosmic ray spectrum is strongly dominated by neutrons, and it is very difficult to reach high cosmic ray energies because nuclei will mostly disintegrate. The neutrino flux is extremely high, dominated by photo-meson production off nucleons.
We have discussed how these scenarios quantitatively depend on the shell parameters, and how our results can be easily translated into different astrophysical models. We note that the ejected UHECR spectra do not exhibit a rigidity-dependent maximal energy cut-off if the nuclear cascade is populated, which is contrary to the assumption of many UHECR models in the literature.
Regarding the impact of different photo-meson models on the neutrino flux, we have pointed out that relatively solid predictions can be obtained if the photo-meson production has a significant contribution from nucleon-photon interactions (Populated Cascade and Optically Thick Case), for which the hadronic physics is relatively well known (such as using SOPHIA), whereas the photo-meson production off nuclei (Empty Cascade) requires further study. The expected neutrino flux in the latter case is, however, smaller.
As one of the most interesting implications for GRBs, it has been shown that the expected neutrino fluence, for the same injection luminosity, spectra, and shell parameters, hardly depends on the injection composition. A profound consequence from this observation is that, apart from different cosmic ray propagation effects, the neutrino stacking bounds on longduration GRBs will also apply if the UHECRs are, in fact, nuclei. The main difference between protons and nuclei comes from the differences in the maximal energies at which they are accelerated in the sources and from the interaction lengths of nuclei in the extragalactic space.
The next step was the extrapolation from one GRB shell to an entire population of long-duration GRBs with identical parameters for all zones in the internal shock scenario (often referred to as the "one-zone model"). The emitted UHECRs were propagated to Earth. It was found that the Auger spectrum and composition can be reasonably well described even with a pure injection composition only, were two cases have been discussed: -Mixed Composition Dip Model. The same class of sources is responsible for UHECR data including the ankle. Compared to the conventional proton dip model, where the dip comes from pair production, the same feature is generated by the nucleons from disintegration (in the sources and during propagation), which pile up at energies below the ankle. -Mixed Composition Ankle Model. The UHECR data needs to be described only above the ankle, whereas a different population is expected to dominate at lower energy.
It has been demonstrated that both cases can be reasonably well described within the GRB framework. However, the dip hypothesis puts narrower constraints on the parameters. The protons below the ankle, often generated within the source, enhance the (prompt) neutrino flux significantly. As a consequence, the Mixed Composition Dip Model is in tension with current neutrino data for the chosen set of parameters, while the Mixed Composition Ankle Model is consistent with current experimental observations, meaning that GRBs can be the sources of UHECR nuclei in spite of low predicted neutrino fluxes. While the results have been derived for the injection of 28 Si, the dependence on the injection composition has been illustrated with a novel approach connecting the physics of the source with the physics of propagation. There are several limitations to our fitting procedure. First of all, the overall goodness of fit is poor in the Mixed Composition Dip Model because of the small statistical error bars. How-ever, we have not yet included systematic errors, such as the energy scale uncertainty of the experiment. In addition, it is clear that a mixed injection composition in the source will improve the goodness of fit, which is, however, beyond the scope of this work. In addition, some parameters (such as the spectral index, the shape of the cut-off at the injection and the source evolution) have been fixed, whereas not doing so would also improve the fit. Therefore, the results in this work can be considered as "proof of principle" for the UHECR fit, and as "indicative" for the neutrino bounds.
We conclude that GRBs can be the sources of the UHECRs, as we have self-consistently described the observed spectrum and composition observed by Auger even in a one-zone model with a pure injection composition into the GRB shells. However, the exclusion power of neutrino stacking bounds applies to nuclei as well, which in particularly constrains the hypothesis that the light composition below the ankle comes from the same population. If, however, the UHECR transition to the GRB contribution occurs at the ankle, GRBs can describe cosmic ray and neutrino data at the expense of relatively high baryonic loadings and obtained GRB parameters (either low γ-ray luminosities or large collision radii).
While some of our results can be applied to multi-zone models (such as the behavior of the nuclear cascade), these models typically predict a lower neutrino flux Globus et al. 2015a;Bustamante et al. 2017). A self-consistent picture has been drawn in Globus et al. (2015a,b) for a particular set of parameters, such as a specific injection composition mix. Dedicated parameter space studies in this context will require more efficient computing techniques, especially if the nuclear cascade is to be computed in each collision independently. In this study, we have presented the technology to perform such computations efficiently and in a sufficiently precise way.
The function χ j→i , which depends on the kinematics of the process, describes which (mean) fraction of the parent energy is deposited in the secondary. If a secondary nucleus is produced, a frequently applied assumption is Lorentz factor conservation, which leads to χ j→i ≃ A i /A j .
Using Eq. (A.1) in Eq. (3), one easily finds for the injection of secondary nuclei from beta decays and spontaneous emission: For instance, for the case of β ± decays, A does not change, which means that χ j→i ≃ 1, and M j→i is the branching ratio into that channel.
Refined computations for the neutrino spectrum of beta decays from relativistic ions can be, for example, found in the context of long-baseline experiments, so-called β-beams (see Burguet-Castell et al. 2004). For the neutrino injection, we use Eq. (A.2) with a peak value χ j→ν directly from the peak of the re-distribution function extracted from Burguet-Castell et al. (2004), related to the Q-value (obtained from the mass difference between parent and daughter nucleus).

A.2. Photodisintegration
The interaction rate Γ ′ j for photohadronic interactions in an isotropic target photon field is given by Here n ′ γ (ε ′ ) is the photon density as a function of photon energy ε ′ and the pitch angle between the photon and proton momenta θ ′ Aγ , σ abs j (ǫ r ) is the absorption cross section for species j, and is the photon energy in the parent rest frame (PRF) in the limit β ′ A ≈ 1. We note that ǫ r corresponds to the available center of mass energy √ s of the interaction, as s = m 2 A + 2 m A ǫ r . It is convenient to write the interaction rate in the form of an integral over the photon density n ′ γ (ε ′ ) as where y ′ ≡ (E ′ j ε ′ )/m A corresponds to the "typical" center-ofmass energy, and is an integral over the cross section (which is, by definition, zero below the threshold); it can be interpreted at pitch-angle averaged cross section. Here the idea is that the function f j (y) can be pre-computed (it only depends on the cross section, and the integral takes care of the pitch angle averaging), and that the interaction rate can be obtained in a single integral in Eq. (A.5).
For the secondary nuclei injection, we re-write Eq.
(3) (using Eq. (A.1)) as with χ j→i ≡ A i /A j and a function Here it is taken into account that the secondary multiplicity strongly depends on the center-of-mass energy (or ǫ r ) with precomputed functions g j→i , while the re-injection can be still obtained by a single integral Eq. (A.7). For the case of a constant target photon spectrum (such as in this work), not even a single integral is needed.
The values for f (y ′ ) and g j→i (y ′ ) are produced from the disintegration models (such as TALYS), as described in the main text. We note that σ abs j (ǫ r ) M j→i = σ incl j→i (ǫ r ) (total inclusive cross section), which can be extracted from the nuclear models. Comparing Eq. (A.8) with (A.6), we can read off that the quantity g j→i (y ′ )/ f (y ′ ) describes the secondary multiplicity as a function of y ′ , including the pitch-angle averaging in the isotropic target photon field. We illustrate this quantity for three different examples in Fig. A.1 for y ′ = 50 MeV (which is slightly above the GDR, because the model is more sophisticated there). In all cases, nucleons are produced in each integration process, as well as light nuclei. The residual nuclei tend to be populated along the main diagonal. However, compared to one-dimensional (one isobar) disintegration models, such as the Puget-Stecker-Bredekamp model (Puget et al. 1976;Khan et al. 2005), the isotope chart will be populated in two dimensions: along the main diagonal, and perpendicular to it in terms of unstable isobars. The importance of using a sufficiently sophisticated disintegration model has been addressed in Boncioli et al. (2017).

A.3. Photo-meson production
The interaction rate for photo-meson production can be obtained as in Eq. (A.5). The main difference to disintegration, however, is that the re-distribution functions of the secondary pions should be taken into account. The approach presented here is a further advancement of Hümmer et al. (2010).
If the re-distribution function for the secondaries (especially the pions) is to be described by Eq. (4), it turns out to be difficult to avoid double integrals in the injection function, such A101, page 26 of 31 as over nucleus and photon energy. These, however, are one of the bottlenecks for efficient computations, and brute-force sampling methods of the interaction models frequently do not take this into account. The idea in Hümmer et al. (2010) for pγ interactions was basically to discretize one of these integrals into a small number of "interaction types" (in their language), where the splitting into t-channel production, resonances, and multipion production was physics-motivated. These interaction types had different characteristics in terms of their multiplicities and inelasticities, and were evaluated similarly to Eq. (A.7). We propose a different method here, which will allow for automatic definitions of the interaction types for many isotopes, and which will even outperform Hümmer et al. (2010) in terms of efficiency and precision while being based on similar principles.
We re-write one of the integrals into an x-integral, and we define x-dependent "interaction types" in the language of Hümmer et al. (2010) by discretizing that integral. We find that for T such interaction types, the injection of secondaries is given by wherex = log 10 (x), and where we define a new re-distribution function Here we identified σ abs j dn j→i /dx with the (differential) inclusive cross section dσ incl j→i /dx. The idea is similar to the photodisintegration: Eq. (A.9) only contains a single integral, to be summed over several discrete valuesx k . The information on h can be directly compiled from the inclusive cross sections once an appropriate splitting in terms of x is defined.
Comparing Eqs. (A.10) to (A.8) (and the corresponding reinjection functions), we note that g j→i (y ′ ) = T k=1 ∆x k h j→i (x k , y ′ ), (A.11) which implies that the re-distribution function for the secondaries has to add up to yield the secondary multiplicity, and the interaction types should be chosen accordingly. The simplest example, frequently used in the literature, is to pick T = 1 and dn j→i /dx(x 1 , y ′ ) = M j→i withx 1 = log 10 χ j→i , in which case Eq. (A.9) reduces to Eq. (A.7). It is clear that using several such interaction types with different values ofx will lead to more precise results, at the expense of computation time linearly growing with T . Compared to Hümmer et al. (2010), one can identify their interaction types M IT j→i f IT (y ′ ) ↔ ∆x k h j→i (x k , y ′ ) with ours (compare Eq. (A.9) to Eq. (28) in Hümmer et al. 2010). Our splitting into interaction types with correspondingx-values is performed as a function of the average interaction energy y ′ , which allows for energy-dependent multiplicities. Since the contributing physical processes change as a function of interaction energy, both approaches lead, in principle, to similar results.
Let us test this for the example of pγ interactions using SOPHIA (Mucke et al. 2000). First of all, we use the mean value theorem on Eq. (A.10), which yields h j→i (y ′ ) ≃ f j→i (y ′ ) dn j→i dx (x k , ǫ r ), (A.12) with an appropriate choice 0 ≤ ǫ r ≤ 2y ′ . It turns out that ǫ r ≃ y ′ is a good approximation (the "typical center-of-mass energy"). In that case the re-distribution function can be easily extracted from SOPHIA without modifications, assuming an isotropic target photon spectrum and appropriate values of ε ′ and E ′ p for a flat target photon spectrum (recall that y ′ = ε ′ E ′ p /m p ). We have tested several different splittings inx, and found that using T = 4 withx 1 = −1.4,x 2 = −1.0,x 3 = −0.6, andx 4 = −0.2 (for fixed ∆x = 0.4) gives reasonable results, while being computationally inexpensive. We illustrate the results for this interaction model and its comparison to SOPHIA in Fig. A.2. The upper left panel illustrates the distribution function dn/dx(x k , y ′ ) for different values ofx. The typical assumption for photohadronic interactions is that pions take about 20% of the proton energy, corresponding to x ≃ 0.2. From the figure, one can easily see that the re-distribution funcion for x = 0.25 peaks at around the ∆-resonance (a few hundred MeV). For higher energies, where multi-pion production dominates, smaller values of x yield larger pion multiplicities. An important result that the figure shows is captured by the y ′ -dependence of the redistribution function, which is different from the usually shown x-dependence (for fixed interaction energy): the functions describe the injection contributions from different (pitch-angle averaged) center-of-mass energies for fixed values ofx, which correspond to the ratio of daughter and parent energies.
The upper right panel of Fig. A.2 shows the contributions from the differentx-types (same color-coding), and the total (light blue). One can easily see that the main contribution comes fromx ≃ 0.25, while low (high) energies receive some contribution from smaller (larger)x-values. The total spectrum matches the one produced with SOPHIA (red dashed) very well and even better than Hummer et al. (2010, dotted, model Sim-B therein). The better relative reproduction of the spectrum is also documented in the lower left panel of Fig. A.2, whereas Hümmer et al. (2010) reproduce the charged pion ratios somewhat better at low energies (see lower right panel). In summary, we obtain a precision somewhat better than Hümmer et al. (2010), with only four instead of 23 interaction types, that is, the new method is more than a factor of five faster. We have also tested alternatives with eight and 17x-values, which, however, do not provide a significant gain in precision.
The implementation of the superposition model for nuclei is described in the main text, wherex A k =x p k − log 10 A can be easily found from thex-value for protons and neutrons. This model does not yet exploit the strength of our method in full, which is to be discussed in future publications.

Appendix B: Effect of minimum photon energy cut-off
In the main text of this study, we assume that the minimal photon energy is low enough such that UHECRs will always find an interaction partner to produce the giant dipole resonance. For UHECR with a γ-factor of 10 10 in the SRF, that means that ε ′ γ,min 10 −3 eV to produce the GDR. If, however, the low-energy photon spectrum cut-off is at higher energy, as it may come from synchrotron self-absorption (Wang et al. 2008;Murase et al. 2008), the disintegration at the highest energies will be dominated by the photo-meson regime. Let us choose one extreme example here, ε ′ γ,min = 0.1ε ′ γ,br ≃ 100 eV, such that photodisintegration will cease about one order of magnitude beyond the break.
A comparison between this case (thick curves) and our standard assumption (no such cut-off, thin curves) is shown in  , it is clear that the ceasing of the disintegration rate leads to slightly higher maximal energies, which are recovered in the upper right and lower left panels for the densities in the source and the ejected spectra, respectively. The effect of the maximal energy on the neutrino spectra (lower right panel) is small because the peak flux, coming from disintegration products (nucleons), is hardly affected. Therefore we expect slight modifications of the UHECR fit in this case, whereas the neutrino flux predictions hardly change. Since the UHECR output is actually somewhat higher, the tension with the prompt neutrino flux, which scales with the baryonic loading, will be somewhat released. approach, however, neglects the fact that the proton (steady) spectrum itself deforms in the case of high optical thicknesses, and that there is substantial feed-down of particles to lower energies by multiple interactions.
In this study, we treat the proton-neutron system as a time-dependent fully-coupled partial-differential-equation system, computing the steady state explicitly, as illustrated in Fig. 1  (left panel). The photohadronic energy losses are treated discretly, that is, the protons and neutrons escape with the photohadronic rate, and are partially re-injected at lower energies. In the optically thin (to photohadronic interactions) regime, this method yields results similar to earlier works. Compare the solid and dashed curves in the upper panels of Fig. C.1, which shows a comparison between the escape spectra computed with the method in Baerwald et al. (2013, dashed curves) and this study (thick solid curves). In these cases, the optical thickness to photohadronic interactions is lower than one, which means that there is hardly any feedback from the neutrons to the protons, and the proton spectrum is hardly affected by photohadronic interactions. There are small differences, which come from the slightly different effective maximal energy (see curves "p steady state", showing N ′ p /t ′ dyn , compared to "initial proton" curves). The reason is that we add an adiabatic cooling term for the protons in this study, which produces exactly the same result as an equivalent escape term for a power-law spectrum (cf., initial proton and steady proton spectra), but modifies the shape of the cut-off (see Eq. (1): a cooling term is sensitive to the derivative of N ′ p at the cut-off).
There are, however, differences in the optically thick regime (see Fig. C.1 for two high luminosity cases (lower panels)). One important difference is that we normalize the injection luminosity with the baryonic loading (see Eq. (11)), as compared to the steady state density in Baerwald et al. (2013, which needs to be powered by a correspondingly higher injection luminosity). The proton steady state is then computed explicitly, and, especially at the highest energies, suppressed by the photohadronic interaction rate. For L = 10 52 erg s −1 , the difference between the old and new methods is moderate because the optical thickness is moderate (around three). The new method implies a slightly lower normalization (because of the injection luminosity normalization) and a slightly stronger depletion of the neutron escape spectrum at higher energies (because of multiple interactions).
For extremely large luminosities L = 10 53 erg s −1 (optical thickness: 36), the protons (and neutrons) from multiple interactions in fact pile up close to the pion production threshold, whereas the escape spectra at the highest energies are heavily suppressed, which leads to softer escape spectra. The neutron density N ′ n follows the proton density N ′ p , because the protonneutron system becomes fully coupled by forward and backward reactions in the case of high optical thicknesses. This means that there is an impact on the ejected cosmic ray flux shape in the optically thick regime, whereas one can easily see from the figure that the neutrino flux is less affected. In that case, the spectral neutrino peak comes from the photo-pion production close to the threshold, where the model-dependent impact is smaller. While the change of normalization reduces the flux somewhat, the additional production channel from neutron interactions, which we take into account here, partially compensates for that.