EDP Sciences
Open Access
Issue
A&A
Volume 611, March 2018
Article Number A101
Number of page(s) 31
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201731337
Published online 18 April 2018

© ESO 2018

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

The origin of ultra-high-energy cosmic rays (UHECRs) at the highest energies ≳1018 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 high-energy 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 gamma-ray 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 (Hümmer et al. 2012; 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 (Baerwald et al. 2015). 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 (Bustamante et al. 2015; 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 (Bustamante et al. 2015). 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 model – especially 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. (2015a,b, 2017) for UHECR nuclei. A GRB proton dip model clearly overproduces the prompt neutrinos (Baerwald et al. 2015), 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.

2 Methods

Our simulations are based on the NeuCosmA code, following the implementation for GRBs in Hümmer et al. (2012); Baerwald et al. (2012, 2013, 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.

2.1 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: (1)

where (with the energy loss rate ), and is the escape rate. The PDE system is to be solved for the differential particle densities in the shock rest frame (SRF); primed quantities refer to the SRF. The coupled system arises because of the injection term (2)

which allows for injection from an acceleration zone (typically a power law), as well as for injection from other species j with the term , such as from photodisintegration or β± decays. We will discuss the relevant interaction processes below. Fully-stripped ions are identified by AZ, where A is the mass number, Z is the charge (or number of protons), and N = AZ is the number of neutrons1.

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 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 acts on species i (and is attached to that PDE), but involves another species’ density , 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. (2012, 2013, 2015), which are effectively performed in the optically thin regime to photohadronic (, ) 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 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 re-inject 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 interactions2 . Although this simple example is only meant for illustration, we note that the optically thick case of photohadronic interactions for proton injection is included in our work; Appendix C gives a more detailed discussion and a comparison to earlier works (accessible after a first read through this section).

Case B)in Fig. 1 refers to nuclei, where only one PDE for one species is illustrated. Similarly, the fully-stripped nuclear isotope will have cooling processes attached to it (synchrotron cooling, adiabatic cooling, and pair production cooling), as well as escape processes (photodisintegration, photo-meson production, and decay, if applicable, which change the species). There may be injection from the acceleration zone, where we only inject a pure composition with one isotope at a time in this study for the sake of simplicity. The motivation for that is two-fold: first of all, it is clear that injecting an isotope cocktail will lead to better fit results. It is therefore interesting to see how far the simplest possible model can describe data. Second, we would like to develop an understanding of the impact of the injection composition, which is difficult if a cocktail of isotopes is injected from the beginning. As injection isotopes, we choose the most abundant stable isotope (A Z) for each element (Z), whenever the results are shown as a function of Z (dark blue isotopes in Fig. 2). The injection from different species can be photohadronic processes (photodisintegration or photo-meson production, depending on the center-of-mass energy), beta decays, or spontaneous decays of proton- or neutron-rich isotopes, as illustrated in Fig. 1. Escape processes, as we define them, not only include escape from the region or escape over the dynamical timescale, but also conversion into a different species (such as by decay or interactions).

thumbnail Fig. 1

Schematic illustration of the partial differential equations (PDEs) for: A) an optically thick coupled proton-neutron system with proton injection (“inject”, from the acceleration zone), and B) an isotope in the nuclear cascade, possibly with injection (from the acceleration zone). Each dashed box corresponds to one species (or one PDE), to which the relevant cooling, escape, and injection processes are “attached”. We note that all disintegration, photo-meson, or decay channels destroying a species (p, n, or A Z) show up as escape terms in the PDE of that species.

Open with DEXTER
thumbnail Fig. 2

Considered 481 isotopes in this work as a function of neutron number N and proton number Z. Here it is assumed that the heaviest stable injection isotope is 56Fe. The color-coding refers to the most abundant stable isotopes (for each element, i.e., equal Z), which are considered as injection isotopes (dark blue), stable isotopes (blue), β± emitters (dark red), possibly followed by the spontaneous emission of nucleons (light red), and spontaneous emitters of nucleons or α particles (white). We note that only the leading (largest branching ratio) processes are shown, and that isotopes with lifetimes longer than one month are regarded as stable in our scheme. Isotopes with rest frame lifetimes τ0 ≤ 10−5 s are marked by “×”, which are potentially interesting for neutrino production in GRBs if they are β± emitters. Isotopes with rest frame lifetimes τ0 ≤ 10−10 s (or unknown lifetimes) are marked by dots, which decay quickly enough even at the highest energies to be integrated out.

