Multi-messenger Light Curves from Gamma-Ray Bursts in the Internal Shock Model

The most energetic particles discovered—ultra-high-energy cosmic rays (UHECRs), above 10⁹ GeV, and high-energy 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. Gamma-ray 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 (Paczynski & Xu 1994; Rees & Mészáros 1994). If the jets contain enough baryons, then protons (Milgrom & Usov 1995; Vietri 1995; Waxman 1995) and nuclei (Murase et al. 2008; 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. 2015a, 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 because 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 gamma-ray 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.

GRBs are conceivably fueled by matter accretion, either in the collapse of a massive star—long-duration bursts, lasting $\gt 2$ s—or in the merging of two neutron stars or a neutron star and a black hole—short-duration bursts, lasting $\lt 2\,{\rm{s}}$. The central engine, likely a newly formed black hole, emits two relativistic matter jets in opposite directions. When one of them points toward Earth, we might detect it as a GRB. We focus on long-duration bursts because they have higher gamma-ray luminosity ($\sim {10}^{52}$ erg s⁻¹), represent $\sim 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.

The GRB prompt emission mechanism has been under debate for years (see reviews by Piran 1999; Mészáros 2006; Kumar & Zhang 2014; Mészáros et al. 2015; Pe'er 2015). In the classical model, prompt gamma-ray emission in the MeV range is attributed to optically thin synchrotron emission from non-thermal electrons accelerated at internal shocks (Rees & Mészáros 1994), possibly supplemented by a thermal component from the fireball. Indeed, there are observational indications of such a thermal component (see, e.g., Pe'er 2012; Guiriec et al. 2011). However, this classical internal shock model has difficulty explaining the data, including the low-energy spectral index (Preece et al. 1998), radiation efficiency (Kumar 1999; Zhang et al. 2007), and spectral energy relations (Amati et al. 2002; Yonetoku et al. 2004).

Among alternative theories, photospheric emission models became 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).

