{'fraction_non_alphanumeric': 0.03767064109210299, 'fraction_numerical': 0.01953981736252343, 'mean_word_length': 4.523641703377386, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-messenger light curves from gamma-ray bursts in the internal shock model

Gamma-ray bursts (GRBs) are promising as sources of neutrinos and cosmic rays. In the internal shock scenario, blobs of plasma emitted from a central engine collide within a relativistic jet and form shocks, leading to particle acceleration and emission. Motivated by present experimental constraints and sensitivities, we improve the predictions of particle emission by investigating time-dependent effects from multiple shocks. We produce synthetic light curves with different variability timescales that stem from properties of the central engine. For individual GRBs, qualitative conclusions about model parameters, neutrino production efficiency, and delays in high-energy gamma rays can be deduced from inspection of the gamma-ray light curves. GRBs with fast time variability without additional prominent pulse structure tend to be efficient neutrino emitters, whereas GRBs with fast variability modulated by a broad pulse structure can be inefficient neutrino emitters and produce delayed high-energy gamma-ray signals. Our results can be applied to quantitative tests of the GRB origin of ultra-high-energy cosmic rays, and have the potential to impact current and future multi-messenger searches.


INTRODUCTION
The most energetic particles discovered -ultra-highenergy cosmic rays (UHECRs), above 10 9 GeV, and highenergy astrophysical neutrinos, above 100 TeV -are detected with regularity, but their sources remain unknown. The extreme physical conditions required to reach these energies restrict the potential source classes. Gammaray bursts (GRBs) are one such class. They are the most luminous transient electromagnetic astrophysical phenomena: they emit gamma rays up to ∼ 100 GeV, concentrated in only a few tens of seconds. Their observed high luminosities and inferred high particle densities make them prime candidate sources of UHECRs and neutrinos.
One of the leading explanations of GRB emission is the fireball model (Rees & Mészáros 1992). In it, a central engine -likely, a black hole -injects plasma blobs with different relativistic speeds into the jet flow. When they collide with one another, they create mildly relativistic shocks in which their kinetic energy is transformed into internal energy that is radiated away as high-energy particles (Rees & Mészáros 1994;Paczynski & Xu 1994). If the jets contain enough baryons, then protons (Milgrom & Usov 1995;Waxman 1995;Vietri 1995) and nuclei Wang et al. 2008) can be accelerated by the shocks into UHECRs which, upon interacting with source photons, create pions and other mesons. Their decays produce PeV gamma rays and neutrinos (Waxman & Bahcall 1997).
In recent years, IceCube has started testing the hypothesis of GRBs as sources of high-energy neutrinos (Abbasi et al. 2012;Aartsen et al. 2015aAartsen et al. , 2016. The simplest versions of the fireball model are now in tension with the data. In these, the neutrino flux prediction for a burst is based on the emission from a single representative plasma blob collision in the jet -this is known as the one-zone approach. While high-luminosity GRBs cannot account for the diffuse high-energy astrophysical neutrino signal detected by IceCube (Tamborra & Ando 2015; Aartsen et al. 2016), they persist as interesting objects because of the following two reasons. First, they can still be the sources of UHECRs. This can happen if UHECR emission and neutrino emission are not tightly correlated (Baerwald et al. 2013), a possibility that we consider in the present work. Second, they are attractive targets for neutrino telescopes, as timing and directional gamma-ray information can be used to reduce neutrino backgrounds. They are ideal sites to look for joint gamma-ray and neutrino signals.
In this paper, we focus on the connection between gamma rays, UHECRs, and neutrinos in individual GRBs. Unlike most of the existing literature, we consider the multi-messenger emission from a multitude of different plasma collisions along the jet, each one occurring under different physical conditions -this is known as the multi-zone approach.
We construct synthetic GRB light curves -curves of photon rate versus time -that successfully reproduce generic features of real ones. We identify the properties of the central engine that lead to such light curves, and study the consequences for different messengers. For example, light curves dominated by broad pulses overlaid with fast variability imply a "disciplined" engine that ejects shells with little spread in speed at any given time. In such cases, the neutrino production efficiency can be low, and high-energy gamma-ray signals can reach Earth delayed with respect to low-energy signals. Conversely, light curves with no broad pulse structure hint at likely efficient neutrino emitters. We also test the robustness of the assumptions going into the minimal neutrino flux estimate found in Bustamante et al. (2015). This paper is organized as follows. In Section 2, we comment on GRB emission models. In Section 3, we introduce our multi-zone simulation and a sample of GRBs computed with it. In Section 4, we compute the gammaray and neutrino light curves for each of them. In Section 5 we calculate their associated quasi-diffuse neutrino fluxes and show that different types of particles are emitted from different regions in the jet. In Section 6, we show that, for some GRBs, delays between the light curves in different energy bands are possible, and why. We summarize and conclude in Section 7. Appendix A contains a detailed presentation of the multi-zone collision model used in our simulations. Appendix C contains a discussion of the impact of alternative assumptions for the collision dynamics.

ONE-ZONE VS. MULTI-ZONE EMISSION
GRBs are conceivably fueled by matter accretion, either in the collapse of a massive star -long-duration bursts, lasting > 2 s -or in the merging of two neutron stars or a neutron star and a black hole -shortduration bursts, lasting < 2 s. The central engine, likely a newly formed black hole, emits two relativistic matter jets in opposite directions. When one of them points towards Earth, we might detect it as a GRB. We focus on long-duration bursts because they have higher gammaray luminosity (∼ 10 52 erg s −1 ), represent ∼ 75% of observed bursts, and their emission mechanism has been more deeply studied.
In the internal shock scenario of the fireball model, internal plasma collisions within the jet account for particle emission during the initial, or prompt, phase of the burst, which typically lasts 10-100 s. Ultimately, the jet reaches the circumburst medium, triggering a late emission phase -the afterglow. We will focus exclusively on the prompt phase, where PeV neutrino emission is well-motivated. This is the energy range where existing water-Cherenkov neutrino telescopes, like IceCube and ANTARES, are most sensitive. However, the internal shock scenario, while successful, is not without issues, which we introduce below.
Among alternative theories, photospheric emission models have become popular after the Fermi satellite was launched (Thompson 1994;Rees & Mészáros 2005;Pe'er et al. 2006;Ioka et al. 2007;Giannios 2008;Beloborodov 2010;Lazzati et al. 2013). In them, the bulk of the prompt emission is attributed to quasi-thermal emission from below the photosphere, i.e., from the region where gamma rays are unable to escape due to high optical depth to electron-photon scattering. Amid different versions of photospheric emission models, dissipative photosphere models have been discussed both theoretically and observationally, due to their appealing implications (Beloborodov 2013;Hascoet et al. 2013;Vurm et al. 2013). Some authors consider the superposition of pure thermal emission to explain the prompt emission spectrum (Lundman et al. 2013).
Even though the classical internal shock model has several theoretical issues, internal shocks may still play an important role in the dissipation of the outflowing kinetic energy. For instance, some photospheric emission models require sub-photospheric dissipation that is often attributed to internal shocks (Rees & Mészáros 2005). Magnetic reconnection may also be driven by shocks beyond the photosphere (Zhang & Yan 2011). Even in the optically-thin internal shock model, the observed spectra can be explained by synchrotron emission with some modifications such as stochastic acceleration of electrons in the shock downstream (Bykov & Mészáros 1996;Murase et al. 2012;Asano & Terasawa 2015).
In the internal shock model, shell collisions inside the jet are responsible for the shape of the GRB light curves. They typically have a fast time variability t v , of ∼ 10-100 ms (Bhat 2013;Golkhou et al. 2015), frequently superposed on top of slower pulse structure (Nakar & Piran 2002b;Ramirez-Ruiz & Fenimore 1999). Since the light curve reflects the particularities of the central emitter and jet propagation (Nakar & Piran 2002a;Zhang et al. 2016;Lopez-Camara et al. 2016), no two are equal, though they can be classified on the basis of their morphology.
One-zone collision models assume average shell properties -such as an average speed or Lorentz factor Γ -derived from gamma-ray observations. Neutrino emission is computed for a single representative collision (Dermer & Atoyan 2003;Guetta et al. 2004) and is scaled by the total number of collisions -roughly, T 90 /t v ∼ 1000, with T 90 ∼ 10 s the burst duration -to yield the flux for the whole burst, implying that all collisions are identical (see, e.g., Asano (2005)).
However, it is unrealistic that all collisions occur at the same radius. The main reason is that all shells would have to be emitted with precisely tuned speeds. A different reason comes from the fact that energy deposited in the afterglow cannot be too large; therefore, a sizable part of the energy must be dissipated during the preceding prompt phase. Because dissipation of kinetic energy into radiation scales with the difference of the Lorentz factors of the colliding shells, this implies that the differences must be large. In fact, there is considerable spread in the inferred values of bulk Lorentz factors of GRBs observed by the Fermi satellite (Ackermann et al. 2012). Assuming for the Lorentz factors a broad or bi-modal distribution leads to efficient dissipation of kinetic energy in the prompt phase (Kumar 1999;Zhang et al. 2007). As a consequence, collisions occur at many positions throughout the jet, from below the photosphere -where gamma rays cannot escape -to the circumburst medium, and that particle emission will be different for each collision.
The impact of a distribution 1 of Lorentz factors or emission radii have been addressed qualitatively  and quantitatively Bustamante et al. 2015;Globus et al. 2015). For example, Bustamante et al. (2015) predicted a minimal quasi-diffuse prompt neutrino flux from super-photospheric collisions, at the level of ∼ 10 −11 GeV cm −2 s −1 sr −1 per flavor, which is presently below the sensitivity of IceCube. In later sections, we discuss the relationship among different types of messengers in multiple-collision models, focusing on how the properties of the central engine impact the GRB light curves and neutrino emission efficiency.
IceCube recently discovered a diffuse astrophysical flux of neutrinos between about 30 TeV and 2 PeV (Aartsen et al. 2013b(Aartsen et al. ,a, 2014b(Aartsen et al. ,a, 2015b. However, no individual neutrinos are associated to known GRBs (Abbasi et al. 2012;Aartsen et al. 2015aAartsen et al. , 2016. The resulting upper limit, extrapolated from a catalog of bursts, is about an order of magnitude below the diffuse flux. Thus, it is unlikely that classical high-luminosity GRBs are the main origin of the observed IceCube neutrinos (Murase & Ioka 2013;Laha et al. 2013;Bustamante et al. 2015). (However, low-luminosity classes of GRBs can describe the diffuse neutrino flux without violating the stacking limits Gupta & Zhang 2007), subject to their assumed source density.) Nevertheless, this does not preclude GRBs from being the dominant sources of UHECRs. So far, results of these observational searches have been interpreted mostly in the context of one-zone models (Waxman & Bahcall 1997;Guetta et al. 2004;Li 2012;Hümmer et al. 2012;He et al. 2012). While the plain one-zone ansatz with commonly assumed parameters favorable for the UHECR explanation is starting to be challenged by Ice-Cube data He et al. 2012), it is not fully ruled out yet.
Another important aspect is that the relationship between neutrino and cosmic-ray emission depends on the UHECR escape mechanism (Baerwald et al. 2013). It has been shown that, depending on the parameter values, a fraction of cosmic rays must directly escape from the sources without neutrino production (Baerwald et al. 2013). Earlier approaches focused on the emission of neutrons only (Ahlers et al. 2011), which leads to a stringent relation between the neutrino and cosmic-ray fluxes. Consequently, bounds on GRB neutrinos translate into bounds on the GRB parameter space; see Baerwald et al. (2015) for the bounds in the one-zone model. Further quantitative tests are necessary to test the possibility that UHECRs mainly come from GRBs, especially taking into account the effects from multiple collisions.

SIMULATING COLLISIONS IN A GRB JET
We study the prompt emission phase of the GRB in the internal shock scenario by simulating collisions of many propagating plasma blobs. Our simulations are based on Kobayashi et al. (1997); Daigne & Mochkovitch (1998). See Appendix A for details. Figure 1 illustrates the evolution of the blobs along the jet. To simplify the computations, we take them to be spherical shells. The simulation starts with the central engine emitting N sh ∼ 1000 shells. Each one propagates at a different relativistic speed; typically, they have Lorentz factors Γ ∼ 100. After some time, they catch up to one another, collide inelastically, and merge into new shells. At collisionless shocks generated during the collisions, a fraction of the kinetic energy of the colliding shells is used for particle acceleration and creation, and becomes the internal energy of the newly created shells. The new shells cool instantly by emitting high-energy particles, continue propagating in the jet, and may collide again. Each simulation in this work comprises about 1000 collisions, as almost all initial shells collide.
Non-thermal electrons receive a fraction e of the internal energy of the new shell, protons receive a fraction p , and magnetic fields receive a fraction B . Because electrons cool fast via synchrotron emission, e can also be regarded as the energy fraction of radiated gamma rays. We assume that there is energy equipartition between electrons and magnetic fields, and ten times more energy in protons (Abbasi et al. 2010;Hümmer et al. 2012), i.e., e = B = 1/12 and p / e = 10; see Section A.4. Thanks to the fast cooling, p / e can be regarded as the non-thermal baryon loading factor, defined as the ratio of cosmic-ray energy to radiation energy .
The volume of the shell (see Section A.3) grows with the distance r to the central emitter ∝ r 2 until the radial expansion of a shell becomes important at the shell spreading radius. Except in collisions, the number of particles in a shell is conserved. So, as a shell propagates and expands the density of particles in it falls ∝ r −2 . This affects where in the jet different types of particles are produced and emitted: neutrinos come predominantly from close to the photosphere, where densities are high; ultrahigh-energy cosmic rays come from intermediate collision radii; and high-energy gamma rays escape abundantly from large collision radii .
In a shell collision, if both protons and electrons are accelerated up to sufficiently high energies, pγ interactions create pions via the ∆ + (1232) resonance and other processes. High-energy gamma rays are produced in π 0 → γγ (and by synchrotron and inverse-Compton radiation by electrons), while high-energy neutrinos are produced in π + → µ + ν µ →ν µ e + ν e ν µ (and, via additional production channels, by its charge-conjugated version). The pion production efficiency is, on average, 20%; this is the fraction of proton energy received by pions. It is distributed roughly evenly among the pion decay products, so each neutrino carries ∼ 5% of the proton energy. Therefore, production of PeV neutrinos requires 20 PeV protons.
At the source, for the protons we assume a power-law spectrum ∝ E −2 p exp −E p /E p,max , with E p the proton energy, as expected from Fermi acceleration (primed quantities are in the shock rest frame). The maximum proton energy E p,max is computed by balancing the acceleration rate, and the synchrotron, adiabatic, and photohadronic energy-loss rates (Baerwald et al. 2013).
For the target photons, we assume a broken power-law spectrum, Eq. (A12), which resembles the gamma-ray spectrum observed at Earth. This assumption is certainly justified beyond the photosphere, from where the photons can escape. Note that we neither generate the photon spectrum from first principles, nor explain its origin, which typically includes radiation processes such as inverse Compton scattering, synchrotron emission, and bremsstrahlung. The Band function can be reproduced by invoking various effects, such as slow heating of electrons (Bykov & Mészáros 1996;Murase et al. 2012;Asano & Terasawa 2015). In this sense, our approach is modelindependent; this comes at the expense of insight on the radiation processes at work and their detailed individual implications on secondary production. To calculate the secondary particle production, we use the state-of-the-art NeuCosmA  software (see Section A.5).
UHECRs are emitted via two mechanisms: either pγ interactions transform protons into neutrons which escape the merged shell, and later beta-decay back into protons (Mannheim et al. 2001), or protons directly leak out of the shell without interacting. The latter situation occurs if the proton Larmor radius -determined by the physical conditions in the shell -exceeds the shell width at the highest energies; see Baerwald et al. (2013) for a detailed discussion. This direct escape produces a hard proton spectrum (it is dubbed a "high-pass filter" in Globus et al. (2015)). Direct proton escape dominates over neutron escape only when the photon densities are low. Therefore, neutrino production associated to neutron escape is higher than that associated to direct proton escape (Mannheim et al. 2001).
Gamma rays are produced throughout the jet. However, they only escape if they are produced above the "photosphere", i.e., the radius below which Thomson scattering occurs frequently enough to effectively trap them (see Section A.5). We adopt the pragmatic point of view that solid predictions about GRB dynamics and secondary production must be grounded in gamma-ray observations, which come exclusively from above the photosphere. The drawback of this assumption is that we cannot accurately calculate secondary production below the photosphere, since the sub-photospheric photons are subject to significant thermalization and the observed spectrum becomes quasi-thermal (Beloborodov 2013;Hascoet et al. 2013;Vurm et al. 2013). For our calculations of the gamma-ray flux, we will use only super-photospheric collisions. For neutrinos, we will additionally estimate possible sub-photospheric contributions. Photomeson production is more efficient below the photosphere, so that neutrino emission around the photosphere can be dominant, especially in photospheric emission models (Murase 2008;Wang & Dai 2009;Bartos et al. 2013;Bustamante et al. 2015).
For super-photospheric collisions, the gamma rays are still attenuated in energy and produce particle cascades, via electron-positron pair-production processes inside shells (see below). After gamma rays leave the source, they scatter off cosmological photon fields -microwave, optical, infrared -and their energies are degraded down to a few hundred GeV or lower. UHECRs also lose energy through photohadronic interactions and inelastic scattering off cosmological photons. Neutrinos are affected only by the adiabatic cosmological expansion, which changes their energies in, at most, a factor of a few.
We normalize the gamma-ray emission of our simulated bursts to typical time-integrated GRB gamma-ray energies, E iso γ,norm = 10 53 erg (in the source reference frame). This is equated to the sum energy emitted by all collisions, sub-and super-photospheric (Eq. (A11)). In our simulated bursts, about 50% of the energy is liberated in super-photospheric collisions.

SYNTHETIC LIGHT CURVES
GRB light curves are highly irregular. They are extremely variable, and a fraction of GRBs have minimum variability timescales of ms (Golkhou et al. 2015 Figure 2. Selected GRB light curves, detected by BATSE (Paciesas et al. 1999) in the > 20 keV range (channels 1-4); from left to right: GRB920513, GRB931008, GRB940210, and GRB920627. The implied time resolution is 64 ms.
typically have a duration of ∼ 1 s with an asymmetric structure (Nakar & Piran 2002b;Ramirez-Ruiz & Fenimore 1999), and some GRBs show distinct quiescent time between the pulses. In the internal shock scenario, the shape of the GRB light curve is the result of shell collisions (Kobayashi et al. 1997;Daigne & Mochkovitch 1998;Spada et al. 2000;Kobayashi & Sari 2001;Beloborodov 2000;Daigne & Mochkovitch 2003;Aoi et al. 2010). When shells collide, a gamma-ray pulse and a neutrino pulse are emitted. The superposition of all pulses emitted during the evolution of the burst makes up the light curve. Different features in the light curves are due to differences in the distribution of collision radii, as determined by the initial speeds of the shells, imprinted on them by the central emitter. See Appendix A.4 for an explanation of how we construct the light curves.
The simulation parameters from our benchmark model from Bustamante et al. (2015) are re-branded here as "GRB 1". Note that we slightly adjusted the collision model 2 ; therefore, the results are not exactly the same as in Bustamante et al. (2015). Most notably, collisions occur on average at slightly lower radii now. The initial values of shell Lorentz factors, Γ k,0 , are randomly sampled from a log-normal distribution defined by the characteristic value Γ 0 and the amplitude of fluctuations A Γ : where the random variable x follows a Gaussian distribution, P (x)dx = (2π) −1 e −x 2 /2 dx. When A Γ < 1, the mean value Γ ≈ Γ 0 and the variance ∆Γ ≈ A Γ Γ 0 . When A Γ > 1, the mean value and variance are significantly affected by fluctuations. Large A Γ 1 are typically required for efficient conversion of kinetic energy to radiated energy (Kobayashi & Sari 2001), since the latter is proportional to the difference in speeds between two colliding shells. In addition, A Γ 0.1 is necessary for the number of collisions to be large, i.e., N coll ∼ N sh . The light curve for GRB 1 is in Fig. 4. It is dominated by fast variability, which is determined by the ratio of T 90 to N coll , on the order of tens of ms, which is typical for GRBs. Figure 2 shows selected real gamma-ray light curves that have slower structure overlaid with fast variability. In the internal shock scenario, the slower structure can be produced by modifying the behavior of the engine. We can change the widths of the initial shells, the separations between them or the spread A Γ , emit shells intermittently, ramp up or down the Lorentz factors during shell emission, or make them oscillate (see also Daigne & Mochkovitch (1998)). Table 1 lists our sample simulations, GRBs 1-6, where we have implemented these options We have chosen parameter values such that the associated synthetic light curves reproduce features similar to Fig. 2. Figure 3 shows the randomly sampled initial values of the Lorentz factors Γ k,0 of the shells in GRBs 1-6, as a function of initial radius r k,0 . See Table 1 for descriptions of the underlying distributions. Figure 4 shows the synthetic gamma-ray and neutrino light curves for GRBs 1-6. GRB 1 has fast time variability without prominent features. GRB 2 has a speedup in Γ during shell emission; a single pulse is overlaid with fast variability. For a slowdown, the pulse occurs earlier and the light curve is time-inverted. However, the efficiency is much lower in this case, as the fast shells are emitted first. GRB 3 has three pulses with linear slowdown. The second and third pulses collide with slow shells from the preceding pulse and therefore contribute more strongly to the light curve than the first pulse. GRBs 4 and 6 have more oscillating periods. GRB 5 is oscillating as well, but its lower amplitude increases linearly during emission. Comparison with Fig. 2 reveals that GRB 3 was inspired by GRB931008; GRBs 4 and 6, by case GRB920627; and GRB 5, by cases GRB920513 and GRB940210. Appendix B contains four more simulation examples, with different engine assumptions but similar behavior to that of GRBs 1-6, showing that these are representative.
Most GRB light curves detected by Fermi (Ajello et al. 2013;von Kienlin et al. 2014) do not have prominent features like the ones in Fig. 2. The reason why the latter are featured in the literature more often than simpler single-pulse or fast-variability light curves is possibly a selection effect (i.e., they are more interesting to show and study). From that perspective, it is conceivable that GRBs 3-6 are not "typical" GRBs.
The light curves of GRBs 3 and 5 are qualitatively different from the others in one key aspect: they have a Table 1 Description of simulated GRBs 1-6 Oscillating Γ (four periods) with distribution widening up; in addition, engine intermittent: Nup pulses followed by N down pulses; corresponds to increasing Lorentz factor during uptime Note.
-Common values for all models: N sh = 1000, δteng = 10 −2 s, d = l = c · δteng, r min = 10 3 km, r dec = 5.5 · 10 11 km, z = 2, e = B = 1/12, p = 5/6, η = 1.0 (acceleration efficiency (Baerwald et al. 2013)). See Table 3 for an explanation of each parameter. The period Tp for the oscillating cases refers to Γ changing between Γ 0,1 and Γ 0,2 , and A Γ changing between A Γ,1 and A Γ,2 ; Tp is a fraction of the total number of emitted shells. This means that Γ and A Γ change between first and second value with a factor sin 2 (k/(N sh · Tp) · π), where k is the index of shell (1 ≤ k ≤ N sh ). For GRB 5, the lower amplitude increases from Γ 0,1 to Γ 0,2 linearly with k.  Table 1 for descriptions of the underlying distributions of initial Lorentz factors in each simulation.
dominant, broad pulse structure overlaid with fast time variability, whereas in the other bursts the fast component is more relevant. This feature can be traced back to the input parameters in Table 1: GRBs 3 and 5 have a "disciplined" central engine that emits shells within a narrow Γ distribution (A Γ = 0.1), while the average Γ changes slowly. Therefore, most collisions occur at larger radii compared to the cases where the spread in Γ is larger. We will see that this also affects the neutrino production efficiency in the case of GRB 5. The internal shock model has been invoked as a successful model to explain irregular features of the GRB light curve (Kobayashi et al. 1997;Daigne & Mochkovitch 1998;Spada et al. 2000;Kobayashi & Sari 2001;Beloborodov 2000;Daigne & Mochkovitch 2003;Aoi et al. 2010). However, because the classical internal shock model is being challenged, light curve predictions of alternative models have also been extensively investigated. In particular, it has been shown that turbulence or magnetic reconnection models can better explain the observed structure of the GRB light curve (Narayan & Kumar 2009;Lazar et al. 2009;Beniamini & Granot 2016). Some recent studies suggested that GRB light curves may consist of the superposition of slow and fast components, as inferred by a gradual depletion of the fast component at low energies (Vetere et al. 2006). Sub-structures of the observed GRB pulses can be easily accounted for by assuming relativistic motions in the bulk of a relativistic jet. Although their origins are unclear at present, such relativistic motions could be realized in the ICMART model . Compared to these explanations, our pulse structure comes from the properties of the engine rather than the jet. Figure 4 shows there is no linear correlation between the heights of gamma-ray and neutrino pulses. This is because the height of a neutrino pulse is, via the pion production efficiency, more sensitive to the collision radius than the height of a gamma-ray pulse, i.e., it depends on the proton and photon densities, which drop ∝ R 2 C , where R C is the collision radius.
There are no long time delays between gamma-ray and neutrino pulses: they are within T 90 (see Eq. (A21)) of each other. There may be, however, short delays. For example, in GRB 3, the neutrino peak corresponding to the first large gamma-ray peak is suppressed. So the first gamma-ray detection will have occurred ∼ 20 s, before the neutrino instrument triggers. In GRB 4, the quiescent periods of gamma-ray and neutrino emission coincide, which may be exploited by neutrino telescopes to set further time window cuts.
In GRB 5, the fast time variability gives rise to neutrino spikes, whereas the longer pulses seen in gamma rays are hardly present in neutrinos. Indeed, the rise time of the spike of particle emission associated to one collision depends strongly on Γ (see Eq. (A16)) and, consequently, it is ∝ R C . Therefore, faster rises -sharper structures -are created by collisions at smaller radii, where the neutrino production efficiency is higher. Table 2 lists the parameters output by GRBs 1-6. GRBs 1, 2, 4, and 5 have time variabilities of ∼ 50 ms and durations of ∼ 50 s. In GRB 3, the duration (33 s) and time variability (36 ms) are smaller, as mainly two of the three peaks contribute. In GRB 6, pulses are separated by a downtime and, therefore, duration and variability are larger (98 s and 97 ms, respectively).
The dissipation efficiency ε of a burst in Table 2 is the ratio between total energy dissipated by all types of particles in super-photospheric collisions and total kinetic energy available at the start of the simulation. Most Note. -The parameters are: variability timescale (tv), total energy emitted as gamma rays (E iso γ,tot ), as protons (E iso p,tot ), as neutrinos (E iso ν,tot ), ratio between neutrino and gamma-ray energies (E iso ν,tot /E iso γ,tot ), and overall emission efficiency ε. Energies are computed using super-photospheric collisions only. Only protons in the UHECR energy range, above > 10 10 GeV, are counted. The efficiency ε is defined as the ratio of total energy dissipated by all types of particles (gamma rays, protons, neutrinos) to the total kinetic energy initially available.
simulations 3 have high ε ≈ 11-27%, which reasonably agrees with previous work (Beloborodov 2000). The efficiency is lower in cases with narrow Γ distribution, as expected.
The gamma-ray emission efficiency ε γ ≈ e ε is about 10 times smaller since most of the dissipation energy is assumed to be carried by protons rather than electrons.
Here, e is the fraction of energy in electrons and photons; see Section A.4. Such a value for the radiation efficiency may be too small compared to ones preferred by observations. However, if the internal energy is carried by thermal protons or confined cosmic rays, it is natural to expect the reconversion of the internal energy into the kinetic energy; see Appendix C. It has been shown that this effect increases the gamma-ray emission efficiency, represented by the ratio of prompt gamma-ray energy to afterglow kinetic energy, calculated in an approach where shells reflects off each other after colliding, i.e., collisions are not perfectly inelastic (Kobayashi & Sari 2001). Recent results on afterglow modeling also suggest a small value of the gamma-ray emission efficiency (Beniamini et al. 2015).

MULTI-MESSENGER EMISSION
Here we discuss the emission of multiple messengers and their relation.

Weak vs. strong neutrino emitters
The time-integrated neutrino fluence of a simulated burst, for a baryon-rich jet, roughly scales as (Bustamante et al. 2015) assuming a fixed photon break energy (see Appendix A.4). The first factor gives the fraction of collisions with high optical depth τ pγ 1 to pγ interactions, ε is the energy dissipation efficiency, and E iso γ,tot is the total energy emitted as gamma rays in super-photospheric collisions. Unlike one-zone predictions (Waxman & Bahcall 1997;Guetta et al. 2004;Li 2012;Hümmer et al. 2012;He et al. 2012), the fluence does not depend on the average Lorentz factor of the shells. Figure 5 shows τ pγ as a function of R C for collisions in GRBs 1 and 5. In GRB 1, a strong neutrino emitter, about 60 in 1000 collision occurred above the photosphere -so that they were optically thin to Thomson scattering -and were still optically thick to photomeson production. The other simulations with broad Γ distributions have similar results. Therefore, the first factor in Eq. (2) is ∼ 0.05 for strong neutrino emitters. The energy dissipation efficiency in Table 2 lies around ε = 0.2 for these GRBs. As a result, for fixed E iso γ,tot , the quasidiffuse neutrino flux that we infer below  is relatively robust.
The situation is different for GRBs 3 and 5. They have lower efficiencies ε ≈ 10%. More importantly, they have no optically thick collisions close to the photosphere; see Fig. 5, bottom, for GRB 5. The reason is that they emit shells with a variable but narrow Lorentz factor distribution (see Fig. 3) which tends to induce collisions at larger radii, where photon densities are low and, therefore, so is neutrino production efficiency. In particular, this makes GRB 5 our weakest neutrino emitter, i.e., it has the lowest ratio of emitted neutrino energy to gamma-ray energy beyond the photosphere.
While the same effect should also make GRB 3 a weak neutrino emitter, its neutrino flux is still on a level with our other examples. The reason for this comes from its very specific initial shell setup. It consists of three narrow pulses, each with decreasing Γ. The collisions are therefore dominated by the first, fast shells of a pulse running into the preceding pulse -these shells have largely different Lorentz factors; in particular, the differences are larger than in GRB 5. Most of the neutrino emission comes from these first collisions, which happen below and slightly above the photosphere.
We saw that GRBs with light curves dominated by fast variability are likely to be strong neutrino emitters. The reverse conclusion does however not hold. While both GRBs 3 and 5 have gamma-ray light curves with broad pulses overlaid with fast variability, only GRB 5 is a weak neutrino emitter. Therefore, in the multi-zone internal shock model, we can tell, by inspection of the gamma-ray light curve alone, whether or not a GRB is likely to be a strong neutrino emitter. Conclusions about weak neutrino emitters require a closer inspection of the specific light curve morphology.  The average fluence per collision drops stronger than approximately ∝ R −2 C . In principle, the fluence from sub-photospheric collisions is high, due to the high extrapolated photon density; however, we do not use those collisions in our flux calculations.

Quasi-diffuse neutrino flux
We derive the all-sky quasi-diffuse ν µ +ν µ flux J νµ associated to a particular simulated GRB by scaling its fluence by the rate of long-duration GRBs,Ṅ = 667 per year, i.e., J νµ = F νµ ·Ṅ · (4π) −1 . Since this flux does not contain contributions from sub-photospheric collisions, it is effectively a lower limit on the prompt GRB neutrino flux. For our original benchmark, GRB 1, Bustamante et al. (2015) found that the flux is robust against variations in burst parameters like δt eng (see Table 3) and N sh . We will discuss below how it depends on underlying assumptions, and what the corresponding fluxes are for GRBs 2-6. Figure 7 shows the fluxes for GRBs 1-6. For all but GRB 5, the flux is E 2 J νµ ≈ 2 · 10 −11 GeV cm −2 s −1 sr −1 around 1 PeV -this is close to the value found for GRB Figure 7. All-sky quasi-diffuse νµ +νµ fluxes in our simulated multi-zone GRBs 1-6. Numerical one-zone predictions ) are included for comparison; they are calculated using the average burst parameters computed as in Section A.4. The shaded regions give the potential contribution from sub-photospheric collisions. The dominant contributions from individual collisions are shown as thin curves, corresponding to cases where the optical depth to photohadronic interactions τpγ is larger (red/light) or smaller (blue/dark) than unity. The IceCube 2016 upper limit was calculated using their latest reported detector effective area and exposure in a stacked GRB search using tracks coming from the Northern Hemisphere (Aartsen et al. 2016). Bustamante et al. (2015). For GRB 5, the flux is somewhat lower, as it has fewer optically thick collisions, which is in agreement with Eq. (2). (For GRB 3, the same could be expected, but instead it has a higher neutrino flux, as explained in Section 5.1.) Most real light curves lack the non-trivial features seen in GRB 5 (and GRB 3), and are therefore likely strong neutrino emitters. So it is conceivable that most bursts are instead like GRBs 1, 2, 4, and 6, and that the quasi-diffuse flux lies indeed at the level predicted using those bursts.

in
IceCube has searched for correlations between neutrino arrival directions and positions of known GRBs (Abbasi et al. 2012;Aartsen et al. 2015aAartsen et al. , 2016. No significant signal from GRBs has been found, in consistency with our prediction. Fig. 7 includes the differential upper bound from Abbasi et al. (2012).
One-zone and multi-zone models have similar (average) burst parameters and compute the flux of secondaries similarly. However, Fig. 7 shows that fluxes calculated with the multi-zone model are typically lower than with the one-zone model ) (see dashed curves in Fig. 7). The reason is that all shells are assumed to have the same collision radius in the one-zone model, which tends to be underestimated: since the neutrino production efficiency decreases non-linearly with the collision radius, the average value of the collision radius is, in general, not representative for the neutrino production. Nevertheless, we could define a single effective collision radius for neutrino production in the one-zone model that is different from radius for gamma rays; this would be done by folding in the production efficiency calculated with the fraction of collisions occurring at that radius. The average or representative jet parameterssuch as the typical collision radius -are derived from gamma-ray observations (see Section A.4), which implies that parameters representative for gamma rays may not be representative for the other messengers.

Cosmic rays
From Table 2, we can see that all of the GRBs 1-6 are relatively efficient cosmic-ray emitters, although the required energy output per GRB, of at least 10 53 erg in the discussed energy range, should likely be a factor of a few larger to explain UHECR observations. 4 Within the presented model, it is conceivable that GRBs are the sources of UHECRs.
The connection between cosmic rays and neutrinos depends on how UHECRs escape the shells. Photohadronic interactions will transform protons into neutrons; neutrinos will also be produced. If all cosmic rays escape as neutrons ("neutron escape") the connection is strong (Mannheim et al. 2001;Ahlers et al. 2011): one neutrino of each flavor is expected per observed UHECR proton. This possibility can be clearly ruled out (Abbasi et al. 2012;Baerwald et al. 2015). However, at the highest energies, when the proton Larmor radius exceeds the shell width, protons can directly escape the shells without interacting in it, which leads to a hard spectrum ("direct escape"). In addition, diffusion may lead to escape depending on properties of magnetic fields.
In each merged shell, one or another escape component dominates 5 , depending on the properties of the shell (Baerwald et al. 2013). If we consider escape processes other than neutron escape, which is implicitly assumed in most of the previous literature (Waxman & Bahcall 1997;Dermer & Atoyan 2003;Guetta et al. 2004;, the latest IceCube data cannot exclude GRBs as the sources of UHECRs even in a one-zone model, but constraints on the average shell parameters can be obtained from the efficient neutrino production . Figure 6, bottom, shows the maximum proton energy E p,max to which protons are accelerated in collisions in GRB 2, as a function of collision radius. Below the photosphere (black boxes), proton synchrotron losses dominate and E p,max increases with R C . Around R C 10 8.5 -10 10 km, protons reach 10 10 GeV and higher. This is where UHECRs are emitted. At large R C , falling magnetic fields yield lower acceleration rates and energies. Neutron emission is correlated with efficient neutrino production, since neutrons and charged pions are produced together in pγ interactions. However, this occurs only in a few collisions, in a narrow range of low collision radii, where proton and photon densities are high; see red dots. In effect, cosmic ray emission via neutron escape is limited by Eq. (2). Most collisions occur at larger radii, so that the average collision radius for CR emission tends to be higher than for neutrinos (blue circles). There, particle densities are low enough for direct proton escape to dominate, without associated neutrino production. Given that only few collisions are neutron-dominated, the pure neutron escape assumptions for GRBs (Ahlers et al. 2011) cannot be justified. However, a quantitative statement requires further study beyond the scope of this work, as it depends on the relative contribution between neutron-dominated and direct escape-dominated shells. For a discussion of UHECR nuclei, see Bustamante et al. (2015); Globus et al. (2015). Figure 8 shows a key feature of the multi-zone GRB model that is not captured by the one-zone model: that neutrinos, gamma rays, and UHECR protons are emitted from different regions of the jet . This holds regardless of the difference in burst parameters among GRBs 1-6. Neutrinos are produced close to the photosphere, as discussed above. UHECRs tend to be produced at somewhat larger radii. At low radii, UHECRs escape as neutrons; at larger radii, most UHECRs escape directly as protons. Gamma rays tend to come from even larger radii. While their production is more evenly distributed in collision radius, at low radii, pair-production (γγ → e + e − ) drives their energy down, so high-energy gamma rays mainly come from large radii.

DELAYED HIGH-ENERGY GAMMA RAYS
The detection of time delays between gamma-ray signals in different energy bands can provide insight into the dynamics of the GRB central engine and jet.
The maximum energy with which a photon can escape the shell where it is created is limited by the optical depth τ γγ to pair production via γγ → e + e − . A photon of a certain energy escapes only if τ γγ < 1. Close to the central engine, photon density and, therefore, optical depth, are high. Bustamante et al. (2015) showed that only low-energy gamma rays escape from that region. Higher-energy gamma rays escape at larger radii. In this section, we explore whether the different shell opacities lead to time delays between gamma-ray energy bands. Figure 9 -top row, left panel -shows that, in GRB 1, t obs is quite uniformly distributed in R C , especially between 10 8 and 10 10 km, i.e., above the photosphere. The central panel shows that for many collisions in this range the maximum gamma-ray energy is limited by pair production, while, from 10 9 km up, an increasing number of collisions is unaffected by it. However, the right panel shows that no separation exists between the arrival times of gamma rays limited and not limited by pair production. In other words, at any time during the burst, low-and high-energy gamma rays arrive indistinctly at Earth from everywhere inside the jet.
For GRB 5 (Fig. 9, bottom row), the situation is different. The average t obs increases with R C between 10 9 and 10 11 km, i.e., some early ( 20 s) gamma-ray detections come from mid-range radii, while all late detections tend come from large radii. As a result, some early gamma rays have lower energies, in the upper Fermi-LAT and lower CTA bands, limited by pair production. Later detections, coming from larger radii, will have consistently higher energies, not limited by pair production. There is also an impact on the neutrino light curve: Fig. 4 shows that the neutrino flux is much lower for later collisions, which come from larger radii, where neutrino production is inefficient. This behavior is characteristic of bursts with narrow Γ distribution, where collisions tend to occur at large radii and late in the jet evolution. Figure 10 shows the gamma-ray light curves in different energy bands for GRBs 3 and 5, our two simulations with narrow Γ distributions. To produce them, we have assumed that the power-law photon density in the source extends to high enough energies 6 . GRB 3 has a delayed start of a few seconds in the LAT band compared to the first detected peak in GBM, and of ∼ 10 s in the higherenergy CTA band. These delays depend strongly on the energy threshold of the instrument. In GRB 5, the LAT peak follows the GBM peak after ∼ 2 s, and the signal in CTA grows to comparable levels ∼ 10 s later. For this GRB, the suppression affects mainly the first peak in the light curve (the overall normalization of the light curves is arbitrary, but the relative normalization among the different energy bands is fixed). The early suppression of high-energy emission is consistent with observations; Figure 8.
Energy output as a function of collision radius in neutrinos, UHECR protons (Ep > 10 10 GeV), and gamma rays. The approximate photospheric and (assumed) circumburst radius are marked, as well as the UHECR escape regions where either neutron escape or direct escape dominates.
Since bursts with time delays between energy bands have narrow Γ distributions, they are associated to light curves with broad pulses overlaid with fast variability and possibly weak neutrino emitters (see Section 4). It is possible, in principle, to use the observation of delays in population studies to find how common these GRBs actually are, which affects the neutrino flux from the full GRB population.
To summarize, our predictions are: 1. In GRBs with light curves that have broader pulses overlaid with fast variability, time delays in different wavelength bands are possible.
3. If such delays are observed, the GRB can be a weak neutrino emitter.
4. The fundamental reason for these apparent delays is an early suppression -rather than an actual delay -of high-energy gamma rays coming from smaller R C , where the radiation densities are higher and the gamma rays cannot escape. This ef-fect is only observable if observation time and R C in the collisions are correlated.
An example of a GRB that could match these predictions is GRB 080916C (Abdo et al. 2009). These predictions are a qualitative summary based on examples GRB 1-6, which we believe to be representative of a larger set of examples that we have produced. Quantitative results depend on the chosen parameter values, which reflects the observation that no two light curves are identical. Fig. 9 captures the central feature we observe in all examples that exhibit a lag: for these, there is a correlation between t obs and R C (within the range of R C values in which the maximal photon energy depends on R C ). This correlation can be traced back to the engine emitting shells in a relatively narrow range of values of Γ. The lags in Fig. 10 are clearly visible, although one can see some differences concerning their interpretation. For example, for GRB 3, the first tall peak in the CTA band is clearly delayed by several seconds with respect to the Fermi-GBM band, whereas, for GRB 5, the precise size of the lag depends on the instrument threshold of CTA.

SUMMARY AND CONCLUSIONS
We have simulated the evolution of a baryon-rich GRB jet in the internal shock scenario of the fireball model, by keeping track of all relativistic plasma shells that propagate in it, of the collisions between shells, and of the gamma rays, protons, and neutrinos emitted at the shocks produced during the collisions. Unlike traditional one-zone models that extrapolate the behavior of the whole burst from a single representative collision, our multi-zone simulations consider many such collisions, each happening under different physical conditions.
We have demonstrated that it is possible to generate the various features observed in GRB light curves by varying the behavior of the central emitter, in particular, the initial speeds with which it puts out shells (see Section 4). The initial speeds determine where collisions between shells will occur during the jet evolution. In this approach, one can relate the properties of the central emitter to observables.
If the initial distribution of shell speeds is "disciplined" or narrow -even if the average speed changes with time -the associated gamma-ray light curve will have one or more broad pulses overlaid with fast time variability. The associated neutrino flux can be low, because collisions tend to occur at large collision radii. In addition, there is a correlation between observation time and collision radius, which implies that early-time collisions occur at low radii, where radiation densities are high, and high-energy gamma-ray signals are suppressed. As a consequence, we expect delays in the light curves in the Fermi-LAT energy band with respect to the ones in the Fermi-GBM band of a few seconds, and in the CTA band with respect to the Fermi-GBM band of order ten seconds.
If the distribution of speeds is broader, collisions occur over a wide range of collision radii and the light curve is dominated by fast time variability. In this case, neutrino production is always efficient, because, typically, several collisions occur where the radiation densities are high. In this case, we do not expect delays between gammaray wavelength bands, because there is not a strong correlation between observation time and collision radius. Inspection of many GRB light curves reveals that most are actually simple, while the ones typically presented in the literature tend to be the more interesting cases with non-trivial features. This means the class of bursts with broad speed distributions may be more representative of the "typical" GRB.
We have also qualitatively confirmed the findings from Bustamante et al. (2015) for very different parameter sets. Notably, we have shown that, regardless of the initial speed distribution of the shells, different messengers come from different regions of the same GRB: neutrinos predominantly come from regions close to the center, UHECR protons come from intermediate regions, and high-energy gamma rays tend to come from regions further out from the center, where photon densities are low enough that their energies are not limited by pair production on source photons. We have also confirmed the minimal predicted neutrino flux around ∼ 2 · 10 −11 GeV cm −2 s −1 sr −1 per flavor around 1 PeV, as long as a significant fraction of the GRBs has broad initial shell speed distributions, which explain observations better; see, e.g., our examples GRB 1, 2, 4 and 6.
Our main results and conclusions are based exclusively on collisions that occur above the photosphere, where photons are not trapped by Thomson scattering. This allows us to assume that the photon spectra in the shells have the same form as the observed spectra at Earth.
Our results could be exploited in targeted GRB neutrino searches, such as the ones performed by IceCube, to cull a smaller catalog of GRBs that are specially likely to be strong neutrino sources. The non-observation of neutrinos from bursts in such a catalog -where associated backgrounds are smaller -could result in stronger upper bounds on prompt high-energy GRB neutrino emission.
Our results could be also tested in different gamma-ray wavelength bands: we do not expect significant delays in GRBs with fast time variability without underlying pulse structure.
Due to their high luminosity, short duration, and the high angular precision with which they are detected, GRBs are arguably our best chance at finding a highenergy neutrino counterpart to electromagnetic emission. The study presented here is a step towards a necessary, realistic multi-messenger understanding of GRBs. The observation of these neutrinos would be a smokinggun signature of high baryonic loading, and thus of the paradigm that GRBs could be the sources of the UHE-CRs.
The authors thank John Beacom, Denise Boncioli, Valerie Connaughton,Željka Bošnjak, and Eli Waxman for useful discussion and comments on the manuscript. MB acknowledges the hospitality of DESY Zeuthen during the completion of this work. 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). MB acknowledges support from NSF Grant PHY-1404311. The work of KM is supported by NSF Grant PHY-1620777.

A.1. Overview
The central object in a GRB (e.g., the black hole created by a collapsing massive star) emits collimated, relativistic jets that are rich in baryons, which we assume to be protons. When one of the jets points towards Earth, the gamma-ray emission from it may be detected as a GRB. We simulate the particle emission from this jet.
Since we are located inside the opening angle of the jet, we cannot distinguish between the emission geometry of a collimated jet and that of spherical, isotropic emission. For simplicity, we develop our formalism under such an isotropically-equivalent scenario. In it, the central engine emits spherical plasma shells that propagate outwards along the jet at relativistic speeds. We simulate their propagation and collisions between them, which produce high-energy gamma rays, protons, and neutrinos.
The simulation covers only the coasting phase of the GRB, during which shells propagate at constant speed, except when they collide. In the preceding, acceleration phase, they reached their maximum individual speeds, limited by the available kinetic energy and their masses. The coasting phase ends when the shells reach the circumburst medium; there, they decelerate, and might produce an afterglow. The acceleration and deceleration phases are not part of the simulation.
In our simulation, shells propagate in one dimension. At any time during the simulation, the k-th shell is characterized by four basic parameters: r k , the shell radius, as measured from the emitter (i.e., the position of the shell inside the jet); l k , the shell width; Γ k , the shell bulk Lorentz factor; and m k , the shell mass.
When two shells collide, they merge into a new shell, with width, speed, and mass calculated from the properties of the shells that collided. The new shell continues propagating in the jet flow and may collide again. Collisions are inelastic. The new shell cools instantly by radiating away its internal energy via particle emission. We compute collisions numerically, following Kobayashi et al. (1997); Daigne & Mochkovitch (1998), as detailed below. Table 3 describes all the relevant simulation parameters. Unless otherwise noted, quantities therein and in the text are expressed in the source reference frame.

A.2. Burst initialization
In the simulation, before shell propagation starts, the central engine has already emitted N sh shells. Each one is described by the initial tuple r k,0 , l k,0 , Γ k,0 , E iso kin,0 . Shells closer to the engine are labeled with higher indices. The behavior of the engine is described by two timescales 7 : an "uptime", δt eng , during which it emits one shell, followed by a "downtime", ∆t eng , during which it is inactive. The former determines the initial shell width, l = c·δt eng , which we assume to be common for all shells, and the latter determines the initial separation between consecutive shells, d = c · ∆t eng . Thus, each of the initial shells is located at position r k,0 = r N sh + (N sh − k) (l + d), where r N sh is the distance from the innermost shell to the emitter, which is an input parameter of the simulation. Results do not depend on r N sh strongly, unless its value is too large; see Table 3.
We choose the values of l and d to reproduce the timescale of pulses in observed light curves (Nakar & Piran 2002a). If t v is the GRB variability timescale, i.e., the characteristic duration of peaks in the light curve, and t q is the characteristic quiescent time between consecutive peaks, we expect that, roughly, l ≈ c · t v and d ≈ c · t q (Kobayashi et al. 1997;Aoi et al. 2010). In the internal shock model, d and l should be comparable. The simulations in Kobayashi et al. (1997) set d = l, while Aoi et al. (2010) set d = 5l. In our simulations, we chose δt eng = ∆t eng , such that d = l = c · δt eng .
The variability timescale of a simulated burst is not an input parameter of the simulation, but a result of it. For our choices of simulation parameter values, we find that the variability timescale, obtained from the post-simulation synthetic light curve (see Section A.4), is close to the input value of δt eng .
The initial values of the shell Lorentz factors follow a pre-defined distribution. In the benchmark scenario GRB 1, it is a log-normal distribution; see Section 4. Table 1 describes the distributions used in GRBs 2-6.
There are two typical schemes to assign initial masses m k,0 to the shells: the equal-mass assumption, i.e., m k,0 = m for all k; and the equal-energy assumption, i.e., m k,0 = E iso kin,0 / Γ k,0 c 2 , with E iso kin,0 the initial bulk kinetic energy, assumed common to all shells. Our simulation uses the latter, since it appears to match observations more closely (Nakar & Piran 2002a).

A.3. Burst evolution
We simulate the coasting phase of the jet, during which shell speeds do not change while they propagate and expand. In our simplified treatment, the shell width and mass also stay constant 8 . Therefore, the shell volume V iso,k = 4πr 2 k l k Initial radius of the k-th shell, Internal km Initial bulk Lorentz factor of the k-th shell Internal -Sampled from the pre-defined Γ distribution m k,0 Initial mass of the k-th shell Internal GeV Note. -All quantities are expressed in the source (engine) frame, except for t obs , which is in the observer's frame.
grows ∝ r 2 k and mass density ρ k = m k /V iso,k decreases ∝ r −2 k . Since the shell mass and Lorentz factor are constant during propagation, its bulk kinetic energy E iso kin,k = Γ k m k c 2 is constant as well. Speed, width, and mass change only in collisions.
At the start of the simulation, we calculate the collision time for all pairs of neighboring shells, i.e., where d k,k+1 ≡ r k − r k+1 − l k+1 is the separation between shells k and k + 1. The time interval until the next collision occurs is the minimum of these times, i.e., We increase the simulation time to t → t + ∆t next . The collision radius R C is set to the radius of the innermost colliding shell. Light emitted from this collision will be detected by a distant observer at time with D (z) the light-travel distance to the emitter with redshift z. These equations satisfy well-known relations in the internal shock scenario, ∆t next ≈ 2Γ 2 (d/c) and R C ≈ 2Γ 2 d. All shells propagate to their new positions r k → r k + cβ k ∆t next , the time interval for the next collision is calculated, and the process is repeated. (The term D (z) /c is just an offset: it will disappear when, in the simulation output, the first emission is set to start at t obs = 0.) We distinguish between the shell bulk kinetic energy, E iso kin,k , and its internal energy, E iso int,k . The former is related to the motion of the compact shell, measured in the source rest frame, while the latter is the aggregated kinetic energy Figure 12. Time evolution, in the source frame, of average shell mass m / m 0 , standard deviation of the Lorentz factor σ(Γ)/σ(Γ 0 ), and total internal energy of a burst, E iso int,tot /E iso int,tot,0 . The ranges are from our numerical results from a 1000 simulations run with random setups for the parameter values N sh = 10000, Γ 0 = 100, A Γ = 0.2, δteng = 10 −3 s, d = l, z = 2, and E iso kin,0 = 10 52 erg. Solid and dashed lines come from Fig. 1 in Beloborodov (2000) for the same parameter set and refer to numerical calculations and analytical estimates respectively.
factor 9 σ(Γ)/σ(Γ 0 ), and total available internal energy of the burst E iso int,tot /E iso int,tot,0 . The latter is calculated directly as E iso int,tot = Γ CM m i (β 2 i /2), whereβ is the speed in the CM-frame. 10 The numerical results of our simulation match the analytical power-law estimates from Beloborodov (2000), which assume that fluctuations in the initial Lorentz factors are small, i.e., A Γ 1. Deviations occur at late times, when the number of remaining shells is low and the analytical predictions are no longer applicable. This late deviation depends strongly on the random initial setup, so we show ranges obtained after running 1000 different simulations.
Changing the collision dynamics can affect the burst efficiency -i.e., the fraction of kinetic energy dissipated as secondaries -and the radii where most of the gamma-ray energy is dissipated. Changes could include re-converting a fraction of collision energy into kinetic energy, partially inelastic collisions, or even fully penetrating shells. We explore simple modifications of our collision model in Appendix B.
A.4. Gamma-ray observables The internal energy of a merged shell is split among electrons, magnetic field, and protons. They receive a fraction e , B , and p , respectively. We assume energy equipartition between electrons and photons and fix e = B = 1/12 and p = 5/6, since this yields the frequently used value of baryonic loading 1/f e = p / e = 10 (Waxman & Bahcall 1999;Abbasi et al. 2012). Thus, the k-th collision dissipates an energy E γ,k = e E coll,k as gamma rays, and energy E γ,p = p E coll,k as protons, and supports a magnetic energy density of U B = B E coll,k /V iso . The latter translates into a magnetic field intensity, in the shock rest frame, of We normalize these individual energies by requiring that the total dissipated energy as gamma rays, an experimentally accessible quantity, matches a given value E iso γ,norm = 10 53 erg. This also fixes the energy in protons and magnetic field. Our simulation does not generate the photon spectrum in the shell. Instead, as in Aoi et al. (2010); Abbasi et al. (2012), we assume that its shape is that of observed GRB spectra. We parametrize the spectrum as a broken power of a gas, E iso int,tot , and its volume is given by E iso int,tot ∝ V −(γ−1) iso ∝ r −2(γ−1) , where γ is the adiabatic index. For the relativistic gas in a shell, γ = 4/3 and E iso int,tot ∝ r −2/3 ∝ t −2/3 . law, i.e., Primed quantities are in the shock rest frame. We fix 11 α γ = 1, β γ = 2, and ε γ,break = 1 keV. The photon spectrum in each shell is normalized via where the minimum and maximum energies are ε γ,min = 0.2 eV and ε γ,max = 1 PeV, respectively (Aoi et al. 2010). Pair production via γγ → e + e − may limit the maximum energy of escaping photons; see Fig. 9, right column. Each collision emits a gamma-ray pulse. The superposition of all pulses, propagated to Earth, is the light curve of the burst; see Section 4. Following Kobayashi et al. (1997), we parametrize the luminosity of the pulse from the k-th collision (in the observer's frame) as a peaked profile, with a fast rise and exponential decay ("FRED"), i.e., where the emission timescale, i.e., the time at which the reverse shock crosses the fast shell, is and the "rise time", is the time elapsed since the start of the emission until the peak luminosity is reached. For an illustration, see Fig. 1 in Kobayashi et al. (1997). The peak luminosity is The time t in the source frame is related to t obs through t = 2Γ 2 m,k t obs 1 + z .
Hence, the synthetic light curve L γ is Unless noted otherwise, we only show the light curves for collisions beyond the photosphere (see Section 4), while we use all collisions for the normalization in Eq. (A11). For the examples we show, the fraction of energy dissipated below the photosphere is around 50%, so that the total super-photospheric energy output in gamma rays is around half of E iso γ,norm . In Section 6 we showed light curves in different energy bands. Each band spans the range E band γ,low , E band γ,high , with E band γ,low = 10 −6 , 10 −1 , 10 2 GeV, and E band γ,high = 10 −2 , 10 2 , 10 6 GeV, for Fermi-GBM, Fermi-LAT, and CTA, respectively.
11 It is uncertain how ε γ,break changes with R C , since the scaling expected in the internal shock model has a problem (Daigne & Mochkovitch 2003), which is why we prefer to set the photon spectra to observations. However, depending on the origin of the prompt gamma-ray emission -e.g., synchrotron radiation, inverse Compton scattering -one can implement model-specific assumptions, as in ;    Bošnjak et al. (2008). Detailed information on mildly relativistic collisionless physics, such as the injection Lorentz factor (that depends on the number fraction of accelerated particles (Eichler & Waxman 2005)), is also necessary. However, such an improvement will not change our conclusions in Section 7. Note.
Following Baerwald et al. (2013), UHECRs escape from each shell as neutrons, produced in photohadronic interactions ("neutron escape"), and as protons that leak out of the shell when their Larmor radius exceeds the shell thickness ("direct escape"). At the highest energies, protons can always leak out -provided their maximum energy is limited by adiabatic cooling. However, which escape component dominates in each shell depends on the optical depth to photohadronic interactions of the shell in question. Fig. 13 is a schematic illustration of the components contributing to UHECR emission depending on the optical depth of the merged shell. (If magnetic fields decay fast enough in the bulk of the shell, the direct escape component may be larger. However, this possibility is not contemplated in our framework.) Our main results, e.g., our light curves and neutrino spectra, are based only on shell collisions that occurred above the photosphere. Below it, Thomson scattering off electrons keeps the photons trapped in the shell. Since our adopted photon spectra are chosen to reproduce observed gamma-ray spectra, we cannot accurately calculate secondary production below the photosphere, where the photon spectra might be different. We mark "sub-photospheric" collisions clearly as such (see, e.g., Figs. 5 and 6) and exclude them from our flux calculations. Excluding sub-photospheric collisions does not qualitatively change the shape of the light curves. However, in cases with broader pulse structures, the onset of each pulse is usually dominated by sub-photospheric collisions. Excluding these collisions slightly delays the onset of each peak.
The optical depth to Thomson scattering is calculate from shell properties. Since shells are, on average, electrically neutral, the electron density is equal to the proton density, i.e., where m p is the proton mass. This assumes that electron-positron pair production does not increase the electron and positron densities significantly. The optical depth to Thomson scattering is then τ Th,k 1 n e,k σ Th l k , with σ Th ≈ 66.52 fm 2 the Thomson cross section. A collision is sub-photospheric if τ Th > 1. For a burst with E iso γ,tot 10 53 erg, e = 1/12, and a dissipation efficiency of ε = 25% (such as GRB 1), the initial kinetic energy per shell is about 10 51.6 erg if 1000 collisions occur. This yields m k,0 = E iso kin,0 /Γ k,0 c 2 10 49 erg for Γ k,0 ∼ 500. From Eqs. (A28) and (A29), the photospheric radius of the k-th shell is R ph,k 1.8 · 10 8 km m k 10 49 erg 1/2 .
(A30) Bustamante et al. (2015) showed that, since the photohadronic optical depth scales similarly to Eq. (A29) (replacing σ Th → σ pγ ), the pion production efficiency at the photosphere is independent of the isotropic volume and only weakly dependent on Γ k . The total neutrino flux of the burst is dominated by a few collisions that occur just above the photosphere, where the pion production efficiency is highest.

B. ADDITIONAL SIMULATIONS
Here we include four additional examples of GRB simulations, GRBs 7-10, that complement the ones showed in the main text. Table 4 describes the simulations. Figure 14 shows the corresponding light curves.
The light curves of GRB 7 and GRB 8 are similar to that of GRB 1; see Fig. 4. The light curves of GRB 9 and GRB 10 are similar to that of GRB 4 and GRB 6.
The similarities in the light curves exist in spite of fundamental differences between the behavior of the engine. This illustrates our point that the qualitative behavior of the examples shown in this work are representative of a larger class of models. Figure 14. Synthetic gamma-ray and neutrino light curves for the simulated GRBs 7-10, from collisions beyond the photosphere. Photon counts are in arbitrary units, obtained by multiplying the flux times a factor of 10 6 GeV −1 cm 2 s. Figure 15. Three modifications to the canonical collision model described in Appendix A and applied to GRB 1. The modified case with no multiplication collisions (green dotted) uses the same per-collision normalization as the reference case (orange, solid). The other two modifications -50% of internal energy reconverted to kinetic energy (blue, dot-dashed) and 50% radiated with 50% reconverted (red, dashed) -are normalized so in each case the burst yields 10 53 erg in gamma rays when adding sub-photospheric and super-photospheric collisions.

C. ALTERNATIVE COLLISION DYNAMICS
Here we discuss the impact of modifications to our canonical collision model, which is used in the main text and described in Appendix A. We focus on alternative scenarios that can be easily implemented in our framework; that is, we assume that, in each collision, the colliding shells merge and do not consider the case in which they reflect off each other, as in Kobayashi & Sari (2001).
One extreme modification is to remove colliding shells from the system after they collide and radiate, which means that multiple collisions are not allowed. This makes simulation results insensitive to details of how shells are treated after colliding. However, removing the shells will modify the whole system, since, in the canonical collision model, collisions among old shells, and between young and old shells, occurred relatively early on. Figure 15 shows the effect on the quasi-diffuse neutrino flux of applying this modification to a simulation that has the same parameters as our reference case, GRB 1. The modified case is labeled "no multiple collisions". The collision energies were normalized using E iso γ,norm = 10 53 erg. The number of collisions is reduced to 499, about half that of the reference GRB 1 model. The modification results in a higher neutrino flux, because, by forbidding multiple-time collisions, most collisions -all of them first-time -occur at low radii, around 10 8.5 km. As we normalize to the same total energy, the per-collision normalization is slightly higher, which further increases the neutrino flux.
Another modification is to assume that only a fraction η of the internal energy in each collision, Eq. (A4), is attributed to the non-thermal spectra of the secondaries, while 1 − η is instantly reconverted into kinetic energy of the merged shell. This fraction 1 − η can, for instance, describe a fraction of thermal protons not directly participating in the prompt emission. For simplicity, we assume that the extra kinetic energy translates into an instantaneous increase of the Lorentz factor the merged shell after cooling: Γ m = [(1 − η)E iso coll + E iso kin,m ]/m m . Yet another modification is linked to our assumption that protons can only directly escape the merged shells -and not leave them by diffusion (Baerwald et al. 2013) or other processes -which implies that a substantial fraction of the non-thermal baryonic energy will remain trapped by the magnetic fields and eventually reconverted into kinetic energy. The typical fraction ξ of electromagnetic energy and non-thermal baryonic energy which is actually radiated for this escape process is 40-50%, estimated from energy partition and from the proportion of baryonic energy in the UHECR energy range compared to the full energy range. To account for this, we consider the amount of internal energy used for computing the secondary production is still given by Eq. (A4), but we modify the dynamics so that a fraction of internal energy is reconverted into kinetic energy. This case does not include a fraction of energy going into thermal protons, unlike the previous case. Figure 15 shows the result of both of these modifications to the kinetic and radiated energy, for the case η = ξ = 0.5. In both cases, the burst was normalized to yield 10 53 erg in gamma rays when adding sub-photospheric and superphotospheric collisions. The number of collisions is similar to that of the reference GRB 1 model. Because the Lorentz factors of the merged shells are higher due to the increased kinetic energy, collisions occur further out in the jet, where particle densities are lower and neutrino production is less efficient. As a result, the neutrino flux associated to these two modifications is slightly lower than the one associated to the reference case, especially if only a fraction η of the energy is radiated. Therefore, the minimal super-photospheric flux prediction of ∼ 10 −11 GeV cm −2 s −1 s −1 holds.
The neutrino flux scales with the fraction η going into the non-thermal spectra (and magnetic field), which means that, for η = 0.1, it would be about one order of magnitude lower than our nominal case. On the other hand, the result is rather insensitive to the fraction of reconverted non-thermal energy 1 − ξ. This means that for the combined case -a fraction η into non-thermal spectra and a fraction 1 − ξ of non-thermal energy reconverted -we expect that the result is dominated by the effect of η.
Finally, the three modifications to the collision dynamics that we have explored do not affect our conclusion about the distribution of particle emission with collision radii: neutrinos still come from low radii, UHECR protons come from intermediate radii, and gamma rays come from large radii.