Open with DEXTER

2.2 Hadronic radiation processes

The injection rate in Eq. (2) of particles of species i and energy Ei from the interaction or decay of the parent particle j can be approximated by (3)

where is the interaction rate, is the fraction of primary energy taken by the secondary, and is the center-of-mass energy (if applicable). In fact, for photo-meson production, one has to integrate over d nji∕dx within the photon energy integral, which is “hidden” in the interaction rate here; see Appendix A. The function (4)

describes the distribution of secondaries of type i per final state energy interval dx, the probability distribution is normalized ∫ pji(x)dx = 1, and the integral over that function gives the average number of secondaries Mji 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.

2.2.1 Beta decays and spontaneous emission

As a starting point, we choose all potentially relevant known isotopes with A < 56 and 56Fe, 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 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 , which corresponds to neutrino energies ≲104 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. 2011, 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 at the highest energies (where ). 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.

2.2.2 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’ restframe). 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 . 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 NA, and ZA if a proton interacts, respectively. The inelastic cross section is assumed to scale like σ = . 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” A2∕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 4He (MacCormick et al. 1997). At higher energies the photon interacts in a more hadron-like way, where the A2∕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 photo-meson 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 photo-meson 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.

2.2.3 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 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 ). 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.

2.3 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 ai1 →… → iNb, where the isotopes i1iN 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 and 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 photo-meson production, where the interaction types (-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 56Fe), 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 56Fe 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.

thumbnail Fig. 3

Automatic isotope selection result as input for the PDE system setup for the injection of 56Fe, tracing recursively all possible disintegration, decay, and mixed paths. The isotopes are selected if their secondary multiplicity M (identified with the pitch-angle averaged number of secondary nuclei Mgji(y)∕fj(y), see Appendix A) exceeds a certain threshold, as indicated by the legend, for any 8 MeV ≲ y ≲ 125 MeV. The threshold M > 0.01 has been used for the computations in this work.

Open with DEXTER

2.4 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 (5)

where is a normalization factor and (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 (6)

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 ΔTtv 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. (2015, 2017).

The normalization of the photon spectrum in Eq. (5) is obtained from (7)

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) often used for analytical computations (assuming that tv 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 in Eq. (7).

2.5 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 (8) and the variability timescale is assumed to be indicative for the shell width (9)

    All collisions are assumed to collide at the same radius in the one-zone model. As a consequence, Eq. (7) strongly depends on Γ and tv , and the pion production efficiency becomes fLγ/(Γ4tvϵγ,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), Rm collision radius, (width of merged shell, related to intermittent timescale of the central emitter), and (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 one-zone 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 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 Rph, 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 (indicative for the shell width in Eq. (9)) and a slow one (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 , which means that it is suppressed by the factor compared to the standard internal shock model.

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 tv = 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.

2.6 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: (10)

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 is determined automatically by balancing the acceleration rate (where 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 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 (11)

to determine in Eq. (10), with the nuclear loading ξi, and 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.

2.7 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, , νμ , ), 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 ) 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 tdyn) the steady state is reached. For the integration, we re-parameterize the PDE system in terms of 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 detail3. 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 escapes over the dynamical timescale. We note that fesc = 1 if the Larmor radius reaches (or exceeds) the shell width, which corresponds to the free-streaming limit. It can be shown that fesc is invariant even if the shell expands (Baerwald et al. 2013). The escape fraction can be translated into an effective escape rate (for ). It is interesting to compare that rate to the escape from the whole shell: in this case , where D is the (spatial) diffusion coefficient. This assumption yields the same result if (Bohm diffusion). 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 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.

2.8 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 low-energy 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 log10(Einj∕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).

2.9 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.

3 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 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 56Fe only, which is the heaviest commonly discussed injection isotope. We will show continuous parameter space scans in the following sections.

3.1 Empty cascade

We define the Empty Cascade as the case that is optically thin to 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 are shown in the upper left panel of Fig. 4 for one example. By comparing and with (see arrows), one can clearly see that the source is optically thin to 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 56Fe is indeed hardly modified and extends to the mentioned maximal energy. Secondary 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 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 56Fe 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 4He) 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 Fig. C.1 for the magnitude of the effect).