Magnetic reconnection models have also been largely of interest (Mészáros & Rees 1997; Lyutikov & Blandford 2003; Zhang & Yan 2011; Bošnjak & Kumar 2012; McKinney & Uzdensky 2012), due to certain advantages (Kumar & Zhang 2014), though detailed microphysics is much more uncertain. For example, they can explain the high polarization observed in some GRBs (Gotz et al. 2009; Yonetoku et al. 2012) and the absence of bright thermal emission (Zhang & Pe'er 2009; Gao & Zhang 2015).

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}_{{\rm{v}}}$, of ∼10–100 ms  (Bhat 2013; Golkhou et al. 2015), frequently superposed on top of slower pulse structure (Ramirez-Ruiz & Fenimore 1999; Nakar & Piran 2002b). Since the light curve reflects the particularities of the central emitter and jet propagation (Nakar & Piran 2002a; Lopez-Camara et al. 2016; Zhang 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 $\langle {\rm{\Gamma }}\rangle$—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}_{{\rm{v}}}\sim 1000$, with ${T}_{90}\sim 10\,{\rm{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⁷ of Lorentz factors or emission radii have been addressed qualitatively (Murase & Nagataki 2006) and quantitatively (Guetta et al. 2001a, 2001b; 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 $\sim {10}^{-11}\,\mathrm{GeV}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{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. 2013a, 2013b, 2014a, 2014b, 2015b, 2015c). However, no individual neutrinos are associated to known GRBs (Abbasi et al. 2012; Aartsen et al. 2015a, 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 (Laha et al. 2013; Murase & Ioka 2013; Bustamante et al. 2015). (However, low-luminosity classes of GRBs can describe the diffuse neutrino flux without violating the stacking limits (Murase et al. 2006; 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; He et al. 2012; Hümmer et al. 2012; Li 2012). While the plain one-zone ansatz with commonly assumed parameters with values that are favorable for the UHECR explanation, is starting to be challenged by IceCube data (He et al. 2012; Hümmer 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.

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) and 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}_{\mathrm{sh}}\sim 1000$ shells. Each one propagates at a different relativistic speed; typically, they have Lorentz factors of ${\rm{\Gamma }}\sim 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.

Zoom In Zoom Out Reset image size

Non-thermal electrons receive a fraction ${\epsilon }_{e}$ of the internal energy of the new shell, protons receive a fraction ${\epsilon }_{p}$, and magnetic fields receive a fraction ${\epsilon }_{B}$. Because electrons cool quickly via synchrotron emission, ${\epsilon }_{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 10 times more energy in protons (Abbasi et al. 2010; Hümmer et al. 2012), i.e., ${\epsilon }_{e}={\epsilon }_{B}=1/12$ and ${\epsilon }_{p}/{\epsilon }_{e}=10;$ see Appendix A.4. Thanks to the fast cooling, ${\epsilon }_{p}/{\epsilon }_{e}$ can be regarded as the non-thermal baryon loading factor, defined as the ratio of cosmic-ray energy to radiation energy (Murase & Nagataki 2006).

The volume of the shell (see Appendix A.3) grows with the distance r to the central emitter ∝ r² 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. Thus, as a shell propagates and expands, the density of particles in it falls ∝r⁻². 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; ultra-high-energy cosmic rays come from intermediate collision radii; and high-energy gamma-rays escape abundantly from large collision radii (Bustamante et al. 2015).

In a shell collision, if both protons and electrons are accelerated up to sufficiently high energies, $p\gamma$ interactions create pions via the ${{\rm{\Delta }}}^{+}(1232)$ resonance and other processes. High-energy gamma-rays are produced in ${\pi }^{0}\to \gamma \gamma$ (and by synchrotron and inverse-Compton radiation by electrons), while high-energy neutrinos are produced in ${\pi }^{+}\to {\mu }^{+}{\nu }_{\mu }\to {\bar{\nu }}_{\mu }{e}^{+}{\nu }_{e}{\nu }_{\mu }$ (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 $\sim 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 $\propto {E}_{p}^{\prime -2}\exp (-{E}_{p}^{\prime }/{E}_{p,\max }^{\prime })$, with ${E}_{p}^{\prime }$ the proton energy, as expected from Fermi acceleration (primed quantities are in the shock rest frame). The maximum proton energy ${E}_{p,\max }^{\prime }$ 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, Equation (14), 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 the slow heating of electrons (Bykov & Mészáros 1996; Murase et al. 2012; Asano & Terasawa 2015). In this sense, our approach is model-independent; 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 (Hümmer et al. 2012) software (see Appendix A.5).

UHECRs are emitted via two mechanisms: either $p\gamma$ 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 with neutron escape is higher than that associated with 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 Appendix 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; Murase 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 by, 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}_{\gamma ,\mathrm{norm}}^{\mathrm{iso}}={10}^{53}$ erg (in the source reference frame). This is equated to the sum energy emitted by all collisions, sub- and super-photospheric (Equation (13)). In our simulated bursts, about $50 \%$ of the energy is liberated in super-photospheric collisions.

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). The time resolution of the detector imposes a lower bound on the variability timescale that can be inferred. Pulses typically have a duration of ∼1 s with an asymmetric structure (Ramirez-Ruiz & Fenimore 1999; Nakar & Piran 2002b), 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; Beloborodov 2000; Spada et al. 2000; Kobayashi & Sari 2001; 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;⁸  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, ${{\rm{\Gamma }}}_{k,0}$, are randomly sampled from a log-normal distribution defined by the characteristic value ${{\rm{\Gamma }}}_{0}$ and the amplitude of fluctuations A_Γ:

$$\begin{eqnarray}&&\mathrm{ln}\left(\displaystyle \frac{{{\rm{\Gamma }}}_{k,0}-1}{{{\rm{\Gamma }}}_{0}-1}\right)={A}_{{\rm{\Gamma }}}\cdot x,\end{eqnarray} \tag{ 1 }$$

where the random variable x follows a Gaussian distribution, $P(x){dx}={(2\pi )}^{-1}{e}^{-{x}^{2}/2}{dx}$. When ${A}_{{\rm{\Gamma }}}\lt 1$, the mean value $\langle {\rm{\Gamma }}\rangle \approx {{\rm{\Gamma }}}_{0}$ and the variance ${\rm{\Delta }}{\rm{\Gamma }}\approx {A}_{{\rm{\Gamma }}}{{\rm{\Gamma }}}_{0}$. When ${A}_{{\rm{\Gamma }}}\gt 1$, the mean value and variance are significantly affected by fluctuations. Large ${A}_{{\rm{\Gamma }}}\gtrsim 1$ are typically required for the 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}_{{\rm{\Gamma }}}\gtrsim 0.1$ is necessary for the number of collisions to be large, i.e., ${N}_{\mathrm{coll}}\sim {N}_{\mathrm{sh}}$. The light curve for GRB 1 is in Figure 4. It is dominated by fast variability, which is determined by the ratio of T₉₀ to ${N}_{\mathrm{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).

Zoom In Zoom Out Reset image size

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

Table 1.  Description of Simulated GRBs 1–6

Model	${{\rm{\Gamma }}}_{\mathrm{0,1}}$	${A}_{{\rm{\Gamma }},1}$	${{\rm{\Gamma }}}_{\mathrm{0,2}}$	${A}_{{\rm{\Gamma }},2}$	Tp	${N}_{\mathrm{up}}$	${N}_{\mathrm{down}}$	${E}_{\gamma ,\mathrm{norm}}^{\mathrm{iso}}$ (erg)	Description
1	500	1.0	⋯	⋯	⋯	⋯	⋯	1053	Fixed Γ and AΓ; benchmark from Bustamante et al. (2015)
2	500	1.0	50	0.1	⋯	⋯	⋯	1053	Linear speedup of Γ
3	50	0.1	500	0.1	0.34	⋯	⋯	1053	Sawtooth Γ (linear slowdown three times) with narrow distribution
4	50	0.1	500	1.0	0.2	⋯	⋯	1053	Oscillating Γ (five periods) with increasing distribution width
5	50	0.1	500	0.1	0.2	⋯	⋯	1053	Oscillating Γ (four periods) with lower amplitude increasing and narrow distribution
6	50	0.1	500	1.0	1/8	250	250	1053	Oscillating Γ (four periods) with distribution widening up; in addition, engine intermittent: ${N}_{\mathrm{up}}$ pulses followed by ${N}_{\mathrm{down}}$ pulses; corresponds to increasing Lorentz factor during uptime

Note. Common values for all models: ${N}_{\mathrm{sh}}=1000$, $\delta {t}_{\mathrm{eng}}={10}^{-2}$ s, $d=l=c\cdot \delta {t}_{\mathrm{eng}}$, ${r}_{\min }={10}^{3}$ km, ${r}_{\mathrm{dec}}=5.5\cdot {10}^{11}$ km, z = 2, ${\epsilon }_{e}={\epsilon }_{B}=1/12$, ${\epsilon }_{p}=5/6$, and $\eta =1.0$ (acceleration efficiency; Baerwald et al. 2013). See Table 3 for an explanation of each parameter. The period T_p for the oscillating cases refers to Γ changing between ${{\rm{\Gamma }}}_{\mathrm{0,1}}$ and ${{\rm{\Gamma }}}_{\mathrm{0,2}}$, and A_Γ changing between ${A}_{{\rm{\Gamma }},1}$ and ${A}_{{\rm{\Gamma }},2};$ T_p is a fraction of the total number of emitted shells. This means that Γ and ${A}_{{\rm{\Gamma }}}$ change between first and second value with a factor ${\sin }^{2}(k/({N}_{\mathrm{sh}}\cdot {T}_{p})\cdot \pi )$, where k is the index of the shell ($1\leqslant k\leqslant {N}_{\mathrm{sh}}$). For GRB 5, the lower amplitude increases from ${{\rm{\Gamma }}}_{\mathrm{0,1}}$ to ${{\rm{\Gamma }}}_{\mathrm{0,2}}$ linearly with k.

Download table as:  ASCIITypeset image

Figure 3 shows the randomly sampled initial values of the Lorentz factors ${{\rm{\Gamma }}}_{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.

Zoom In Zoom Out Reset image size

Figure 4 shows that 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 because 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 Figure 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.

Zoom In Zoom Out Reset image size

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 Figure 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 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}_{{\rm{\Gamma }}}=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; Beloborodov 2000; Spada et al. 2000; Kobayashi & Sari 2001; 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 (Lazar et al. 2009; Narayan & Kumar 2009; Zhang & Zhang 2014; 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 (Zhang & Zhang 2014). Compared to these explanations, our pulse structure comes from the properties of the engine rather than the jet.

Figure 4 shows that 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 $\propto {R}_{{\rm{C}}}^{2}$, where ${R}_{{\rm{C}}}$ is the collision radius.

There are no long time delays between gamma-ray and neutrino pulses: they are within T₉₀ (see Equation (23)) 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 with one collision depends strongly on Γ (see Equation (18)) and, consequently, it is $\propto {R}_{{\rm{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 $\sim 50\,\mathrm{ms}$ and durations of $\sim 50\,{\rm{s}}$. In GRB 3, the duration ($33\,{\rm{s}}$) and time variability ($36\,\mathrm{ms}$) are smaller because mainly two of the three peaks contribute. In GRB 6, pulses are separated by a downtime and, therefore, duration and variability are larger ($98\,{\rm{s}}$ and $97\,\mathrm{ms}$, respectively).

Table 2.  Parameters Output by Simulated GRBs 1–6

Model	${N}_{\mathrm{coll}}$	${t}_{{\rm{v}}}$ (ms)	T90 (s)	${E}_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}$ (erg)	${E}_{p,\mathrm{tot}}^{\mathrm{iso}}$ (erg)	${E}_{\nu ,\mathrm{tot}}^{\mathrm{iso}}$ (erg)	${E}_{\nu ,\mathrm{tot}}^{\mathrm{iso}}/{E}_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}$ $( \% )$	ε $( \%$)
1	987	53	54.0	$5.2\cdot {10}^{52}$	$6.2\cdot {10}^{52}$	$1.4\cdot {10}^{52}$	26.9	26.8
2	999	47	47.0	$6.5\cdot {10}^{52}$	$4.4\cdot {10}^{52}$	$9.3\cdot {10}^{51}$	14.3	19.6
3	951	33	35.5	$6.2\cdot {10}^{52}$	$5.0\cdot {10}^{52}$	$7.6\cdot {10}^{51}$	12.3	10.5
4	987	52	52.6	$4.6\cdot {10}^{52}$	$4.0\cdot {10}^{52}$	$9.4\cdot {10}^{51}$	20.4	21.7
5	990	57	57.9	$8.9\cdot {10}^{52}$	$1.7\cdot {10}^{52}$	$6.6\cdot {10}^{50}$	0.7	10.6
6	985	97	98.0	$6.1\cdot {10}^{52}$	$4.2\cdot {10}^{52}$	$1.1\cdot {10}^{52}$	18.0	23.6

Note. The parameters are variability timescale (${t}_{{\rm{v}}}$), total energy emitted as gamma-rays (${E}_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}$), as protons (${E}_{p,\mathrm{tot}}^{\mathrm{iso}}$), and as neutrinos (${E}_{\nu ,\mathrm{tot}}^{\mathrm{iso}}$), ratio between neutrino and gamma-ray energies (${E}_{\nu ,\mathrm{tot}}^{\mathrm{iso}}/{E}_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}$), and overall emission efficiency ε. Energies are computed using super-photospheric collisions only. Only protons in the UHECR energy range, above $\gt {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.

Download table as:  ASCIITypeset image

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 simulations⁹ have high $\varepsilon \approx 11 \% \mbox{--}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 ${\varepsilon }_{\gamma }\approx {\epsilon }_{e}\varepsilon$ is about 10 times smaller since most of the dissipation energy is assumed to be carried by protons rather than electrons. Here, ${\epsilon }_{e}$ is the fraction of energy in electrons and photons; see Appendix 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).

Here we discuss the emission of multiple messengers and their relation.

5.1. Weak versus Strong Neutrino Emitters

The time-integrated neutrino fluence of a simulated burst, for a baryon-rich jet, roughly scales as (Bustamante et al. 2015)

$$\begin{eqnarray}&&{{ \mathcal F }}_{\nu }\propto \displaystyle \frac{{N}_{\mathrm{coll}}({\tau }_{p\gamma }\gtrsim 1)}{{N}_{\mathrm{coll}}}\cdot \varepsilon \cdot {E}_{\gamma ,\mathrm{tot}}^{\mathrm{iso}},\end{eqnarray} \tag{ 2 }$$

assuming a fixed photon break energy (see Appendix A.4). The first factor gives the fraction of collisions with high optical depth ${\tau }_{p\gamma }\gtrsim 1$ to $p\gamma$ interactions, ε is the energy dissipation efficiency, and ${E}_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}$ is the total energy emitted as gamma-rays in super-photospheric collisions. Unlike one-zone predictions (Waxman & Bahcall 1997; Guetta et al. 2004; He et al. 2012; Hümmer et al. 2012; Li 2012), the fluence does not depend on the average Lorentz factor of the shells.

Figure 5 shows ${\tau }_{p\gamma }$ as a function of ${R}_{{\rm{C}}}$ for collisions in GRBs 1 and 5. In GRB 1, a strong neutrino emitter, about 60 in 1000 collisions 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 Equation (2) is ∼0.05 for strong neutrino emitters. The energy dissipation efficiency in Table 2 lies around $\varepsilon =0.2$ for these GRBs. As a result, for fixed ${E}_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}$, the quasi-diffuse neutrino flux that we infer below (Bustamante et al. 2015) is relatively robust.

Zoom In Zoom Out Reset image size

The situation is different for GRBs 3 and 5. They have lower efficiencies $\varepsilon \approx 10 \%$. More importantly, they have no optically thick collisions close to the photosphere; see the bottom panel of Figure 5 for GRB 5. The reason is that they emit shells with a variable but narrow Lorentz factor distribution (see Figure 3) that 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. However, the reverse conclusion does 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.

Figure 6 shows, in the top panel, the neutrino fluence for GRB 2. The neutrino fluence follows the behavior of ${\tau }_{p\gamma }$ from Figure 5 (which is shown there for different examples). The average fluence per collision drops stronger than approximately $\propto {R}_{{\rm{C}}}^{-2}$. 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.

Zoom In Zoom Out Reset image size

5.2. Quasi-diffuse Neutrino Flux

We derive the all-sky quasi-diffuse ${\nu }_{\mu }+{\bar{\nu }}_{\mu }$ flux ${J}_{{\nu }_{\mu }}$ associated with a particular simulated GRB by scaling its fluence by the rate of long-duration GRBs, $\dot{N}=667$ per year, i.e., ${J}_{{\nu }_{\mu }}={{ \mathcal F }}_{{\nu }_{\mu }}\cdot \dot{N}\cdot {(4\pi )}^{-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 $\delta {t}_{\mathrm{eng}}$ (see Table 3) and ${N}_{\mathrm{sh}}$. We will discuss below how it depends on underlying assumptions, and what the corresponding fluxes are for GRBs 2–6.

Table 3.  Main Parameters of the Burst Simulation

Parameter	Description	Type	Units	Notes
Burst initialization
${N}_{\mathrm{sh}}$	Initial number of shells	Input	⋯	⋯
$\delta {t}_{\mathrm{eng}}$	Uptime of the engine	Input	s	⋯
${\rm{\Delta }}{t}_{\mathrm{eng}}$	Downtime of the engine	Input	s	⋯
${r}_{{N}_{\mathrm{sh}}}$	Distance from innermost shell to central emitter	Input	km	Default: 103 km
${r}_{\mathrm{dec}}$	Deceleration radius (circumburst medium starts)	Input	km	Default: $5.5\cdot {10}^{11}$ km
AΓ	Fluctuation amplitude of ${{\rm{\Gamma }}}_{k,0}$ distribution	Input	⋯	⋯
${E}_{\mathrm{kin},0}^{\mathrm{iso}}$	Initial bulk kinetic energy of shells	Input	erg	Common to all initial shells
z	Redshift of the emitter	Input	⋯	Used to calculate ${t}_{\mathrm{obs}}$, Equation (5): z = 2 by default
l	Initial shell width	Internal	km	Common to all initial shells: $l=c\cdot \delta {t}_{\mathrm{eng}}$
d	Initial separation between consecutive shells	Internal	km	Common to all initial shell pairs: d = l by default
${r}_{k,0}$	Initial radius of the kth shell,	Internal	km	${r}_{k,0}={r}_{{N}_{\mathrm{sh}}}+({N}_{\mathrm{sh}}-k)(l+d)$
${{\rm{\Gamma }}}_{k,0}$	Initial bulk Lorentz factor of the kth shell	Internal	⋯	Sampled from the pre-defined Γ distribution
${m}_{k,0}$	Initial mass of the kth shell	Internal	GeV	${m}_{k,0}={E}_{\mathrm{kin},0}^{\mathrm{iso}}/({{\rm{\Gamma }}}_{k,0}{c}^{2})$
Burst evolution
t	Time in the source frame	Internal	s	⋯
rk	Radius of the kth shell	Internal	km	Grows as ${r}_{k}={r}_{k,0}+c{\beta }_{k}t$
lk	Width of the kth shell	Internal	km	Changes only in collisions
${{\rm{\Gamma }}}_{k}$	Bulk Lorentz factor of the kth shell	Internal	⋯	Changes only in collisions
mk	Mass of the kth shell	Internal	GeV	Changes only in collisions
${\beta }_{k}$	Bulk speed of the kth shell	Internal	⋯	${\beta }_{k}=\sqrt{1-{{\rm{\Gamma }}}_{k}^{-2}}$
${V}_{\mathrm{iso},k}$	Isotropically equivalent volume of the kth shell	Internal	km3	${V}_{\mathrm{iso},k}=4\pi {r}_{k}^{2}{l}_{k}$
${E}_{\mathrm{kin},k}^{\mathrm{iso}}$	Bulk kinetic energy of the kth shell	Internal	erg	Changes only in collisions
${\rho }_{k}$	Mass density of the kth shell	Internal	GeV km−3	${\rho }_{k}={m}_{k}/{V}_{\mathrm{iso},k}$
Shell collisions
${m}_{{\rm{f}}({\rm{s}})}$	Mass of the fast (slow) colliding shell	Internal	GeV	⋯
${{\rm{\Gamma }}}_{{\rm{f}}({\rm{s}})}$	Bulk Lorentz factor of the fast (slow) shell	Internal	⋯	⋯
${{\rm{\Gamma }}}_{\mathrm{fs}(\mathrm{rs})}$	Lorentz factor of the forward (reverse) shock	Internal	⋯	See Equation (9)
${\beta }_{\mathrm{fs}(\mathrm{rs})}$	Speed of the forward (reverse) shock	Internal	⋯	${\beta }_{\mathrm{fs}(\mathrm{rs})}=\sqrt{1-{{\rm{\Gamma }}}_{\mathrm{fs}(\mathrm{rs})}^{-2}}$
${\beta }_{{\rm{m}}}$	Bulk speed of the merged shell	Internal	⋯	${\beta }_{{\rm{m}}}=\sqrt{1-{{\rm{\Gamma }}}_{{\rm{m}}}^{-2}}$
${\rho }_{{\rm{m}}}$	Mass density of the merged shell	Internal	GeV km−3	See Equation (10)
${t}_{\mathrm{coll}}$	Collision time (source frame)	Internal	s	⋯
${N}_{\mathrm{coll}}$	Total number of collisions	Output	⋯	⋯
${t}_{\mathrm{obs}}$	Collision time (observer's frame)	Output	s	See Equation (5)
${{\rm{\Gamma }}}_{{\rm{m}}}$	Bulk Lorentz factor of the merged shell	Output	⋯	See Equation (7)
${E}_{\mathrm{coll}}^{\mathrm{iso}}$	Internal energy liberated in the collision	Output	erg	See Equation (6)
${l}_{{\rm{m}}}$	Width of the merged shell	Output	km	See Equation (8)
RC	Collision radius	Output	km	⋯
$\delta {t}_{\mathrm{em}}$	Time for reverse shock to cross fast shell	Output	s	See Equation (17)

Note. All quantities are expressed in the source (engine) frame, except for ${t}_{\mathrm{obs}}$, which is in the observer's frame.

Download table as:  ASCIITypeset image

Figure 7 shows the fluxes for GRBs 1–6. For all but GRB 5, the flux is ${E}^{2}\,{J}_{{\nu }_{\mu }}\approx 2\cdot {10}^{-11}\,\mathrm{GeV}\ \,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}$ around 1 PeV—this is close to the value found for GRB 1 in Bustamante et al. (2015). For GRB 5, the flux is somewhat lower because it has fewer optically thick collisions, which is in agreement with Equation (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.

Zoom In Zoom Out Reset image size

IceCube has searched for correlations between neutrino arrival directions and positions of known GRBs (Abbasi et al. 2012; Aartsen et al. 2015a, 2016). No significant signal from GRBs has been found, in consistency with our prediction. Figure 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, Figure 7 shows that fluxes calculated with the multi-zone model are typically lower than with the one-zone model (Hümmer et al. 2012; see dashed curves in Figure 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 nonlinearly with the collision radius, the average value of the collision radius is, in general, not representative of the neutrino production. Nevertheless, we could define a single effective collision radius for neutrino production in the one-zone model that is different from the 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 parameters—such as the typical collision radius—are derived from gamma-ray observations (see Appendix A.4), which implies that parameters representative of gamma-rays may not be representative of the other messengers.

5.3. Cosmic Rays

5.4. Multi-messenger Emission from Different Radii

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 (Bustamante et al. 2015). 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 ($\gamma \gamma \to {e}^{+}{e}^{-}$) drives their energy down, so high-energy gamma-rays mainly come from large radii.

Zoom In Zoom Out Reset image size

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 ${\tau }_{\gamma \gamma }$ to pair production via $\gamma \gamma \to {e}^{+}{e}^{-}$. A photon of a certain energy escapes only if ${\tau }_{\gamma \gamma }\lt 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}_{\mathrm{obs}}$ is quite uniformly distributed in ${R}_{{\rm{C}}}$, especially between 10⁸ and ${10}^{10}\,\mathrm{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⁹ 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.

Zoom In Zoom Out Reset image size

For GRB 5 (Figure 9, bottom row), the situation is different. The average ${t}_{\mathrm{obs}}$ increases with ${R}_{{\rm{C}}}$ between 10⁹ and ${10}^{11}\,\mathrm{km}$, i.e., some early ($\lesssim 20$ s) gamma-ray detections come from mid-range radii, while all late detections tend to 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: Figure 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.¹² 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 higher-energy 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; see, e.g., Castignani et al. (2014).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since bursts with time delays between energy bands have narrow Γ distributions, they are associated with 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 as follows.

• 1.  
  In GRBs with light curves that have broader pulses overlaid with fast variability, time delays in different wavelength bands are possible.
• 2.  
  Compared to detection in the GBM energy band (10⁻⁶–10⁻² GeV), typical delays are a few seconds long in the LAT band (10⁻¹–10² GeV) and ∼10 s in the CTA band (10²–10⁶ GeV).
• 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}_{{\rm{C}}}$, where the radiation densities are higher and the gamma-rays cannot escape. This effect is only observable if observation time and ${R}_{{\rm{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. Figure 9 captures the central feature we observe in all examples that exhibit a lag: for these, there is a correlation between ${t}_{\mathrm{obs}}$ and ${R}_{{\rm{C}}}$ (within the range of ${R}_{{\rm{C}}}$ values in which the maximal photon energy depends on ${R}_{{\rm{C}}}$). This correlation can be traced back to the engine emitting shells in a relatively narrow range of values of Γ. The lags in Figure 10 are clearly visible, though 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.

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 the order of 10 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 gamma-ray 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 $\sim 2\cdot {10}^{-11}\,\mathrm{GeV}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{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 also be 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 high-energy neutrino counterpart to electromagnetic emission. The study presented here is a step toward a necessary, realistic multi-messenger understanding of GRBs. The observation of these neutrinos would be a smoking-gun signature of high baryonic loading, and thus of the paradigm that GRBs could be the sources of the UHECRs.

The authors thank John Beacom, Denise Boncioli, Valerie Connaughton, Željka Bošnjak, and Eli Waxman for useful discussion and comments on the manuscript. M.B. 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). M.B. acknowledges support from NSF Grant PHY-1404311. The work of K.M. is supported by NSF Grant PHY-1620777.

Here we describe our GRB jet simulation, based on Kobayashi et al. (1997) and Daigne & Mochkovitch (1998).

A.1. Overview

A.2. Burst Initialization

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.¹⁴ Therefore, the shell volume ${V}_{\mathrm{iso},k}=4\pi {r}_{k}^{2}{l}_{k}$ grows $\propto {r}_{k}^{2}$ and mass density ${\rho }_{k}={m}_{k}/{V}_{\mathrm{iso},k}$ decreases $\propto {r}_{k}^{-2}$. Since the shell mass and Lorentz factor are constant during propagation, its bulk kinetic energy ${E}_{\mathrm{kin},k}^{\mathrm{iso}}={{\rm{\Gamma }}}_{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.,

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{k,k+1}=\displaystyle \frac{{d}_{k,k+1}}{c({\beta }_{k+1}-{\beta }_{k})}\,,\end{eqnarray} \tag{ 3 }$$

where ${d}_{k,k+1}\equiv {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.,

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{\mathrm{next}}=\min [{\rm{\Delta }}{t}_{k,k+1}].\end{eqnarray} \tag{ 4 }$$

We increase the simulation time to $t\to t+{\rm{\Delta }}{t}_{\mathrm{next}}$. The collision radius ${R}_{{\rm{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

$$\begin{eqnarray}&&{t}_{\mathrm{obs}}=\left(\displaystyle \frac{D(z)-{R}_{{\rm{C}}}}{c}+t\right)(1+z),\end{eqnarray} \tag{ 5 }$$

where $D(z)$ is the light-travel distance to the emitter with redshift z. These equations satisfy well-known relations in the internal shock scenario, ${\rm{\Delta }}{t}_{\mathrm{next}}\approx 2{{\rm{\Gamma }}}^{2}(d/c)$ and ${R}_{{\rm{C}}}\approx 2{{\rm{\Gamma }}}^{2}d$. All shells propagate to their new positions ${r}_{k}\to {r}_{k}+c{\beta }_{k}{\rm{\Delta }}{t}_{\mathrm{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}_{\mathrm{obs}}=0$.)

We distinguish between the shell bulk kinetic energy, ${E}_{\mathrm{kin},k}^{\mathrm{iso}}$, and its internal energy, ${E}_{\mathrm{int},k}^{\mathrm{iso}}$. 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 of particles moving randomly inside the shell, measured in the shock rest frame.

In a collision, the kinetic energy of the two colliding shells is used partly as bulk kinetic energy for the new shell and partly as its internal energy. For simplicity, we assume that the new shell immediately cools by prompt particle emission; see Kobayashi & Sari (2001) for alternative treatments. While collision details depend on modeling of hydrodynamical properties (Daigne & Mochkovitch 2003), here we adopt the simple collision prescription for the relativistic limit introduced in Kobayashi et al. (1997), which we outline below. Figure 11 is a schematic illustration of the collision process.

In the collision of a slow (s) and a fast (f) shell, the internal energy of the merged (m) shell is the difference of kinetic energy before and after the collision, i.e.,

$$\begin{eqnarray}&&{E}_{\mathrm{coll}}^{\mathrm{iso}}=({{\rm{\Gamma }}}_{{\rm{f}}}{m}_{{\rm{f}}}+{{\rm{\Gamma }}}_{{\rm{s}}}{m}_{{\rm{s}}}){c}^{2}-{{\rm{\Gamma }}}_{{\rm{m}}}({m}_{{\rm{s}}}+{m}_{{\rm{f}}}){c}^{2}\,.\end{eqnarray} \tag{ 6 }$$

We assume that this amount of internal energy of the merged shell is radiated away as secondary particles. From momentum and energy conservation, and assuming ${{\rm{\Gamma }}}_{{\rm{f}}},{{\rm{\Gamma }}}_{{\rm{s}}}\gg 1$, the Lorentz factor of the merged shell is

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{m}}}\simeq \sqrt{\displaystyle \frac{{{\rm{\Gamma }}}_{{\rm{f}}}{m}_{{\rm{f}}}+{{\rm{\Gamma }}}_{{\rm{s}}}{m}_{{\rm{s}}}}{{m}_{{\rm{f}}}/{{\rm{\Gamma }}}_{{\rm{f}}}+{m}_{{\rm{s}}}/{{\rm{\Gamma }}}_{{\rm{s}}}}},\end{eqnarray} \tag{ 7 }$$

which reduces to ${{\rm{\Gamma }}}_{{\rm{m}}}\simeq \sqrt{{{\rm{\Gamma }}}_{{\rm{f}}}{{\rm{\Gamma }}}_{{\rm{s}}}}$ if ${m}_{{\rm{f}}}\simeq {m}_{{\rm{s}}}$. Its width is given by (Kobayashi et al. 1997; Aoi et al. 2010)

$$\begin{eqnarray}&&{l}_{{\rm{m}}}\simeq {l}_{{\rm{s}}}\displaystyle \frac{{\beta }_{\mathrm{fs}}-{\beta }_{{\rm{m}}}}{{\beta }_{\mathrm{fs}}-{\beta }_{{\rm{s}}}}+{l}_{{\rm{f}}}\displaystyle \frac{{\beta }_{{\rm{m}}}-{\beta }_{\mathrm{rs}}}{{\beta }_{{\rm{f}}}-{\beta }_{\mathrm{rs}}},\end{eqnarray} \tag{ 8 }$$

where ${\beta }_{\mathrm{fs}(\mathrm{rs})}=\sqrt{1-{{\rm{\Gamma }}}_{\mathrm{fs}(\mathrm{rs})}^{-2}}$ is the speed of the forward (reverse) shock, whose Lorentz factor is

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{fs}(\mathrm{rs})}={{\rm{\Gamma }}}_{{\rm{m}}}\sqrt{\displaystyle \frac{1+2{{\rm{\Gamma }}}_{{\rm{m}}}/{{\rm{\Gamma }}}_{{\rm{s}}({\rm{f}})}}{2+{{\rm{\Gamma }}}_{{\rm{m}}}/{{\rm{\Gamma }}}_{{\rm{s}}({\rm{f}})}}}\,.\end{eqnarray} \tag{ 9 }$$

The volume of the new shell is ${V}_{\mathrm{iso},{\rm{m}}}=4\pi {R}_{{\rm{C}}}^{2}{l}_{{\rm{m}}}$. The density is different between the shocked faster shell and shocked slower shell. The new shell has an average density obtained from the assumption of an inelastic collision ${m}_{m}\simeq {m}_{{\rm{f}}}+{m}_{{\rm{s}}}$ already implied in Equation (6), i.e.,

$$\begin{eqnarray}&&{\rho }_{{\rm{m}}}\simeq \displaystyle \frac{{l}_{{\rm{f}}}\cdot {\rho }_{{\rm{f}}}+{l}_{{\rm{s}}}\cdot {\rho }_{{\rm{s}}}}{{l}_{{\rm{m}}}}.\end{eqnarray} \tag{ 10 }$$

Therefore, its mass is ${m}_{{\rm{m}}}={V}_{\mathrm{iso},{\rm{m}}}{\rho }_{{\rm{m}}}$, and its kinetic energy is ${E}_{\mathrm{kin},{\rm{m}}}^{\mathrm{iso}}={{\rm{\Gamma }}}_{{\rm{m}}}{m}_{{\rm{m}}}{c}^{2}$. After the collision, the original fast shell is removed from the simulation and the new shell replaces the former slow shell. It is then propagated with the remaining shells in the jet.

If a shell reaches the circumburst medium, where it decelerates, it is removed from the simulation. Following Rees & Mészáros (1992), we assume ${r}_{\mathrm{dec}}=5.5\cdot {10}^{11}$ km for the radius at which this happens (see, e.g., Equation (15) in Mészáros 2006).

The simulation finishes when all shells have reached the circumburst medium, all shells have merged into a single remaining shell, or all remaining shells are ordered outward with increasing Lorentz factor, so that no more collisions are possible. The output lists ${N}_{\mathrm{coll}}$ collisions,

$$\begin{eqnarray}&&({t}_{\mathrm{obs},k},{R}_{{\rm{C}},k},{l}_{{\rm{m}},k},{{\rm{\Gamma }}}_{{\rm{m}},k},{E}_{\mathrm{coll},k}^{\mathrm{iso}}),\end{eqnarray} \tag{ 11 }$$

where $1\leqslant k\leqslant {N}_{\mathrm{coll}}$. The minimum ${t}_{\mathrm{obs},k}$ is taken to be the start of the observation time of the burst and is set to zero. Collisions are arranged so that ${t}_{\mathrm{obs},1}=0\leqslant {t}_{\mathrm{obs},2}\leqslant \ldots \leqslant \,{t}_{\mathrm{obs},{N}_{\mathrm{coll}}}$.

Figure 12 shows the time evolution (in the source frame) of macroscopic burst parameters in one of our simulations: average shell mass $\langle m\rangle /\langle {m}_{0}\rangle$ (subscripts of zero indicates values at simulation start), standard deviation of the Lorentz factor¹⁵ $\sigma ({\rm{\Gamma }})/\sigma ({{\rm{\Gamma }}}_{0})$, and total available internal energy of the burst ${E}_{\mathrm{int},\mathrm{tot}}^{\mathrm{iso}}/{E}_{\mathrm{int},\mathrm{tot},0}^{\mathrm{iso}}$. The latter is calculated directly as ${E}_{\mathrm{int},\mathrm{tot}}^{\mathrm{iso}}={{\rm{\Gamma }}}_{\mathrm{CM}}\sum {m}_{i}({\tilde{\beta }}_{i}^{2}/2)$, where $\tilde{\beta }$ is the speed in the CM-frame.¹⁶ 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}_{{\rm{\Gamma }}}\ll 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.

Zoom In Zoom Out Reset image size

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 ${\epsilon }_{e}$, ${\epsilon }_{B}$, and ${\epsilon }_{p}$, respectively. We assume energy equipartition between electrons and photons and fix ${\epsilon }_{e}={\epsilon }_{B}=1/12$ and ${\epsilon }_{p}=5/6$, since this yields the frequently used value of baryonic loading $1/{f}_{e}={\epsilon }_{p}/{\epsilon }_{e}=10$ (Waxman & Bahcall 1999; Abbasi et al. 2012). Thus, the kth collision dissipates an energy ${E}_{\gamma ,k}={\epsilon }_{e}{E}_{\mathrm{coll},k}$ as gamma-rays, and energy ${E}_{\gamma ,p}={\epsilon }_{p}{E}_{\mathrm{coll},k}$ as protons, and supports a magnetic energy density of ${U}_{B}={\epsilon }_{B}{E}_{\mathrm{coll},k}/{V}_{\mathrm{iso}}$. The latter translates into a magnetic field intensity, in the shock rest frame, of

$$\begin{eqnarray}{B}^{\prime } & \simeq & 44.7{\left(\displaystyle \frac{{{\rm{\Gamma }}}_{m}}{{10}^{2.5}}\right)}^{-1}{\left(\displaystyle \frac{{\epsilon }_{B}}{{\epsilon }_{e}}\right)}^{\tfrac{1}{2}}{\left(\displaystyle \frac{{E}_{\gamma ,k}}{{10}^{50}\mathrm{erg}}\right)}^{\tfrac{1}{2}}\\ & & \times {\left(\displaystyle \frac{{R}_{{\rm{C}}}}{{10}^{9}\mathrm{km}}\right)}^{-1}{\left(\displaystyle \frac{{l}_{{\rm{m}}}}{{10}^{3}\mathrm{km}}\right)}^{-\tfrac{1}{2}}\,\mathrm{kG}\,.\end{eqnarray} \tag{ 12 }$$

We normalize these individual energies by requiring that the total dissipated energy as gamma-rays, an experimentally accessible quantity,

$$\begin{eqnarray}&&{E}_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}\equiv \displaystyle \sum _{k=1}^{{N}_{\mathrm{coll}}}\,{E}_{\gamma ,k}^{\mathrm{iso}},\end{eqnarray} \tag{ 13 }$$

matches a given value ${E}_{\gamma ,\mathrm{norm}}^{\mathrm{iso}}={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 law, i.e.,

$$\begin{eqnarray}{n}_{\gamma }^{\prime }({\varepsilon }^{\prime })\equiv \displaystyle \frac{{{dn}}_{\gamma }^{\prime }}{d{\varepsilon }^{\prime }}={C}_{\gamma }^{\prime }\left\{\begin{array}{ll}{({\varepsilon }^{\prime }/{\varepsilon }_{\gamma ,\mathrm{break}}^{\prime })}^{-{\alpha }_{\gamma }}, & {\varepsilon }^{\prime }\lt {\varepsilon }_{\gamma ,\mathrm{break}}^{\prime }\\ {({\varepsilon }^{\prime }/{\varepsilon }_{\gamma ,\mathrm{break}}^{\prime })}^{-\beta \gamma }, & {\varepsilon }^{\prime }\geqslant {\varepsilon }_{\gamma ,\mathrm{break}}^{\prime }\end{array}\right..\end{eqnarray} \tag{ 14 }$$

Primed quantities are in the shock rest frame. We fix¹⁷ ${\alpha }_{\gamma }=1$, ${\beta }_{\gamma }=2$, and ${\varepsilon }_{\gamma ,\mathrm{break}}^{\prime }=1\,\mathrm{keV}$.

The photon spectrum in each shell is normalized via

$$\begin{eqnarray}&&{V}_{\mathrm{iso},k}^{\prime }{\int }_{{\varepsilon }_{\gamma ,\min }^{\prime }}^{{\varepsilon }_{\gamma ,\max }^{\prime }}d{\varepsilon }^{\prime }{\varepsilon }^{\prime }{n}_{\gamma }^{\prime }({\varepsilon }^{\prime })=\displaystyle \frac{{E}_{\gamma ,k}^{\mathrm{iso}}}{{{\rm{\Gamma }}}_{{\rm{m}},k}},\end{eqnarray} \tag{ 15 }$$

where the minimum and maximum energies are ${\varepsilon }_{\gamma ,\min }^{\prime }=0.2\,\mathrm{eV}$ and ${\varepsilon }_{\gamma ,\max }^{\prime }=1\,\mathrm{PeV}$, respectively (Aoi et al. 2010). Pair production via $\gamma \gamma \to {e}^{+}{e}^{-}$ may limit the maximum energy of escaping photons; see Figure 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 kth collision (in the observer's frame) as a peaked profile, with a fast rise and exponential decay ("FRED"), i.e.,

$$\begin{eqnarray}{L}_{\gamma ,k}({t}_{\mathrm{obs}})=\left\{\begin{array}{ll}0, & {t}_{\mathrm{obs}}\lt 0\\ {h}_{k}\left[1-\tfrac{1}{{(1+{ct}/{R}_{{\rm{C}},k})}^{2}}\right], & 0\leqslant {t}_{\mathrm{obs}}\lt {t}_{\mathrm{rise},k}\\ {h}_{k}\left\{\tfrac{1}{{[1+(t-\delta {t}_{\mathrm{em},k})c/{R}_{{\rm{C}},k}]}^{2}}\right. & \\ \left.\quad -\tfrac{1}{{(1+{ct}/{R}_{{\rm{C}},k})}^{2}}\right\}, & {t}_{\mathrm{obs}}\geqslant {t}_{\mathrm{rise},k}\end{array}\right.\end{eqnarray} \tag{ 16 }$$

where the emission timescale, i.e., the time at which the reverse shock crosses the fast shell, is

$$\begin{eqnarray}&&\delta {t}_{\mathrm{em},k}\equiv \displaystyle \frac{{l}_{{\rm{m}},k}}{c({\beta }_{{\rm{f}}}-{\beta }_{\mathrm{rs}})}\end{eqnarray} \tag{ 17 }$$

and the "rise time,"

$$\begin{eqnarray}&&{t}_{\mathrm{rise},k}\equiv \displaystyle \frac{\delta {t}_{\mathrm{em},k}}{2{{\rm{\Gamma }}}_{{\rm{m}},k}^{2}}(1+z),\end{eqnarray} \tag{ 18 }$$

is the time elapsed since the start of the emission until the peak luminosity is reached. For an illustration, see Figure 1 in Kobayashi et al. (1997). The peak luminosity is

$$\begin{eqnarray}&&{h}_{k}=\displaystyle \frac{{E}_{\gamma ,k}^{\mathrm{iso}}}{1+z}\displaystyle \frac{1}{{t}_{\mathrm{rise},k}}.\end{eqnarray} \tag{ 19 }$$

The time t in the source frame is related to ${t}_{\mathrm{obs}}$ through

$$\begin{eqnarray}&&t=\displaystyle \frac{2{{\rm{\Gamma }}}_{{\rm{m}},k}^{2}{t}_{\mathrm{obs}}}{1+z}\ .\end{eqnarray} \tag{ 20 }$$

Hence, the synthetic light curve L_γ is

$$\begin{eqnarray}&&{L}_{\gamma }({t}_{\mathrm{obs}})=\displaystyle \sum _{k=1}^{{N}_{\mathrm{coll}}}\,{L}_{\gamma ,k}({t}_{\mathrm{obs}}).\end{eqnarray} \tag{ 21 }$$

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 Equation (13). In our simulations, 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}_{\gamma ,\mathrm{norm}}^{\mathrm{iso}}$.

In Section 6, we showed light curves in different energy bands. Each band spans the range $[{E}_{\gamma ,\mathrm{low}}^{\mathrm{band}},{E}_{\gamma ,\mathrm{high}}^{\mathrm{band}}]$, with ${E}_{\gamma ,\mathrm{low}}^{\mathrm{band}}={10}^{-6}$, 10⁻¹, 10² GeV, and ${E}_{\gamma ,\mathrm{high}}^{\mathrm{band}}={10}^{-2}$, 10², 10⁶ GeV, for Fermi-GBM, Fermi-LAT, and CTA, respectively. To compute the gamma-ray contribution of the kth collision to each band, we calculate the fraction

$$\begin{eqnarray}&&{f}_{k}^{\mathrm{band}}\\ \quad =\left\{\begin{array}{ll}0, & {E}_{\gamma ,\max ,k}\lt {E}_{\gamma ,\mathrm{low}}^{\mathrm{band}}\\ 1, & {E}_{\gamma ,\max ,k}\gt {E}_{\gamma ,\mathrm{high}}^{\mathrm{band}}\\ \tfrac{\mathrm{log}({E}_{\gamma ,\max ,k})-\mathrm{log}({E}_{\gamma ,\mathrm{low}}^{\mathrm{band}})}{\mathrm{log}({E}_{\gamma ,\mathrm{high}}^{\mathrm{band}})-\mathrm{log}({E}_{\gamma ,\mathrm{low}}^{\mathrm{band}})} & ,\ \ {E}_{\gamma ,\mathrm{low}}^{\mathrm{band}}\leqslant {E}_{\gamma ,\max ,k}\leqslant {E}_{\gamma ,\mathrm{high}}^{\mathrm{band}}\end{array}\right.,\end{eqnarray} \tag{ 22 }$$

where ${E}_{\gamma ,\max ,k}$ is the maximum gamma-ray energy emitted by the collision. The light curve for each band is then simply calculated using Equation (21), replacing ${L}_{\gamma ,k}\to {f}_{k}^{\mathrm{band}}{L}_{\gamma ,k}$. Note that this simplified treatment assumes that the product of instrument response times flux is approximately constant within the anticipated energy band—which is typically a good estimate within the energy bands the instruments are designed for. A more detailed model for the instrument response or a different shape of the target photon spectrum at TeV energies—which we assumed to be a power law beyond the break—will not affect the qualitative shape of the light curves, but may slightly change the relative power in different energy bands or light-curve peaks.

The burst duration and variability timescale ${t}_{{\rm{v}}}$ are derived from the light curve. For the duration, we use T₉₀, the time elapsed between the detection of 5% and 95% of the total gamma-ray energy, i.e.,

$$\begin{eqnarray}&&{T}_{90}\equiv {t}_{95}-{t}_{5},\end{eqnarray} \tag{ 23 }$$

where

$$\begin{eqnarray}&&{\int }_{0}^{{t}_{f}}{L}_{\gamma }({t}_{\mathrm{obs}}){{dt}}_{\mathrm{obs}}=f\,\displaystyle \frac{{E}_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}}{1+z}\end{eqnarray} \tag{ 24 }$$

and f = 0.05 or 0.95. The variability timescale is estimated as

$$\begin{eqnarray}&&{t}_{{\rm{v}}}={T}_{90}/{N}_{\mathrm{coll}}.\end{eqnarray} \tag{ 25 }$$

This procedure yields values of ${t}_{{\rm{v}}}$ close to $\delta {t}_{\mathrm{eng}}$.

In some cases (e.g., Figure 7), we have compared simulation results to "standard" estimators from the one-zone model:

$$\begin{eqnarray}&&{t}_{{\rm{v}}}^{\mathrm{std}}=\displaystyle \frac{d+l}{c}(1+z),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{T}^{\mathrm{std}}={N}_{\mathrm{coll}}\,{t}_{{\rm{v}}}^{\mathrm{std}},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{m}}}^{\mathrm{std}}=\langle {{\rm{\Gamma }}}_{{\rm{m}}}\rangle ,\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{R}_{{\rm{C}}}^{\mathrm{std}}=2{({{\rm{\Gamma }}}_{{\rm{m}}}^{\mathrm{std}})}^{2}\displaystyle \frac{{{ct}}_{{\rm{v}}}^{\mathrm{std}}}{1+z}.\end{eqnarray} \tag{ 29 }$$

A.5. Neutrinos and Cosmic-rays

In analogy to gamma-rays, the proton and neutrino spectra of the complete burst are obtained by summing over the spectra emitted by all the individual collisions.

We compute secondary particle production using the NeuCosmA software (Baerwald et al. 2012; Hümmer et al. 2012; Baerwald et al. 2013). This assumes a proton density ${n}_{p}^{\prime }\propto {E}_{p}^{\prime -2}\exp (-{E}_{p}^{\prime }/{E}_{p,\max }^{\prime })$, with the maximum proton energy ${E}_{p,\max }^{\prime }$ obtained from balancing the acceleration rate (with perfect efficiency) with synchrotron, adiabatic, and photohadronic energy loss rates (Baerwald et al. 2013). The proton density is normalized like the photon density, Equation (15), but replacing ${E}_{\gamma ,k}^{\mathrm{iso}}\to {E}_{\gamma ,k}^{\mathrm{iso}}/{f}_{e}$. Secondary pion, muon, and kaon spectra, and, consequently, neutrino spectra are computed as in Hümmer et al. (2012), which includes magnetic field effects on the secondaries, state-of-the-art normalization of the spectra, helicity-dependent muon decays, and flavor mixing.

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

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{n}_{e,k}^{\prime }\simeq \displaystyle \frac{{m}_{k}}{{m}_{p}\,{V}_{\mathrm{iso},k}^{\prime }},\end{eqnarray} \tag{ 30 }$$

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

$$\begin{eqnarray}&&{\tau }_{\mathrm{Th},k}^{\prime }\simeq \displaystyle \frac{1}{{n}_{e,k}^{\prime }\,{\sigma }_{\mathrm{Th}}\,{l}_{k}^{\prime }},\end{eqnarray} \tag{ 31 }$$

with ${\sigma }_{\mathrm{Th}}\approx 66.52$ fm² the Thomson cross-section. A collision is sub-photospheric if ${\tau }_{\mathrm{Th}}^{\prime }\gt 1$.

For a burst with ${E}_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}\simeq {10}^{53}\,\mathrm{erg}$, ${\epsilon }_{e}=1/12$, and a dissipation efficiency of $\varepsilon =25 \%$ (such as GRB 1), the initial kinetic energy per shell is about ${10}^{51.6}\,\mathrm{erg}$ if 1000 collisions occur. This yields ${m}_{k,0}={E}_{\mathrm{kin},0}^{\mathrm{iso}}/{{\rm{\Gamma }}}_{k,0}{c}^{2}\simeq {10}^{49}\,\mathrm{erg}$ for ${{\rm{\Gamma }}}_{k,0}\sim 500$. From Equations (30) and (31), the photospheric radius of the kth shell is

$$\begin{eqnarray}&&{R}_{\mathrm{ph},k}\simeq 1.8\cdot {10}^{8}\,\mathrm{km}{\left(\displaystyle \frac{{m}_{k}}{{10}^{49}\mathrm{erg}}\right)}^{1/2}.\end{eqnarray} \tag{ 32 }$$

Bustamante et al. (2015) showed that, since the photohadronic optical depth scales similarly to Equation (31) (replacing ${\sigma }_{\mathrm{Th}}\to {\sigma }_{p\gamma }$), the pion production efficiency at the photosphere is independent of the isotropic volume and only weakly dependent on ${{\rm{\Gamma }}}_{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.

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.

Zoom In Zoom Out Reset image size

Table 4.  Description of Simulated GRBs 7–10

Model	${{\rm{\Gamma }}}_{\mathrm{0,1}}$	${A}_{{\rm{\Gamma }},1}$	${{\rm{\Gamma }}}_{\mathrm{0,2}}$	${A}_{{\rm{\Gamma }},2}$	Tp	${N}_{\mathrm{up}}$	${N}_{\mathrm{down}}$	${E}_{\gamma ,\mathrm{norm}}^{\mathrm{iso}}$ (erg)	Description
7	5	⋯	1500	⋯	⋯	⋯	⋯	1053	Box-like Γ distribution
8	50	0.1	500	1.0	⋯	⋯	⋯	1053	Linear speedup of Γ
9	500	1.0	⋯	⋯	⋯	100	100	1052	100 pulses separated by 100 pulse-times
10	500	1.0	⋯	⋯	⋯	300	350	1052	300 pulses separated by 350 pulse-times

Note. Common values for all models: ${N}_{\mathrm{sh}}=1000$, $\delta {t}_{\mathrm{eng}}={10}^{-2}$ s, $d=l=c\cdot \delta {t}_{\mathrm{eng}}$, ${r}_{\min }={10}^{3}$ km, ${r}_{\mathrm{dec}}=5.5\cdot {10}^{11}$ km, z = 2, ${\epsilon }_{e}={\epsilon }_{B}=1/12$, ${\epsilon }_{p}=5/6$, $\eta =1.0$ (acceleration efficiency, see Baerwald et al. 2013). See Table 3 for an explanation of each parameter.

Download table as:  ASCIITypeset image

The light curves of GRB 7 and GRB 8 are similar to that of GRB 1; see Figure 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.

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}_{\gamma ,\mathrm{norm}}^{\mathrm{iso}}={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.

Zoom In Zoom Out Reset image size

Another modification is to assume that only a fraction η of the internal energy in each collision, Equation (6), is attributed to the non-thermal spectra of the secondaries, while $1-\eta$ is instantly reconverted into kinetic energy of the merged shell. This fraction $1-\eta$ 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: ${{\rm{\Gamma }}}_{{\rm{m}}}=[(1-\eta ){E}_{\mathrm{coll}}^{\mathrm{iso}}+{E}_{\mathrm{kin},{\rm{m}}}^{\mathrm{iso}}]/{m}_{{\rm{m}}}$.

However, 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 be 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 that the amount of internal energy used for computing the secondary production is still given by Equation (6), 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 $\eta =\xi =0.5$. In both cases, the burst was normalized to yield 10⁵³ erg in gamma-rays when adding sub-photospheric and super-photospheric 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 with these two modifications is slightly lower than the one associated with the reference case, especially if only a fraction η of the energy is radiated. Therefore, the minimal super-photospheric flux prediction of $\sim {10}^{-11}$ GeV cm⁻² s⁻¹ s⁻¹ holds.

The neutrino flux scales with the fraction η going into the non-thermal spectra (and magnetic field), which means that, for $\eta =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-\xi$. This means that for the combined case—a fraction η into non-thermal spectra and a fraction $1-\xi$ 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.