Impact of the Collision Model on the Multi-messenger Emission from Gamma-Ray Burst Internal Shocks

We discuss the production of multiple astrophysical messengers (neutrinos, cosmic rays, gamma-rays) in the Gamma-Ray Burst (GRB) internal shock scenario, focusing on the impact of the collision dynamics between two shells on the fireball evolution. In addition to the inelastic case, in which plasma shells merge when they collide, we study the Ultra Efficient Shock scenario, in which a fraction of the internal energy is re-converted into kinetic energy and, consequently, the two shells survive and remain in the system. We find that in all cases, a quasi-diffuse neutrino flux from GRBs at the level of 10−11– (per flavor) is expected for protons and a baryonic loading of 10, which is potentially within the reach of IceCube-Gen2. The highest impact of the collision model for multi-messenger production is observed for the Ultra Efficient Shock scenario, that promises high conversion efficiencies from kinetic to radiated energy. However, the assumption that the plasma shells separate after a collision and survive as separate shells within the fireball is found to be justified too rarely in a multicollision model that uses hydrodynamical simulations with the PLUTO code for individual shell collisions.


Introduction
Gamma-ray bursts (GRBs) have been proposed as plausible candidates for the origin of the Ultra-High Energy Cosmic Rays (UHECRs) and neutrinos, invoking photohadronic interactions in the fireball scenario (Waxman & Bahcall 1997). It is however evident from GRB stacking searches that GRBs cannot be the dominant source of the observed diffuse astrophysical neutrino flux (Abbasi et al. 2012;Aartsen et al. 2017). In radiation models of high-luminosity GRBs where all emission regions look alike (one-zone models), the nominal predictions for the neutrino flux are in tension with these stacking searches in spite of recently improved neutrino-flux estimates (He et al. 2012;Hummer et al. 2012;Li 2012). Dedicated scans of the source parameters (including the baryonic loading, the luminosity in protons versus gammarays) and fits to the UHECR spectrum and composition data confirm that the simple one-zone emission picture is in tension with neutrino data for most of the parameter space Biehl et al. 2018).
As a possible solution, multicollision models have been proposed, which tend to predict a lower neutrino flux (Bustamante et al. , 2017Globus et al. 2015) for the same baryonic loading. These models also demonstrate that the different messengers (neutrinos, cosmic rays, gamma-rays) are predominantly emitted at different radii within the GRB jet . Also, the stochasticity pattern and the time delays between different energy bands in the light curves can contain information on the neutrino production efficiency (Bustamante et al. 2017).
As specific implementations of the fireball model (Rees & Meszaros 1992, 1994, most GRB multicollision models invoke emission from internal shocks (Kobayashi et al. 1997;Daigne & Mochkovitch 1998), in which a set of plasma shells is emitted from an intermittent central engine with a specific distribution of Lorentz factors. The emitted shells can have, for instance, equal masses or equal energies, which are related through the relation = G E m kin . Because of the different velocities, the shells will eventually catch up with each other and collide inelastically, forming a merged shell that continues to propagate through the jet. In each shell collision, shocks form that accelerate particles, converting internal energy (from the inelastic collision) into nonthermal radiation.
In principle, multicollision models provide a natural explanation for the the fast time variability and variety of light-curve shapes observed in GRBs; (see, e.g., Kobayashi et al. 1997;Daigne & Mochkovitch 1998;Beloborodov 2000;Spada et al. 2000;Kobayashi & Sari 2001;Daigne & Mochkovitch 2003). A fundamental problem is, however, the moderate dissipation efficiency of kinetic energy into nonthermal particles. On one hand, the afterglow emission caused by the shells running into the circumburst medium after the prompt phase might be too high if a large amount of energy remains in the jet. On the other hand, extremely powerful engines are required to reach the observed gamma-ray luminosities. The Ultra Efficient Shock scenario has been proposed as a potential solution to this problem (Kobayashi & Sari 2001), in which a fraction of the internal energy is reconverted into kinetic energy. Consequently, the two shells survive and remain in the system after the collision, resulting in swift thermalization. In this scenario, a high overall dissipation efficiency can be achieved, even if the efficiency of each individual two-shell collision is low. A key ingredient of the scenario is the assumption that two shells remain after each collision, whereas hydrodynamical studies and analytical estimates (Kino et al. 2004) demonstrate that this assumption is dependent on the collision parameters and that most collisions likely result in only one merged shell.
In this work, we revisit the multicollision multi-messenger models in Bustamante et al. (2015Bustamante et al. ( , 2017 from the point of view of the collision dynamics and study the impact on the production of multiple messengers. We employ the methods from Baerwald et al. (2012), Hummer et al. (2012), and Biehl et al. (2018) for the radiation calculations and use broken power-law target photon spectra that resemble observations. After recapping the previous approach, i.e., the fully inelastic case with equal-mass injection in Section 2, we scan the parameter space for the two-shell collision for configurations with two post-collision shells and scrutinize the plausibility of this assumption in the fireball evolution. The multi-messenger production is studied in Section 4 for alternative collision models: (A) a version of the reference model, in which a fraction of the internal energy is converted into adiabatic expansion of the merged shell; (B) the ideal Ultra Efficient scenario according to Kobayashi & Sari (2001); and (C) a "hybrid" model, where PLUTO is used to determine the fate of each two-shell collision individually. We discuss the impact on the observables (light curves and neutrino fluxes) in Section 5 and present our discussion and conclusions in Sections 6 and 7, respectively.

Methods and Reference Model
As the reference case, we pick the collision model from Bustamante et al. (2017), based on Kobayashi et al. (1997), and recapitulate the main ideas and mathematical expressions before entering the discussion of alternatives in the next sections.
The relativistic outflow of the jet is discretized as a onedimensional sequence of spherical plasma shells with different velocities and masses. When a faster shell catches up with a slower one, they interact, and collisionless shocks convert kinetic energy into internal energy (which will be radiated as nonthermal particles). During the interaction, the two shells merge into a single new shell, which will continue to propagate as a part of the jet (see Bustamante et al. 2017 and first subsection for a detailed description). The internal energy remaining from the collision is converted into radiation through the interaction of accelerated particles, based on the model from an updated version of the NEUCOSMA software (Biehl et al. 2018). 4 In this section, the "GRB1" in Bustamante et al. (2017) acts as a reference parameters set but with an equal-mass instead of equal-energy setup, resulting in a few quantitative differences, which are discussed in Appendix C.1. It will become clear in the following sections why this choice allows for a better comparison with the alternative models. For the sake of simplicity, we focus on the proton-only case in this study. Changing the composition would not affect the level of the neutrino flux, as it mainly depends on the energy injected per nucleon. However, the maximal energy of neutrino flux is lower for heavier injection, as the maximal rigidity is reduced by Z/A. See Biehl et al. (2018) for a detailed discussion of nuclear injection in GRBs.

Collision Model
A plasma shell is characterized by its mass (m), width (l), and the Lorentz factor (Γ) with respect to the engine frame. When a fast ("rapid," index r) shell collides with a slow (index s) shell, forward and reverse shocks develop. They propagate from the contact discontinuity ("CD," which is the surface separating the initial slow from the fast shell) through the joint density profile initially compressing the matter (see Figure 2 for an illustration of the setup). After the shock waves reach the edges of the density profile, two rarefaction waves develop that propagate back towards the CD. They decompress the matter and reconvert internal energy into kinetic energy. The kinetic properties of the kth shell are determined by the Lorentz factor Γ k , the mass m k or the kinetic energy E k kin, , the radius R k , and the width l k . The other properties, such as the speed β k , the volume p = V R l 4 k k k 2 , and the density ρ k , can be derived. In the following, all unprimed quantities (unless explicitly given in the observer's frame) are given in the engine frame, and all primed quantities refer to the shock rest (or merged shell) frame.
The simulation of the reference model initially contains 1000 shells with their Lorentz factors Γ k ʼs sampled from a lognormal distribution with Γ 0 =500 and A=1.0. We further assume that the spatial separation for all shells is equal to their width , equivalent to the temporal separation d = t 0.01 s k for relativistic shells) and that the innermost shell starts at a radius R min . The initial distribution of kinetic energies E kin,k is a free model choice, and typical assumptions are equal shell masses m k , energies E kin,k or densities ρ k . As previously mentioned, the reference model uses the equal-mass case. From this initial setup, the system evolves until all shells have reached the circumburst medium at R max .
For the collision between a rapid and a slow shell, the properties for the merged shell (index m) are determined from simple analytical expressions following Kobayashi et al. (1997). The Lorentz factor Γ m follows from the conservation of momentum as (using the approximation β=1-1/(2Γ 2 ) that holds for Γ?1). The internal energy available for radiation is then the difference of kinetic energies before and after the (inelastic) collision: We define the efficiency of each collision as the ratio between the dissipated energy and the produced internal energy: The update impacts the description of the optically thick (to photohadronic interactions) case; see Appendix C of Biehl et al. (2018) for details. This leads to slightly lower predicted cosmic-ray fluxes at the highest energies. The baryonic loading is defined with respect to the injection luminosity (previously, the steady-state density), resulting in a more transparent evaluation of the energy budget.
The idealized assumption of η=1 in the Reference model will be modified for the alternative models in Section 4. The Lorentz factors of the forward (fs) and the reverse shock (rs) are determined from The timescale of emission t em is estimated from the time taken by the reverse shock to cross the rapid shell This timescale is observed Doppler-boosted in the observer's frame. The width of the compressed merged shell is computed from For Equations (5)- (7), we follow the description given in Kobayashi et al. (1997). This formulation might not be valid in the nonrelativistic regime; however, collisions in the nonrelativistic regime dissipate little energy and, therefore, contribute only marginally to the multi-messenger emission. Even though t em and l m might be slightly misestimated for those collisions, the impact on our main results is therefore small.

Radiation Model
The collision parameters derived in the previous section are now used for the computation of the energy density and the emission spectra in the merged shell (primed frame). It is assumed that fractions of the injected luminosity are distributed into protons (ò p ), electrons (ò e ), and the magnetic field (ò B ), such that ò p +ò e +ò B =1. Normally, it is assumed that the nonthermal electrons loose energy quickly via synchrotron radiation, implying that gammarays will carry a comparable amount of energy ò γ ;ò e . We assume a baryonic loading of ò p /ò γ =10. It is convenient to relate these quantities to the (observed gamma-rays), consequently ò γ =ò B =1/12 and ò p =10/12.
We assume a power-law spectrum motivated by Fermi shock acceleration ( ( ) ,max ) for the injected proton component. The maximal energy ¢ E p,max is found by comparing the timescales of efficient acceleration ( ¢ R L is the Larmor radius) and photohadronic and adiabatic energy losses (for a detailed discussion of the maximum energies, see also Samuelsson et al. 2019). The proton spectrum is normalized to the available luminosity fraction (ò p ) as outlined above.
For the target photon spectrum, we do not perform explicit self-consistent radiation calculations but instead assume a shape motivated by observations-a broken power law with spectral indices −1 and −2 below and above the break energy  e ¢ 1 keV break , respectively (Gruber et al. 2014). The maximal photon energy is assumed to be limited to e ¢ = 10 GeV max 6 and reduced in case γγ-annihilation sets in at lower energies. The normalization is performed equivalent to the proton case with the fraction ò γ . With these power-law indices, the energy densities depend (at most) logarithmically on the maximal and minimal energies, reducing the impact of our choice of e ¢ max . Successful modeling of GRB prompt spectra within the internal shock model has been performed in, e.g., Bosnjak et al. (2009), Daigne et al. (2011), and Bošnjak & Daigne (2014. Realistic values for the synchrotron peak energy and spectral indices can be achieved, for instance, by allowing a small fraction of electrons to be accelerated to high energies (as in Eichler & Waxman 2005;Spitkovsky 2008). However, an explicit modeling of photon synchrotron spectra is not within the scope of this work. The coupled proton-neutron system is evolved to the steady state using NEUCOSMA (Biehl et al. 2018), which takes photohadronic, pair production, adiabatic, and synchrotron losses into account. The neutrons escape freely, whereas the protons are assumed to be magnetically confined and to escape only from the boundaries (within their Larmor radius), referred to as "direct escape" (Baerwald et al. 2013). This assumption implies that close to the photosphere, the cosmic-ray emission is dominated by neutron escape and for large collision radii by direct proton escape . The mechanism for cosmic-ray escape is currently in discussion, and it is likely that several competing components contribute; see discussion in Zhang et al. (2018). The implementation for secondary particle emission is described in great detail in Biehl et al. (2018).
The model does not account for the emission from subphotospheric collisions, for which the optical thickness to Thomson scattering is larger than one-resulting in a different, thermalized shape of the target photon spectra. Even if cosmicray acceleration could take place at such low radii, the highradiation densities will prevent the particles from reaching the UHE range. Since the pion production efficiency scales with the density similar to the Thomson optical depth, the neutrino production is most efficient close to the photosphere. There could be a significant contribution from subphotospheric collisions; hence, our neutrino-flux estimate shall be considered as minimal prediction. We chose examples in which the fraction of subphotospheric collisions is small to reduce the impact of this effect.

Discussion of Reference Model
We will use the format of Figure 1 to characterize the behavior of different models and, hence, explain it here in greater detail. Since the initial configuration of the 1000 shells is drawn from the distribution in Equation (1), the model naturally produces stochastic fluctuations in the variables. We compute 100 representations if not otherwise noted. In the figures, we show the average as a solid curve and the region between the edge cases with a shaded band.
The left panel shows the differential energy dissipation in the engine frame for the different messengers: neutrinos, UHECRS above 10 10 GeV, gamma-rays, and gamma-rays above 1 GeV only. In the dark shaded region, all collisions occur below the photosphere, in the light shaded region, the collision may be above or below the photosphere (depending on the other shell parameters). The result is consistent with earlier works (Bustamante et al. 2014(Bustamante et al. , 2017: neutrinos originate near the photosphere due to the high density; UHECRs prefer intermediate collision radii since high magnetic fields are required for efficient acceleration-but not so high that radiation losses (such as synchrotron radiation) limit the maximal energy; gamma-rays trace the region where most energy is dissipated (see middle panel and discussion below); HE gamma-rays above 1 GeV in the source frame prefer larger collision radii where maximal energies are not dominated by losses. The different astrophysical messengers originate from different radii of the same jet, thus drawing attention to the collision model that defines the connections among these different regions. The statistical fluctuations from the initial shell setup (shaded areas) imply some variability in this picture, but they do not change it qualitatively.
The middle panel of Figure 1 shows the distribution of radius and dissipated energy of the individual collisions in arbitrary units. This figure demonstrates that most of the relevant collisions happen above the photosphere and there is a not very noticeable anticorrelation between R C and E diss , confirming the relation to the gamma-ray output in the left panel.
The maximal cosmic-ray energy in the upper right panel of Figure 1 exhibits a similar shape as the UHECR curve in the left panel, which only accounts for cosmic rays with E p,max within the red-shaded area. The lower panel shows the magnetic field B′ (left axis, black curve) and the optical thickness to photohadronic interactions evaluated at (right axis, green curve). Both scale with the collision radius, tg The collisions close to the photosphere are optically thick to photohadronic interactions and are responsible for most of the neutrino production (see left panel).
We define the overall energy-dissipation efficiency from kinetic energy to radiation as In contrast to η, it describes the efficiency of the whole system and is an output instead of an input parameter. For the Reference model, we find ò≈36%, i.e., somewhat higher than for the equal-energy setup in Bustamante et al. (2015Bustamante et al. ( , 2017. The reason for this is that the efficiency for individual collision η is higher in the equal-mass case (Kino et al. 2004). Compared to the equal-energy setup, the Reference model discussed here relatively efficiently converts the kinetic energy drawn from the GRB engine into secondaries. We discuss this issue in greater detail in Section 4 for alternative model assumptions, and we compare to the literature in Section 6.1. The kinetic energies assumed in our models (see Table 2) are at the higher end of observations (Gruber et al. 2014), motivated by the higher baryonic loading required to reach the UHECR energy density observed at Earth.

Probability for Two-shell Final States
The key ingredient of the the Ultra Efficient model by Kobayashi & Sari (2001) is the emergence of two shells after a collision. Here, we study the probability of such events using ranges of collision parameters from stochastic GRB multicollision models with hydrodynamical simulations; see Appendix A for details. We use the PLUTO code (Mignone et al. 2007) with a one-dimensional setup. A study for a twodimensional setup demonstrated qualitatively similar behavior to one dimension (Mimica et al. 2004). Magnetic fields are neglected. The colliding shells are assumed to be cold. For a treatment of arbitrarily hot plasma shells see, e.g., Pe'er et al. (2017), who find a possible suppression of the shock formation and a dependence of the energy per particle after the collision on the pre-collision plasma temperature.
The collision process and the post-collision shell configurations are illustrated in Figure 2. Panel(a) demonstrates the equal-energy case (in the source frame) and panel(b) the equalmasses case. In both panels, the shells are shown at ¢ = t 0, the slow shell is depicted in blue, the fast one in red. The massdensity profiles at ¢ t shock , when both shocks have crossed their respective shells, are shown in purple. In Figure 2(a), the resulting mass-density profile is clearly single-peaked, while in Figure 2(b), two distinct peaks moving at different speeds can be seen. In order to identify two-shell post-collision configurations, we evaluate the mass-density profile at the time when both shocks have crossed the shells and the internal energy is maximal. For a given snapshot in time, the mass density can exhibit multiple peaks; Kino et al. (2004) predict up to three peaks. Theoretical estimates comparing the timescales of the wave propagation (when the shock and rarefaction waves cross the two shells) may result in unrealistic approximations, since double-peaked profiles can rapidly evolve into single-peaked but for a constant mass injection rate. Panel (a): Energy dissipated in neutrinos, high-energy (HE) gamma-rays, gamma-rays, and UHECR as a function of collision radius. We define "HE gamma-rays" as gammarays that have energies above 1 GeV in the source frame, and UHECRs as cosmic rays that have energies above 10 10 GeV. In the dark-gray shaded area, only subphotospheric collisions occur, while in the light-gray shaded area, both sub-and super-photospheric collisions occur; subphotospheric collisions are not included in this model. Panel ones. After ¢ t shock , the density profile continues to evolve since the velocity profile is not uniform across the shell(s) (see the lower panels of Figure 2). Instead, the shell edges move in opposite directions in the CD frame. This leads to a dilution of the density profile as time evolves, invalidating the assumption of a constant shell width for the entire duration of a simulation (see also Pe'er et al. 2017). As discussed below, the dissipation of internal energy reduces this effect by slowing down the thermal expansion that ultimately leads to a washed-out singlepeaked profile.
For the classification of a two-shell collision, we define the relative depth of the dip between two mass peaks where ρ Min is the density at the dip between the two maxima and ρ Max the density at the lower peak (illustration see Figure 2 (b)).
The main impact on the post-collision mass-density profile comes from the pre-collision mass ratio (m r /m s ), the Lorentz factors (Γ r /Γ s ), and the shell widths (l r /l s ). We fix the width ratio l r /l s =1 since in GRB internal shock models (see, e.g., Kobayashi et al. 1997;Daigne & Mochkovitch 1998;Sari 2001 andBosnjak et al. 2009;Globus et al. 2015;Bustamante et al. 2017), plasma shells are ejected at constant time intervals resulting in equal widths. An alternative choice l r =0. 1 l s is discussed in the Appendix. Figure 3 shows the depth d as a function of m r /m s and Γ r /Γ s for the parameter ranges that enclose most of the collisions in our models (see Section 4). As already noticed by Kino et al. (2004), pronounced dips, and, hence, double-peaked density profiles, occur for almost equal shell rest masses and for high Lorentz factor ratios Γ r /Γ s . The latter can be understood from the shock-/rarefaction-wave timescales: if Γ r /Γ s is high, then  ¢ ¢ l l s r . Therefore, the reverse shock takes substantially longer to cross the fast shell than the forward shock to cross the slow one. As a result, the rarefaction wave from the direction of the slow shell has enough time to create a pronounced dip separating the two shells.
The radiation model assumes that a fraction of the internal energy is converted into nonthermal electrons and/or ions that leave the system, which effectively cools the system. To study the impact on the collision dynamics, we introduce a simplified energy-dissipation term in our PLUTO simulations that removes internal energy from the system until a certain threshold has been reached (for more details see Appendix A). The rate of this process is assumed to be proportional to the available internal energy. The impact on the occurrence of two-shell configurations is demonstrated by the contours in Figure 3(b) and (c) for two efficiency choices. The dissipation of internal energy reduces the velocity and amplitude of the rarefaction waves, since those are powered from the available internal energy. Too slow rarefaction waves result in shallower dips in the mass-density profile (see also Appendix B), decreasing d and moving up the contours in Figure 3.
For the Ultra Efficient shock scenario, this result means that the dissipation of internal energy into nonthermal particles significantly reduces the probability for two-shell final states. At the same time, the reduction of internal energy available to drive the rarefaction waves increases the lifetime of a two-shell configuration since the thermal expansion of a the doublepeaked profile slows down. Hence, a higher energy-dissipation rate reduces the emergence probability of two-shell configurations but increases their lifetime. The model "PLUTO" described in the next section takes the first effect more rigorously into account.

Alternative Collision Models
In this section, we define three alternative models that use different assumptions for the collision of two shells and study the impact on the multi-messenger production. The format of the discussion follows Section 2 where the details on the Collision of a rapid (r) and slow (s) shell in frame of the contact discontinuity (CD) in the case of equal initial energies (a) or equal masses (b). The massdensity profiles at t′=0 are in red and blue, and the hatched common profile is a snapshot at ¢ t shock after the collision. The velocity profile β of the hatched surface is shown in the lower panels in units of c. The Lorentz factors for the initial shells are Γ r =500 and Γ s =100, and the widths are equal in the source frame (l s =l r ). The snapshot (after the collision) time ¢ t shock is defined when the reverse and forward shocks have crossed the respective shell. The parameter * is the depth of the density dip (where ρ Max the density at the lower peak), which will be used later as a discriminator between single-and double-peaked profiles.
Reference model can be found, which is used as benchmark case.
The Reduced Efficiency model assumes that in each collision, a fraction η=0.5 of the internal energy in Equation (3) is dissipated as radiation, whereas the remaining kinetic energy goes into the adiabatic expansion of a single merged shell; the Ultra Efficient model assumes partially inelastic collisions in which the remaining internal energy is reconverted into kinetic energy, always resulting in two shells after the collision (Kobayashi & Sari 2001); in the PLUTO model, each two-shell collision is simulated individually with the PLUTO code.
Since each model would produce a different result for the same initial shell setup, we modify the setup of each model to re-produce comparative burst durations, variability times derived from the light curve, and total gamma-ray luminosities. Here, the variability timescale refers to the fastest time variability observed on top of longer-lasting light-curve pulses. By construction, all GRBs are normalized to release the same amount of energy in photons in the optically thin regime. The burst duration is determined by the initial size of the system (the sum of all initial shell widths and separations), which matches among the different models. The time variability primarily depends on the number of collisions for a constant burst duration, which scales with the number of initial shells. However, in the Ultra Efficient model, shells do not merge when colliding, leading to substantially more collisions for the same number of initial shells. We compensate for this by reducing = N 1000 initial shells -125. While the resulting variability timescales are rather small, they do not necessarily translate into observed time variabilities, since these are limited by the instrument's response (i.e., the actual variability could be smaller). We have verified that our results related to the multimessenger production are not qualitatively affected by the choice of 1000 initial of shells; see Appendix C.3. In all cases, we assume an engine with a constant mass outflow and constant up-and downtimes resulting in equal initial shell widths and separations. As previously discussed, a constant mass outflow is more likely to result in a splitting into two shells, which the Ultra Efficient model is based on. An overview of the model parameters is given in Table 1.

Reduced Efficiency Model
The Reduced Efficiency model closely follows the Reference model from Section 2 except that in each collision, only a fraction η=0.5 of the internal energy is dissipated as radiation. In the other models, η necessarily has to be smaller than one, as some energy is re-converted to kinetic energy. We therefore chose this model to quantify the impact of the reduced efficiency η compared to the Reference model (η=1). During the collision, the merged shell gets compressed according to Equation (7); we assume that the emission takes place at this time. The remaining internal energy results in thermal expansion of the shell to the width = + l l l m r s , which is the largest possible value that will never overlap with neighboring shells. The reduction of η translates into a reduced overall efficiency ò;18%, see Table 2. As we choose to normalize the total gamma-ray output to the same value for each model, the initial kinetic energy budget has to be increased. As a consequence, the shells have higher density, and the photosphere moves to somewhat larger radii; see Figure 4, first row.
In comparison to the Reference model, differences are visible mainly at large collision radii; see the solid versus dashed curves in Figure 4. We observe lower magnetic fields, slightly decreased optical depths, higher maximum proton energies, and  (9)) as a function of the pre-collision mass ratio (m m r s ) and the ratio of Lorentz factors (G G r s ). In all cases, we fixed Γ s =100 and l r =l s . The system is evolved until both shocks have crossed the respective shells. For panel (a), energy dissipation in not included, in panel (b) 20% of the (theoretically) produced internal energy is removed from the system, and in panel (c), 50% is removed. The red (black) curve corresponds to the case where the colliding shells share the same initial kinetic energy (rest mass). The dots represent the cases shown in Figure 2. shells is the number of shells produced by a single collision. In each collision, the dissipated energy is given by η≡E diss /E int , and the width of the shell(s) after the collision is given by  Note. The energies n E ,tot iso , g E ,tot iso , and E p,tot iso are the total energies released in neutrinos, gamma-rays, and UHECRs above 10 GeV 10 in the engine frame (where only super-photospheric collisions are taken into account), and E p,max is the maximal cosmic-ray energy achieved in the engine frame. We show the mean value over all statistical realizations as well as the standard deviation obtained from the Monte Carlo simulation. generally more energy released in UHECRs. This effect comes from the thermal expansion of the merged shells that gain about a factor of two in width per collision, whereas in the Reference model, the merged shell is still compressed. Consequently, the radiation density decreases at later stages of the fireball evolution, which is populated by previously collided shells. Since the maximal proton energies are higher and the photohadronic losses lower, the energy dissipated in UHECRs increases.

Ultra Efficient Model
In the Ultra Efficient Shock scenario proposed by Kobayashi & Sari (2001), the fraction 1−η of E int is re-converted to kinetic energy, which causes the merged shell to split into two shells. The dissipated energy E int is calculated from Equation (3)  As in the Reduced Efficiency model, the particle acceleration and emission is assumed to take place in the compressed merged shell, calculated as in Section 2. After the collision, both shells are assumed recover their initial masses and widths; see Kobayashi & Sari (2001) for a more refined discussion of these assumptions.
In the Ultra Efficient model, shells repeatedly bounce off each other, causing the fireball to thermalize and leading to an efficient conversion ò;36% of the kinetic energy into radiation. For the same moderate dissipation per collision η=0.5, this model yields twice the overall efficiency of the Reduced Efficiency model that reaches ò;18%. The consequent reduction of the engine's power requirements leads to a smaller deposition of kinetic energy in the afterglow.
Due to the distribution of fewer initial shells across the same system size (to match the total burst duration), the widths and, more importantly, the separations are increased. Compared to the Reduced Efficiency model, the shells travel longer before collisions appear and thus the bulk of collisions moves out to larger radii. This qualitative change of collision patterns and of average collision radii is visible in the middle panel of the second row in Figure 4 (compared to the upper middle panel). Because the fireball thermalizes, the relative Lorentz factor difference of the colliding shells shrinks, leading to less dissipated energy per collision at late times of the evolution. Therefore, the collisions at large radii insignificantly contribute to the fireball's total emission (left panel of second row). The consequence is a reduced neutrino flux (see Table 2) since neutrinos are predominantly produced in collisions close to the photosphere, where the densities and pγ interaction rates are high (see also discussion in next section).

PLUTO Model
In the PLUTO model, each collision is individually simulated with PLUTO, using η=0.5 and d=0.7 (see Section 3 and Appendix A for more details). We analyze the resulting mass-density profile to obtain the post-collision shell masses, widths, and Lorentz factors by fitting box distributions to the obtained mass-density profile and calculating the weighted average of the Lorentz factor profile. In contrast to the Reduced Efficiency model, the thermal expansion of the shell(s) after the collision is not taken into account.
The result for the PLUTO model is demonstrated in the third row of Figure 4. Note that the ensemble size has been reduced to 20 due to computational complexity. The distribution of collisions resembles the Reference model very closely. Only a small fraction (about 8%) of the collisions results in two postcollision shells that move the result closer to the Ultra Efficient model. The magnetic field strength, the optical depth to pγ reactions, and the maximum proton energy are closer to the Reference model than the Reduced Efficiency model, since shells do not thermally expand after the collision. The overall dissipation efficiency is higher (21% compared to 18%) than in the Reduced Efficiency model; see Table 2.
Despite choosing the initial conditions to enhance the occurrence of two-shell configurations, as in the the Ultra Efficient model, our results indicate that these conditions can only be maintained under special circumstances.

Impact on Observations
So far, we have discussed and compared the evolution of the system, including the production of the astrophysical messengers, in the engine frame. Here, we focus on two observables: the expected quasi-diffuse neutrino flux in comparison with the current stacking limit, and the light curves. If not noted otherwise, we use z=2 in this section, and all quantities are given in the observer's frame.

Neutrino-flux Prediction
A model comparison of the expected quasi-diffuse neutrino fluxes among the different models, including the ranges expected from the ensemble fluctuations, is shown in Figure 5. All predictions are relatively similar, except for that of the Ultra Efficient model, which has a substantial amount of collisions at larger radii, as previously discussed and summarized in Table 2.
The expected neutrino flux is at the level of 10 −11 -----10 GeV cm s sr 10 2 1 1 (per flavor). Note the predicted flux in the equal-mass reference setup is about a factor of two higher than in the equal-energy setup used in Bustamante et al. (2017); see discussion in Appendix C and Figure C1. Since the GRB stacking search is basically background-free due to a probability density function including direction, time window, and energy, we expect an improvement of a factor of seven to ten in the proposed next-generation experiment IceCube-Gen2 (Aartsen et al. 2014). Therefore, if taking the prediction at face value, there is a good chance of finding GRB neutrinos with the next-generation detector.
There are several caveats here: the normalization method has been chosen to be compatible with the IceCube method to compute the stacking limit; namely, the fluence of one GRB is converted into a quasi-diffuse flux with a certain number of GRBs contributing per year (that can be related to the local GRB rate if the cosmological distribution of GRBs is known or measured), the baryonic loading of ten is an ad-hoc choice here, and the cosmological distribution of sources is not taken into account (consistent with the stacking method, because for most GRBs, the redshift is not known). The normalization of the flux prediction is therefore somewhat arbitrary, where some factors (number of GRBs per year and redshift) have been chosen to be consistent with the stacking method and are expected to influence both the stacking limit and flux prediction in a similar way.
The baryonic loading, however, is chosen to be consistent with early estimates of the UHECR energy budget, see, e.g., Waxman & Bahcall (1997), assuming that GRBs are the sources of the UHECRs. Details, however, depend on the UHECR composition (especially in the light of recent Auger measurements; see, e.g., Aab et al. 2017), spectral shape, and the cosmic-ray escape mechanism from the source. Recent approaches therefore fit the cosmic-ray data, obtaining the baryonic loading from a fit (see Baerwald et al. 2015 andBiehl et al. 2018 in the one-zone model for protons and nuclei, respectively, and see Boncioli et al. 2019 for low-luminosity GRBs). A similar approach (without explicit derivation of the baryonic loading and parameter-space scan) has been used for a multicollision model in Globus et al. (2015).
When comparing the energy output of (2-4) ×10 52 erg in UHECRs per GRB ( Table 2, depending on the collision model) with the (3-11) ×10 53 erg to power UHECRs (see in Table 1 in Baerwald et al. 2015, depending on the source evolution), the baryonic loading (times gamma-ray luminosity) in our model will need to be a factor of 7-50 higher if UHECR data are to be fit. This implies that the neutrino limit may be more severe than shown here and that the collision model and escape mechanism in internal shocks can be already tested. A more detailed parameter-space study is beyond the scope of this work and will be presented elsewhere.
If, on the other hand, GRBs are not the dominant source of UHECRs, the GRB neutrino stacking searches allow us to derive an upper bound for the baryonic loading of GRBs and their contribution to the UHECR flux. Independent information can come from the spectral energy distribution (SED), which may be altered in the presence of hadronic processes (for instance, in the Fermi-LAT band; Asano et al. 2011;Wang et al. 2018). It is an open question how large baryonic loadings can be tolerated in self-consistent SED computations of GRBs (see, e.g., Asano & Inoue 2007;Asano et al. 2009;Asano & Mészáros 2014).

Light Curves
Another interesting question concerns the impact of the collision model on GRB light curves that are predicted by the multicollision models. Figure 6 shows two examples generated with the Reference model and the Ultra Efficient model; the Reduced Efficiency and PLUTO model light curves are very similar to those of the Reference model. Recall that the total burst durations are chosen to be equal by construction. We also checked different wavelength bands (not shown here), but we did not observe significant differences (such as time delays) because of the stochasticity of the setup chosen hereconsistent with Bustamante et al. (2017, see discussion therein).
The light curves for the Reference model are composed of numerous thin peaks of comparable height, similar to Kobayashi et al. (1997). In contrast, the Ultra Efficient model produces light curves that are dominated by pronounced, broad pulsesmodulated by some time variability; see also Kobayashi & Sari (2001). The widening of the pulses originates from the wider initial shells and, hence, comes from the properties of the engine. A few collisions with high energy dissipation result in a few pronounced pulses contributing most of the photon counts.
, assuming a rate of  = N 667 identical GRBs per year at z=2 contributing to this flux (see Abbasi et al. 2012;Aartsen et al. 2017). The curves correspond to the average, and the shaded area corresponds to the range of neutrino fluxes obtained in the simulations. The GRB stacking limit, shown for comparison, is taken from the best limit obtained in Aartsen et al. (2017). Despite qualitative differences in the shape of the light curves, the time variability, derived by dividing the total burst duration by the number of collisions, is similar for both models (see Table 2).

Dissipation efficiency
Originally, the Ultra Efficient collision model (Kobayashi & Sari 2001) was introduced to increase the emission efficiency of the prompt fireball phase. This is achieved at the cost of reduced efficiency for individual collisions and a higher count of multiple collisions within the fireball. We indeed find a higher efficiency (see Section 3 and Table 1) but only for a small subset of fine-tuned shell parameters. The more realistic Pluto and the Reduced Efficiency models result in lower efficiencies and, as a consequence, require the initial kinetic energy to be significantly higher for the same gamma-ray output. The Reference model reaches efficiency values comparable to the Ultra Efficient model due to high singlecollision efficiencies under the ideal assumption of η=1. The efficiencies of either model are higher than the 1%-15% discussed, for example, in Pe'er (2015) and in the multimessenger context in Globus et al. (2015). This is related to the combination of a large spread in initial shell Lorentz factors and a equal-mass outflow, two features that are known to increase the efficiency Kino et al. 2004). Those are technical parameters of the internal shock fireball model that, at present, lack clear experimental constraints.
The inferred efficiency of detected GRBs is usually based on observations of the energy content in afterglows. Comparing the inferred efficiency to that of the internal shock model has given rise to criticism. This statement is, however, derived from an incomplete observation of the afterglow spectrum and recent observations of HE emission (Abdalla et al. 2019;Acciari et al. 2019a) support the arguments (Fan & Piran 2006;Beniamini et al. 2015Beniamini et al. , 2016) that these estimates may underestimate the energy content in the afterglow.

Model Assumptions
We discuss our choices of model parameters in the context of previous studies. The choice of 1000 initial shells (leading to approximately 1000 collisions) is at the upper end of the typical range and may lead to a short time variability. The number of collisions (and thus the number of initial shells) and the short time variability are directly related by » N T t coll 90 min,var . When comparing to observations, we find that MacLachlan et al. (2013) show values ranging from » T t 100 90 min,var -1000. In Golkhou & Butler (2014) and Golkhou et al. (2015), the values seem to be more in the range of 10-100, with a significant fraction of values smaller than that. However, this critically depends on t min,var -and in practice, the observable short-time variability is constrained by the instrument's response, limiting the possibilities to rule out very small values. On the same note, there is currently no evidence for a very high initial spread of shell Lorentz factors (as assumed in our log-normal distribution). However, this choice is not uncommon in attempts to reach high dissipation efficiencies (Kobayashi et al. 1997;Kobayashi & Sari 2001;Bustamante et al. 2017). Since, in this paper, we study the impact of modifications to the single-collision model on the multimessenger production and emission, we do not attempt to depart from these previous assumptions for the sake of comparability with the previous results in Bustamante et al. (2017).
Another feature that can be compared to other GRB-related studies is the required engine kinetic energy. For our assumptions, the (isotropic-equivalent) engine power of a few times 10 54 erg is required, which seems to be at the upper end of observations (Gruber et al. 2014). As outlined in Beniamini et al. (2015), the actual kinetic energies may be higher, but the appropriateness of the timescale applied in Beniamini et al. (2015) is disputed. The recent discovery of inverse Compton emission from GRB 190114C also yields values in this ballpark (Acciari et al. 2019b). In general, such energies are unavoidable to power the UHECR flux, see, e.g., Section 2 in Baerwald et al. (2015). At present, there is no direct observational evidence for cosmic-ray emission from GRBs due to the lack of HE neutrino observations in coincidence with detected GRBs. Although, in our model, neutrinos are not expected for the current generation of neutrino detectors; see Figure 5.
In summary, some of our assumptions seem extreme, especially in combination. This has been partially motivated by the comparison to the existing literature and partially by the paradigm that GRBs would be the dominant source of UHECRs.

Summary and Conclusions
We have studied the production of multiple astrophysical messengers (gamma-rays, neutrinos, and UHECR protons) in the GRB internal shock scenario. We have used a multicollision scenario involving a set of plasma shells emitted from the central engine with different Lorentz factors, where we have computed the production/emission of the messengers in each collision individually. The focus of this work has been the impact of the collision model between two shells on the multimessenger production in the entire fireball, where we have tested different scenarios: Reference: Plasma shells merge (inelastic collision) and all internal energy is dissipated in nonthermal radiation. Reduced Efficiency: Plasma shells merge and a fraction of internal energy is dissipated (the rest goes into thermal expansion). Ultra Efficient: Plasma shells do not merge and a part of the internal energy is re-converted into kinetic energy that separates the shells. PLUTO: Fate of plasma shells determined for each collision individually from hydrodynamical simulations.
We have also tested different assumptions for the behavior of the central engine (outflow at equal-mass or equal-energy rate), while fixing the observables of the GRB (gamma-ray luminosity, duration, time variability, observed broken power-law spectrum) for better comparability. Compared to earlier works, we have also studied the impact of ensemble fluctuations from the stochastic setup on the results. For all models, we have recovered the qualitative behavior of the multi-messenger emission. The different astrophysical messengers are emitted from different regions of the same object: neutrinos come from the innermost collision radii where the radiation densities are highest, gamma-rays come from a wide range of collision radii where most energy is dissipated, HE gamma-rays prefer large collision radii because they can escape the optically thin shells, and UHECRs come from intermediate radii requiring balanced magnetic field strengths that sustain acceleration up to high energies without introducing too strong synchrotron losses. Additional contribution to the neutrino flux may come from below the photosphere; given the lack of observational constraints on the photon density and the maximal proton energies, the current modeling framework does not permit us to include this component. Thus, our neutrino-flux estimations are conservative, lower limits in this regard.
Substantial differences in the distribution of collision radii have been found in the Ultra Efficient case, coming with a quantitative impact on, e.g., the neutrino flux, UHECR production radius, and shape of the GRB light curve. In that case, a smaller number of initial shells (with larger widths and inter-shell spacings) bounce frequently off each until the system thermalizes. Consequently, more collisions occur at large radii compared to the other models, while more energy is dissipated in the first collisions. For this case, we have found the highest dissipation efficiency (kinetic energy converted into radiation) of about 40% and about a factor of two lower neutrino flux.
We have scrutinized the assumptions of the Ultra Efficient scenario in hydrodynamical simulations with PLUTO. We have found that the bounce back (instead of merging) of shells only hold for very specific conditions: a collision only results in two shells if energy dissipation is not included, the shells are set up with roughly equal masses, and there is a large ratio between Lorentz factors of the rapid and the slow shell. These optimal conditions occur less likely at later stages of the fireball evolution where multiple preceding collisions gradually thermalize the system. In that case, most collisions result in a single dispersing shell. We have coupled the simulation of the fireball with the explicit PLUTO simulation for the collisions of individual shells, and for this study, we included energy dissipation in the PLUTO simulations. This has demonstrated that the emergence probability of two-shell configurations-as assumed in the Ultra Efficient scenario-is very rare (∼8%), which means that our Reference model result is roughly recovered. As a result, the promise of an increased efficiency (like in the Ultra efficient model) can't be kept, and the efficiency-related problems of the fireball model remain unsolved.
In all cases, the expected neutrino flux has been in the range 10 −11 -----10 GeV cm s sr 10 2 1 1 (per flavor) for a baryonic loading of ten, which is potentially in reach of the nextgeneration neutrino observatories, such as IceCube-Gen2. However, we have also noted that the normalization of this flux is somewhat arbitrary. For example, if GRBs are to be the origin of the UHECRs, the baryonic loading will emerge as a deduced quantity from a fit to the observed UHECR flux and composition. From energy budget considerations, we anticipate that the UHECR output will not be sufficient to power UHECRs with the current parameters and the assumptions for the UHECR escape. Further studies will therefore be needed to establish when and if neutrinos can rule out the UHECR origin from GRBs in the framework of internal shock multicollision models, and to what extent an observation will depend on the UHECR composition, the escape assumptions, and the collision model. factor of the fast shell in the CD frame is above two (G ¢ > 2 r ) andĝ = 5 3 (nonrelativistic) else. The simulations are performed using Cartesian geometry, a linear reconstruction, Hancock time-stepping, and the CD-restoring HLLC approximate Riemann solver (for more details about these methods, see the PLUTO user guide).
The simulation starts at t′=0, the time at which both shells get in contact. The coordinate system is defined such that the boundary between initial shells is at x′=0 for t′=0. The fast shell (x′<0) has a positive velocity and the slow shell (x′>0) is negative. The shells are surrounded by low-density plasma (r r = -10 shell 12 ), moving at the same speed as the neighboring shell. The simulation is stopped after both the reverse and the forward shocks have crossed the shells.
To account for energy dissipation due to acceleration of cosmic rays and/or electrons, we add a term to the differential equation solved by PLUTO. In a very simplified treatment, we assume the dissipation rate of internal energy at a given point in time and space to be proportional to the available internal energy. We withdraw energy from the system until a threshold h E int,th * is reached, setting the dissipation rate to zero afterwards: e l s e A1 diss,total int,th * where ρ is the rest mass density, e is the internal energy per unit mass (thus, ρe is the internal energy density), χ is a freely chosen parameter, and η * controls the total dissipated energy by relating it to E int,th , which is the dissipated energy as predicted by the internal shock model (Equation (3)). We set c h = ¢ t 10 shock * and h h = * in the multicollision modeling.

Appendix B Additional Figures: Hydrodynamic Simulations
As an addition to Section 3, we illustrate the impact of energy dissipation on the evolution of the hydrodynamic system in Figure B1 and discuss the impact of the relative shell widths on the final shell configuration.
In Section 3, we equal the shell widths l r =l s in the source frame. Since the relative shock and rarefaction-wave timescales define the produced mass-density profile (Kino et al. 2004), it is natural to expect an impact of relative shell widths on the final shell configuration. The parameter scan for l r =0.1·l s is shown in Figure B2. In contrast to the findings in Section 3, we observe low Lorentz factor ratios Γ r /Γ s to increase the dip depth d. Without a detailed quantitative discussion, we want to explain this result in a qualitative way, looking at the changing ratios of ρ s /ρ r and ¢ > ¢ l l s r (primed quantities refer to the frame of the CD). While generally for l r =l s , ρ s >ρ r and ¢ < ¢ l l s r , both of these relationships change to the opposite for l r =0.1·l s . As a result, the dip separating the two postcollision shells is created in the slow shell instead of the fast one. While before increasing the ratio of the shock crossing times ¢ ¢ t t rs fs shock, shock, (rs refers to the reverse shock, fs to the forward shock) led to more pronounced dips in the massdensity profile, this will now be the case for ¢ ¢ t t fs rs shock, shock, . As a result, the opposite changes in Γ r /Γ s are favorable for the formation of two, distinct shells, and the plot flips with respect to the y-axis.
This result illustrates the strong dependence of the resulting mass-density profile on the individual shell parameters and the difficulty of finding universally applicable predictions. It supports our approach in the multicollision modeling, where we account for the variety of possible impacts by modeling each collision individually with PLUTO instead of using simplified, universal formulas. Figure B1. Snapshot of the mass-density profile (in the frame of the CD) of a single two-shell collision for three choices of η * in Equation (A1). In this example, we set Γ r /Γ s =6.0 and m r =m s , and the snapshot is taken at ¢ = ¢ = t t 2.8 s shock (which equals the shock crossing time of the reverse shock). For η * =0, both the reverse and the forward shocks have crossed the respective shells. For higher η * , the step in the mass-density profile at the left side of the figure indicates that the reverse shock has not yet completely crossed the slow shell, implying that energy dissipation slows down the shocks. Also, the dip between the two shells, which is clearly visible for η * =0, disappears for higher η * . We conclude that energy dissipation reduces the separation between the shells after the collision.  For completeness, we discuss three additional fireball scenarios: (1) a model corresponding to GRB1 in Bustamante et al. (2017; the Reference model but using a source emitting with a constant luminosity), (2) a PLUTO model, where energy dissipation is switched off in the hydrodynamic modeling and all of the internal energy is converted into radiation, and (3) a model with a reduced number of initial shells, with the same parameters otherwise as the Reference model.

C.1. Equal-energy Case
The change from equal mass to equal energy, as in Bustamante et al. (2017), has a few relevant implications. For comparison, a panel for the equal-energy setup is shown in Figure B3. 5 The main difference is visible in the middle panel: for the equal-energy assumption, the bulk of the collisions occurs at smaller radii; also, a smaller spread in the dissipated energy can be observed. When considering the overall efficiency, the equal-energy scenario is significantly less efficient in converting the fireball kinetic energy into radiation. This effect can be understood by looking at the efficiency of single collisions: if two shells of equal mass collide, their kinetic energy is converted much more efficiently into internal energy than for shells of equal energy (for high Γ r /Γ s , the single-collision efficiency approaches 1 for equal-mass shells and 0.3 for equal-energy shells (Kino et al. 2004); the increased single-collision efficiency is also reflected in a higher overall efficiency; see Table B1). Due to the normalization we apply, a lower overall efficiency translates into an increase of the required engine power and, subsequently, initial shell masses. More massive shells are more dense; therefore, the fireball reaches the optically thin regime at larger radii and fewer collisions at small radii are super-photospheric. Those collisions are efficient neutrino emitters (see Bustamante et al. 2017), which decrease the neutrino flux in the equal-energy scenario with respect to the equal-mass scenario (see also Figure C1). The applied normalization also increases the magnetic field as well as the optical depths (the first one is due to the increase of E int , and the latter is due to the increase of the shell masses). As the optical depth to pγ-reactions limits the maximum proton energy in the collisions, higher values of τ pγ decrease the observed E p,max ; see Figure B3. Comparing the equal-energy case with the models presented in the main text, we conclude that the engine behavior (constant power versus constant mass outflow) has a much higher impact on the observed particle fluxes than the choice of collision model.

C.2. PLUTO without Energy Dissipation
As another addition to the main text, we here discuss a PLUTO model where energy dissipation is not included in the hydrodynamic simulations (h = 0 * in Equation (A1)). This corresponds to a theoretical approach where the production of internal energy is completed before particle acceleration and energy dissipation start. After performing the hydrodynamic simulation, we calculate the internal energy of the merged shell kin,after , where the E kin is the summed kinetic energy of all shells. The dissipated energy is then given by E diss =E int (thus, the efficiency is given by η=1.0). In contrast to the PLUTO model in Section 4, the width of the emitting shell is given by Equation (7). Again, we don't take thermal expansion after the collisions into account since the internal energy that could potentially power this effect goes into radiation.
As far as the Reference model is concerned , we can confirm the expected neutrino fluxes for the same gamma-ray luminosity. With the given normalization, and also the evolution of the magnetic field, the optical thickness to pγ reactions and the maximum proton energies show similar behaviors. This once more strengthens our statement about the validity of the standard collision model. However, an increased amount of collisions (22 ± 1% of all collisions resulted in two post-collision shells) increases the overall efficiency of the fireball by roughly 10%. This is consistent with Kobayashi & Sari (2001) who increase the efficiency of their bursts by the nearly complete thermalization of the system through multiple shell collisions. The result once again confirms that the probability of two post-collision shells increases if energy dissipation is lower or not at all included in the hydrodynamic simulations.

C.3. Reduced Number of Shells
To directly compare with Bustamante et al. (2015Bustamante et al. ( , 2017, we assumed 1000 initial shells throughout the simulations discussed in the main text. However, this assumption results in variability timescales close to the fastest ever observed (see discussion in Section 4). We qualitatively discuss here the impact of a smaller choice of initial shells on the simulations.
Since the rate and spatial distribution of collisions directly impact the light curve, it will look qualitatively different. The shell number can be reduced by either staring with a smaller jet extension or by increasing the separation between shells without changing the total length. The reduction of the initial jet size will result in shorter burst durations, typically dominated by only a few bright pulses. Longer quiescent  phases between single peaks are expected for larger shell widths and separations. Figure C2 shows emission regions of different astrophysical messengers for 100 and 200 initial shells, respectively. Both configurations demonstrate that the production radii of different particles as well as the predicted neutrino fluxes (when normalized to the same gamma-ray luminosity in the optically thin regime) are similar to the Reference model with 1000 initial shells. The large statistical fluctuations for 100 initial shells result in higher sample variabilities of the predicted particle fluxes. Despite wider bands of the neutrino fluxes, and, hence, larger burst-to-burst variations, the model does not exceed the limits by IceCube and the conclusions derived with the choice of 1000 initial shells.