thumbnail Fig. 4

Example for the “Empty Cascade” source class (isotropic luminosity Lγ = 1049 erg s−1 ) for injection of pure 56Fe. 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 56Fe 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 ≃ 108.3 km, Γ = 300, ξFe = 10, keV, and z = 2.

Open with DEXTER

3.2 Populated cascade

We define the Populated Cascade as the case that is optically thick to interactions of the injection isotope that will disintegrate and populate the cascade. At the same time, the source environment is optically thin to 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 interactions. The maximal primary energy is given by 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.

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 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.

thumbnail Fig. 5

Example for the “Populated Cascade” source class (isotropic luminosity Lγ = 1051 erg s−1 ) for injection of pure 56Fe. The caption of Fig. 4 gives details.

Open with DEXTER

3.3 Optically thick case

We define the Optically Thick Case as the one that is optically thick to 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 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 interactions limiting the maximal primary energy.

thumbnail Fig. 6

Example for the “Optically Thick Case” (to nucleons and nuclei) source class (isotropic luminosity Lγ = 1053 erg s−1 ) for injection of pure 56Fe. The caption of Fig. 4 gives details.

Open with DEXTER

4 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 scattering4 .

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 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 interaction rates increase with 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 interaction rate increasingly approaches the 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 tv = 0.01 s fixed (instead of scaling it with the internal shock model relationship Eq. (8) among collision radius, tv , and Γ). We note that the pion production efficiency fLγtvR2 (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 is proportional to the (constant) adiabatic cooling rate, one has for fixed maximal energy. This observation is interesting because the cosmic ray fit is sensitive to the maximal energy: since maximal energy and 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 tv). 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 tv 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 = 1049 erg s−1, R = 108.3 km internal shock model prototype will correspond to a L = 1051 erg s−1, R = 109.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.

thumbnail Fig. 7

Nuclear cascade class as a function of luminosity Lγ and collision radius R in the internal shock scenario. The injection composition is pure 16O (left left) or 56Fe (right panel). The other parameters are the same as in Fig. 4. Black dots indicate the position of the examples shown in Sect. 3. The gray dashed contours display log 10 (Ei,max[GeV]) in the observer’s frame. Below the photosphere (red solid line), gamma-rays cannot leave the source anymore due to electron-positron pair production. We note that in this figure we fix Γ = 300, and scale tvR according to Eq. (8).

Open with DEXTER
thumbnail Fig. 8

Impact of different astrophysical model assumptions (see marked shadings) for fixed tv = 0.01 s. Similar to Fig. 7 for 56Fe, but the internal shock model relationship Eq. (8) among collision radius, tv , 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 tv determines the shell thickness and Γ the energy shift only, that is, the dependence on these parameters is milder (the pion production efficiency scales ∝ tv 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 (taken from Zhang & Kumar 2013, for multi-collision model see Sect. 2.5) are sketched in the figure for the chosen value of tv (fast time variability in light curve). We note that the maximal energies are shown as gray dashed curves as in Fig. 7, but not the transition curves among the cascade regions. Instead, the dashed-dotted line separates the dominated (lower right corner) and adiabatic-loss-dominated (upper left corner) maximal primary energy regions.

Open with DEXTER

5 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 56Fe, 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 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 56Fe, 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 photo-meson 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 σ. 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 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.

thumbnail Fig. 9

Contribution of primary and secondary isotopes to neutrino production for the injection of pure 56Fe. The panels show the total all-flavor neutrino fluence per shell of a GRB as a function of the energy in the observer’s frame. Different luminosities (panels) correspond to the examples for the prototypes in Sect. 3. The different curves show the contribution of the photo-meson production off the primary (56 Fe) and secondary isotopes produced in the photodisintegration, where the contribution from protons and neutrons is separated.

Open with DEXTER
thumbnail Fig. 10

Impact of injection composition: total all-flavor neutrino fluence per shell of a GRB as a function of the energy in the observer’s frame. Different luminosities (panels) correspond to the examples for the prototypes in Sect. 3. The different curves correspond to the injection of pure isotopes: protons, 56Fe, or any other isotope in between (see legend). This comparison relies on equal injection luminosities for the primaries with a baryonic loading (luminosity in nuclei versus photons) of ten, and on equal injection spectra ∝ E−2 .

Open with DEXTER

6 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 ΔTtv 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 cross-section 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 1018 eV (Mixed Composition Dip Model) and above 1019 eV (Mixed Composition Ankle Model) by using the combined spectrum (Valiño et al. 2015) and the shower depth (Xmax) 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 Matteo et al. (2015) 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.

6.1 Mixed Composition Dip model

We first fit the spectrum and composition above 1018 eV, including the ankle region. The result of the fit is shown in Fig. 11, where the as a function of R and Lγ is given. The blue region indicates the current IceCube exclusion region from 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 × 104. 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 the corresponding prompt and cosmogenic neutrino fluxes are shown. The best-fit case (A) is excluded by the 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.

thumbnail Fig. 11

Parameter space scan for Mixed Composition Dip Model in the internal shock scenario. Here 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 1018 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. (2015, 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 ( log10ξ).

Open with DEXTER
thumbnail Fig. 12

Cosmic ray observables obtained with the best-fit parameters in the Mixed Composition Dip Model, corresponding to point A for (log10(R∕km), log10(Lγ∕(erg s−1))) = (8.1,50.5) in Fig. 11. Top: simulated energy spectrum of UHECR (multiplied by E3 ) 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 Xmax 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), 14N (green), and 56Fe (blue)) compositions, compared to data from Porcelli et al. (2015).

Open with DEXTER
thumbnail Fig. 13

Cosmic ray spectra and composition observables obtained for selected points in the parameter space in Fig. 11 of the Mixed Composition Dip model (similar to Fig. 12). Here the points are: A: (log 10 (R∕km), log10(Lγ∕(erg s−1))) = (8.1, 50.5), B: (8.1, 49.5), C: (8.1, 51.7).

Open with DEXTER
thumbnail Fig. 14

Prompt and cosmogenic (muon flavor) neutrino spectra obtained for selected points in the parameter space in 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: (log10(R∕km), log10(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.

Open with DEXTER

6.2 Mixed Composition Ankle Model

The fit is performed above 1019 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 σ(Xmax) is almost flat. The same message is given by the corresponding cosmogenic flux shown in the right panel of Fig. 18.

thumbnail Fig. 15

Parameter space scan for Mixed Composition Ankle Model in the internal shock scenario. Here 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 1019 ev, including a penalty for the overshooting of the flux at lower energies. The caption of Fig. 11 gives details.

Open with DEXTER
thumbnail Fig. 16

Cosmic ray observables obtained with the best-fit parameters in the Mixed Composition Ankle Model, corresponding to point A (log10(R∕km), log10(Lγ∕(erg s−1))) = (8.1, 49) in Fig. 15. The caption of Fig. 12 gives details.

Open with DEXTER
thumbnail Fig. 17

Cosmic ray spectra and composition observables obtained for selected points in the parameter space in Fig. 15 of the Mixed Composition Dip Model (similar to Fig. 16). Here the points are: A: (log 10 (R∕km), log10(Lγ∕(erg s−1))) = (8.1, 49), B: (7.8, 50), and C: (9.6, 51).

Open with DEXTER
thumbnail Fig. 18

Prompt and cosmogenic (muon flavor) neutrino spectra obtained for selected points in the parameter space in Fig. 15 for the Mixed Composition Ankle Model. Here A: (log 10 (R∕km), log10(Lγ∕(erg s−1))) = (8.1, 49), B: (7.8, 50), and C: (9.6, 51).

Open with DEXTER

7 Impact of injection composition

So far, we have tested several selected injection isotopes: 16O and 56 Fe to describe the disintegration in a GRB shell as kind of extreme examples for nuclei, and 28Si 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 28Si 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 Emax. This maximal energy is determined by the equilibrium between acceleration rate and the sum of energy loss rates (see Sect. 2.6).

  • Composition at Emax. 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 Etrans has to be 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 1049 erg s−1Lγ ≲ 1051 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 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 1049 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 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. From the figure, one can see that injection charges between about 7 (nitrogen) and 14 (silicon) may provide reasonable results. For example, for the chosen Z = 14 in the previous section, the other estimators are 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 1049 erg s−1Lγ ≲ 1051 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.

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.

thumbnail 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 Emax and the transition energy Etrans. The maximum energy Emax is determined by acceleration and loss processes.

Open with DEXTER
thumbnail Fig. 20

Qualitative estimators, as defined in the main text, for the UHECR fit as a function of injection Z. The shaded regions illustrate our expectations from cosmic ray data (top left from Heinze et al. 2016; top right and bottom left from Aab et al. 2017; and bottom right from Kampert & Unger 2012). The different curves correspond to different values of Lγ , as illustrated in the legend. In the upper right panel, the maximal energy following the Peters cycle (rigidity dependent maximal energy) is illustrated for two examples. The production radius is fixed to R ≃ 108.3 km in this figure.

Open with DEXTER
thumbnail Fig. 21

Ejected average composition of the cosmic rays leaving the GRB shell as a function of the energy in the observer’s frame. Different luminosities correspond to the examples for different source classes in Sect. 3. Different curves correspond to different pure injection compositions.

Open with DEXTER

8 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 long-duration 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 28Si, 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. However, 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 (Bustamante et al. 2015; 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.

Acknowledgements

We would like to thank Mauricio Bustamante and our colleagues from the Pierre Auger Collaboration for fruitful discussions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant No. 646623).

Appendix A Efficient computation of nuclear processes

This appendix contains details about the computational methods used for the efficient treatment of the interaction (and decay) processes, balancing performance and precision. In this study, we assume that the target photons are isotropically distributed in the SRF.

A.1 Beta decays and spontaneous emission

The simplest approximation for the re-distribution function p in Eq. (4) is the δ-function, (A.1)

The function χji, 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 χjiAiAj.

Using Eq. (A.1) in Eq. (3), one easily finds for the injection of secondary nuclei from beta decays and spontaneous emission: (A.2)

here and χji = AiAj. For instance, for the case of β± decays, A does not change, which means that χji ≃ 1, and Mji 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 for photohadronic interactions in an isotropic target photon field is given by (A.3)

Here is the photon density as a function of photon energy ε and the pitch angle between the photon and proton momenta , is the absorption cross section for species j, and (A.4)

is the photon energy in the parent rest frame (PRF) in the limit . We note that ϵr corresponds to the available center of mass energy of the interaction, as .

It is convenient to write the interaction rate in the form of an integral over the photon density as (A.5)

where corresponds to the “typical” center-of-mass energy, and (A.6)

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 fj (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 (A.7)

with χjiAiAj and a function (A.8)

Here it is taken into account that the secondary multiplicity strongly depends on the center-of-mass energy (or ϵr ) with pre-computed functions gji, 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 gji(y) are produced from the disintegration models (such as TALYS), as described in the main text. We note that (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 gji(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).

thumbnail Fig. A.1

Secondary multiplicities for the disintegration of 16O (left), 28 Si (middle), and 56Fe (right), marked in blue, for the TALYS disintegration model used in this work. The color-coding refers to the pitch-angle averaged number of secondary nuclei gji (y )∕fj(y) produced on average for y = 50 MeV. The total (average) number of secondary nuclei and nucleons produced per interaction is shown in the lower right corner of each panel.

Open with DEXTER

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 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 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 multi-pion 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 (A.9)

where , and where we define a new re-distribution function (A.10)

Here we identified with the (differential)inclusive cross section d. The idea is similar to the photodisintegration: Eq. (A.9) only contains a single integral, to be summed over several discrete values . 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 re-injection functions), we note that (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 with , in which case Eq. (A.9) reduces to Eq. (A.7). It is clear that using several such interaction types with different values of 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 with ours (compare Eq. (A.9) to Eq. (28) in Hümmer et al. 2010). Our splitting into interaction types with corresponding -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 interactions using SOPHIA (Mucke et al. 2000). First of all, we use the mean value theorem on Eq. (A.10), which yields (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 for a flat target photon spectrum (recall that ).

We have tested several different splittings in , and found that using T = 4 with , , , and (for fixed ) 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 for different values of . 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 function 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 re-distribution function, which is different from the usually shown -dependence (for fixed interaction energy): the functions describe the injection contributions from different (pitch-angle averaged) center-of-mass energies for fixed values of which correspond to the ratio of daughter and parent energies.

thumbnail Fig. A.2

Illustration of new interaction model implementation (based on SOPHIA) for photo-meson production of π+ from interactions for GRB benchmark from Hümmer et al. (2010): re-distribution functions (upper left panel), spectral contributions (upper right panel), relative deviation from SOPHIA (lower left panel), and charged pion ratio (lower right panel). (See main text for details.)

Open with DEXTER

The upper right panel of Fig. A.2 shows the contributions from the different -types (same color-coding), and the total (light blue). One can easily see that the main contribution comes from , while low (high) energies receive some contribution from smaller (larger) -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 17 -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, where can be easily found from the -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 1010 in the SRF, that means that 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, , 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 Fig. B.1. From the rate plot (upper left panel), 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.

thumbnail Fig. B.1

Effect of a minimal photon energy cut-off at on the example in the Optically Thick Case shown in Figs. 6 and 9. Thin curves show the original plots, and thick curves the result with cut-off.

Open with DEXTER

Appendix C Optically thick proton-neutron system

In Fig. C.1 we discuss the impact of optical thickness on photohadronic interactions in a coupled proton-neutron system (with proton injection only). An effective treatment of the optically thick case was performed in Baerwald et al. (2013), where the proton spectrum was assumed to be a power law with a cut-off at the maximal energy given by the dominant radiation process (see thick black curves in Fig. C.1). The neutron injection was then computed using the proton and photon densities. Since the neutrons are stopped by photohadronic interactions over the respective mean free path, it was assumed that only neutrons produced at the edges of the shells (within the mean free path of the photohadronic interactions) can escape (see dashed blue curves). This is numerically equivalent to computing the steady state for the neutrons from , assuming that it is dominated by the photohadronic interaction rate at the peak of the spectrum. Then one uses , that is, neutrons from the whole region can escape over the free-streaming timescale we use this in this work. In that case, one has , which reproduces the assumption in Baerwald et al. (2013). This 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.

thumbnail Fig. C.1

Comparison of optically thin (upper row) and thick (lower row) proton-neutron systems for pure proton injection. Ejected particle fluences per shell are shown as a function of the energy in the observer’s frame for different luminosities. The dashed curves refer to the model used in Baerwald et al. (2013), while the solid lines refer to the model used in this work. (See main text for details.)

Open with DEXTER

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 discretely, 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 , 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 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 = 1052 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 = 1053 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 follows the proton density , because the proton-neutron 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.

References


1

We use the same symbol for the neutron number and differential number densities, as these cannot be mixed up and it is conventional notation.

2

The neutron decay lifetime is relatively long compared to the dynamical timescale relevant for GRBs.

3

These assumptions for particle escape are similar to Globus et al. (2015a), where this escape mechanism is called “high-pass filter”.

4

Assuming that the plasma is electron-(fully stripped) isotope-dominated 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 . 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 decreases the photospheric radius by .

5

Since the interaction rates (cf., Eq. (7)) and , one finds the maximal primary energy to be constant along the curves in the 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 .

All Figures

thumbnail Fig. 1

Schematic illustration of the partial differential equations (PDEs) for: A) an optically thick coupled proton-neutron system with proton injection (“inject”, from the acceleration zone), and B) an isotope in the nuclear cascade, possibly with injection (from the acceleration zone). Each dashed box corresponds to one species (or one PDE), to which the relevant cooling, escape, and injection processes are “attached”. We note that all disintegration, photo-meson, or decay channels destroying a species (p, n, or A Z) show up as escape terms in the PDE of that species.

Open with DEXTER
In the text
thumbnail Fig. 2

Considered 481 isotopes in this work as a function of neutron number N and proton number Z. Here it is assumed that the heaviest stable injection isotope is 56Fe. The color-coding refers to the most abundant stable isotopes (for each element, i.e., equal Z), which are considered as injection isotopes (dark blue), stable isotopes (blue), β± emitters (dark red), possibly followed by the spontaneous emission of nucleons (light red), and spontaneous emitters of nucleons or α particles (white). We note that only the leading (largest branching ratio) processes are shown, and that isotopes with lifetimes longer than one month are regarded as stable in our scheme. Isotopes with rest frame lifetimes τ0 ≤ 10−5 s are marked by “×”, which are potentially interesting for neutrino production in GRBs if they are β± emitters. Isotopes with rest frame lifetimes τ0 ≤ 10−10 s (or unknown lifetimes) are marked by dots, which decay quickly enough even at the highest energies to be integrated out.

Open with DEXTER
In the text
thumbnail Fig. 3

Automatic isotope selection result as input for the PDE system setup for the injection of 56Fe, tracing recursively all possible disintegration, decay, and mixed paths. The isotopes are selected if their secondary multiplicity M (identified with the pitch-angle averaged number of secondary nuclei Mgji(y)∕fj(y), see Appendix A) exceeds a certain threshold, as indicated by the legend, for any 8 MeV ≲ y ≲ 125 MeV. The threshold M > 0.01 has been used for the computations in this work.

Open with DEXTER
In the text
thumbnail Fig. 4

Example for the “Empty Cascade” source class (isotropic luminosity Lγ = 1049 erg s−1 ) for injection of pure 56Fe. 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 56Fe 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 ≃ 108.3 km, Γ = 300, ξFe = 10, keV, and z = 2.

Open with DEXTER
In the text
thumbnail Fig. 5

Example for the “Populated Cascade” source class (isotropic luminosity Lγ = 1051 erg s−1 ) for injection of pure 56Fe. The caption of Fig. 4 gives details.

Open with DEXTER
In the text
thumbnail Fig. 6

Example for the “Optically Thick Case” (to nucleons and nuclei) source class (isotropic luminosity Lγ = 1053 erg s−1 ) for injection of pure 56Fe. The caption of Fig. 4 gives details.

Open with DEXTER
In the text
thumbnail Fig. 7

Nuclear cascade class as a function of luminosity Lγ and collision radius R in the internal shock scenario. The injection composition is pure 16O (left left) or 56Fe (right panel). The other parameters are the same as in Fig. 4. Black dots indicate the position of the examples shown in Sect. 3. The gray dashed contours display log 10 (Ei,max[GeV]) in the observer’s frame. Below the photosphere (red solid line), gamma-rays cannot leave the source anymore due to electron-positron pair production. We note that in this figure we fix Γ = 300, and scale tvR according to Eq. (8).

Open with DEXTER
In the text
thumbnail Fig. 8

Impact of different astrophysical model assumptions (see marked shadings) for fixed tv = 0.01 s. Similar to Fig. 7 for 56Fe, but the internal shock model relationship Eq. (8) among collision radius, tv , 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 tv determines the shell thickness and Γ the energy shift only, that is, the dependence on these parameters is milder (the pion production efficiency scales ∝ tv 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 (taken from Zhang & Kumar 2013, for multi-collision model see Sect. 2.5) are sketched in the figure for the chosen value of tv (fast time variability in light curve). We note that the maximal energies are shown as gray dashed curves as in Fig. 7, but not the transition curves among the cascade regions. Instead, the dashed-dotted line separates the dominated (lower right corner) and adiabatic-loss-dominated (upper left corner) maximal primary energy regions.

Open with DEXTER
In the text
thumbnail Fig. 9

Contribution of primary and secondary isotopes to neutrino production for the injection of pure 56Fe. The panels show the total all-flavor neutrino fluence per shell of a GRB as a function of the energy in the observer’s frame. Different luminosities (panels) correspond to the examples for the prototypes in Sect. 3. The different curves show the contribution of the photo-meson production off the primary (56 Fe) and secondary isotopes produced in the photodisintegration, where the contribution from protons and neutrons is separated.

Open with DEXTER
In the text
thumbnail Fig. 10

Impact of injection composition: total all-flavor neutrino fluence per shell of a GRB as a function of the energy in the observer’s frame. Different luminosities (panels) correspond to the examples for the prototypes in Sect. 3. The different curves correspond to the injection of pure isotopes: protons, 56Fe, or any other isotope in between (see legend). This comparison relies on equal injection luminosities for the primaries with a baryonic loading (luminosity in nuclei versus photons) of ten, and on equal injection spectra ∝ E−2 .

Open with DEXTER
In the text
thumbnail Fig. 11

Parameter space scan for Mixed Composition Dip Model in the internal shock scenario. Here 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 1018 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. (2015, 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 ( log10ξ).

Open with DEXTER
In the text
thumbnail Fig. 12

Cosmic ray observables obtained with the best-fit parameters in the Mixed Composition Dip Model, corresponding to point A for (log10(R∕km), log10(Lγ∕(erg s−1))) = (8.1,50.5) in Fig. 11. Top: simulated energy spectrum of UHECR (multiplied by E3 ) 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 Xmax 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), 14N (green), and 56Fe (blue)) compositions, compared to data from Porcelli et al. (2015).

Open with DEXTER
In the text
thumbnail Fig. 13

Cosmic ray spectra and composition observables obtained for selected points in the parameter space in Fig. 11 of the Mixed Composition Dip model (similar to Fig. 12). Here the points are: A: (log 10 (R∕km), log10(Lγ∕(erg s−1))) = (8.1, 50.5), B: (8.1, 49.5), C: (8.1, 51.7).

Open with DEXTER
In the text
thumbnail Fig. 14

Prompt and cosmogenic (muon flavor) neutrino spectra obtained for selected points in the parameter space in 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: (log10(R∕km), log10(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.

Open with DEXTER
In the text
thumbnail Fig. 15

Parameter space scan for Mixed Composition Ankle Model in the internal shock scenario. Here 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 1019 ev, including a penalty for the overshooting of the flux at lower energies. The caption of Fig. 11 gives details.

Open with DEXTER
In the text
thumbnail Fig. 16

Cosmic ray observables obtained with the best-fit parameters in the Mixed Composition Ankle Model, corresponding to point A (log10(R∕km), log10(Lγ∕(erg s−1))) = (8.1, 49) in Fig. 15. The caption of Fig. 12 gives details.

Open with DEXTER
In the text
thumbnail Fig. 17

Cosmic ray spectra and composition observables obtained for selected points in the parameter space in Fig. 15 of the Mixed Composition Dip Model (similar to Fig. 16). Here the points are: A: (log 10 (R∕km), log10(Lγ∕(erg s−1))) = (8.1, 49), B: (7.8, 50), and C: (9.6, 51).

Open with DEXTER
In the text
thumbnail Fig. 18

Prompt and cosmogenic (muon flavor) neutrino spectra obtained for selected points in the parameter space in Fig. 15 for the Mixed Composition Ankle Model. Here A: (log 10 (R∕km), log10(Lγ∕(erg s−1))) = (8.1, 49), B: (7.8, 50), and C: (9.6, 51).

Open with DEXTER
In the text
thumbnail 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 Emax and the transition energy Etrans. The maximum energy Emax is determined by acceleration and loss processes.

Open with DEXTER
In the text
thumbnail Fig. 20

Qualitative estimators, as defined in the main text, for the UHECR fit as a function of injection Z. The shaded regions illustrate our expectations from cosmic ray data (top left from Heinze et al. 2016; top right and bottom left from Aab et al. 2017; and bottom right from Kampert & Unger 2012). The different curves correspond to different values of Lγ , as illustrated in the legend. In the upper right panel, the maximal energy following the Peters cycle (rigidity dependent maximal energy) is illustrated for two examples. The production radius is fixed to R ≃ 108.3 km in this figure.

Open with DEXTER
In the text
thumbnail Fig. 21

Ejected average composition of the cosmic rays leaving the GRB shell as a function of the energy in the observer’s frame. Different luminosities correspond to the examples for different source classes in Sect. 3. Different curves correspond to different pure injection compositions.

Open with DEXTER
In the text
thumbnail Fig. A.1

Secondary multiplicities for the disintegration of 16O (left), 28 Si (middle), and 56Fe (right), marked in blue, for the TALYS disintegration model used in this work. The color-coding refers to the pitch-angle averaged number of secondary nuclei gji (y )∕fj(y) produced on average for y = 50 MeV. The total (average) number of secondary nuclei and nucleons produced per interaction is shown in the lower right corner of each panel.

Open with DEXTER
In the text
thumbnail Fig. A.2

Illustration of new interaction model implementation (based on SOPHIA) for photo-meson production of π+ from interactions for GRB benchmark from Hümmer et al. (2010): re-distribution functions (upper left panel), spectral contributions (upper right panel), relative deviation from SOPHIA (lower left panel), and charged pion ratio (lower right panel). (See main text for details.)

Open with DEXTER
In the text
thumbnail Fig. B.1

Effect of a minimal photon energy cut-off at on the example in the Optically Thick Case shown in Figs. 6 and 9. Thin curves show the original plots, and thick curves the result with cut-off.

Open with DEXTER
In the text
thumbnail Fig. C.1

Comparison of optically thin (upper row) and thick (lower row) proton-neutron systems for pure proton injection. Ejected particle fluences per shell are shown as a function of the energy in the observer’s frame for different luminosities. The dashed curves refer to the model used in Baerwald et al. (2013), while the solid lines refer to the model used in this work. (See main text for details.)

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.