Impact of the collision model on the multi-messenger emission from Gamma-Ray Burst internal shocks

TL;DR

This paper investigates how different collision models in Gamma-Ray Burst internal shocks affect multi-messenger emissions, finding that the Ultra Efficient Shock scenario enhances neutrino production and energy conversion efficiency, with implications for IceCube-Gen2 detection.

Contribution

It introduces a detailed comparison of collision models in GRB internal shocks and their impact on multi-messenger emission predictions, especially highlighting the Ultra Efficient Shock scenario.

Findings

Neutrino flux levels are estimated at 10^{-11} to 10^{-10} GeV cm^{-2} s^{-1} sr^{-1}.

Ultra Efficient Shock scenario shows higher energy conversion efficiency.

Shell separation assumptions are rarely justified in hydrodynamical multi-collision simulations.

Abstract

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

Tables

Table 1. Table 1: Qualitative comparison of the different collision models. Here N inital shells subscript superscript 𝑁 shells inital N^{\mathrm{shells}}_{\mathrm{inital}} is the number of shells ejected by the source, l initial shells subscript superscript 𝑙 shells initial l^{\mathrm{shells}}_{\mathrm{initial}} are the initial shell widths (assumed to be equal for all shells) and N post − coll . shells subscript superscript 𝑁 shells post coll N^{\mathrm{shells}}_{\mathrm{post-coll.}} is the number of shells produced by a single collision. In each collision, the dissipated energy is given by η ≡ E diss / E int 𝜂 subscript 𝐸 diss subscript 𝐸 int \eta\equiv E_{\mathrm{diss}}/E_{\mathrm{int}} , and the width of the shell(s) after the collision is given by l post − coll . shells subscript superscript 𝑙 shells post coll l^{\mathrm{shells}}_{\mathrm{post-coll.}} . In the case of the PLUTO model, the number of post-collision shells as well as their widths are obtained from an analysis of the post-collision mass density profile from the simulation.

\toprule		Initial setup		Single collision result
Model name	$N_{\mathrm{inital}}^{\mathrm{shells}}$	$\frac{l_{\mathrm{initial}}^{\mathrm{shells}}}{c}$	$N_{\mathrm{post}-\mathrm{coll}.}^{\mathrm{shells}}$	$\eta$	$l_{\mathrm{post}-\mathrm{coll}.}^{\mathrm{shells}}$
Reference	1000	0.01s	1	1
Reduced Efficiency	1000	0.01s	1	0.5	$=l_r+l_s$
Ultra Efficient	125	0.08s	2	0.5	$l_r$, $l_s$
PLUTO	1000	0.01s	1-2	0.5

Table 2. Table 2: Results for the different models: ϵ ≡ E diss tot / E kin tot italic-ϵ subscript superscript 𝐸 tot diss subscript superscript 𝐸 tot kin \epsilon\equiv E^{\mathrm{tot}}_{\mathrm{diss}}/E^{\mathrm{tot}}_{\mathrm{kin}} describes the overall dissipation efficiency of the burst, t ν subscript 𝑡 𝜈 t_{\nu} the (observable) time variability of the light curve assuming a redshift of z = 2 𝑧 2 z=2 , N coll subscript 𝑁 coll N_{\mathrm{coll}} the total number of collisions and E kin tot subscript superscript 𝐸 tot kin E^{\mathrm{tot}}_{\mathrm{kin}} the total initial kinetic energy. The energies E ν , tot iso superscript subscript 𝐸 𝜈 tot iso E_{\nu,\mathrm{tot}}^{\mathrm{iso}} , E γ , tot iso superscript subscript 𝐸 𝛾 tot iso E_{\gamma,\mathrm{tot}}^{\mathrm{iso}} and E p , tot iso superscript subscript 𝐸 𝑝 tot iso E_{p,\mathrm{tot}}^{\mathrm{iso}} are the total energies released in neutrinos, gamma-rays and UHECRs above 10 10 ​ GeV superscript 10 10 GeV 10^{10}\,\mathrm{GeV} in the engine frame (where only super-photospheric collisions are taken into account), and E p , max subscript 𝐸 𝑝 max E_{p,\mathrm{max}} is the maximal cosmic ray energy achieved in the engine frame. We show the mean value over all statistical realizations as well as the standard deviation obtained from the Monte Carlo simulation.

\toprule		$ϵ$ [%]	$t_{\nu }$ [ms]	$N_{\mathrm{coll}}$	$E_{\mathrm{kin}}^{\mathrm{tot}}$ [${10}^{54}$ erg]	$E_{p,\mathrm{tot}}^{\mathrm{iso}}/E_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}$	$E_{\nu ,\mathrm{tot}}^{\mathrm{iso}}/E_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}$	$E_{p,\max }$ [${10}^{12}$ GeV]
Reference	35.8 $\pm$ 1.4	55.2 $\pm$ 1.3	970.1 $\pm$ 3.3	1.75 $\pm$ 0.07	0.42 $\pm$ 0.03	0.29 $\pm$ 0.05	1.2 $\pm$ 0.4
Reduced Efficiency	17.9 $\pm$ 0.7	54.8 $\pm$ 1.3	976.0 $\pm$ 3.3	3.50 $\pm$ 0.13	0.56 $\pm$ 0.04	0.24 $\pm$ 0.05	1.2 $\pm$ 0.4
Ultra Efficient	36.0 $\pm$ 4.3	47.5 $\pm$ 10.7	1107 $\pm$ 220	1.76 $\pm$ 0.22	0.62 $\pm$ 0.06	0.14 $\pm$ 0.06	1.2 $\pm$ 0.5
PLUTO	21.2 $\pm$ 1.4	50.2 $\pm$ 1.5	1055 $\pm$ 9	2.95 $\pm$ 0.20	0.96 $\pm$ 0.4	0.18 $\pm$ 0.04	1.4 $\pm$ 0.6

Table 3. Table 3: Same as Tab.  6a , but for an engine ejecting shells with a constant luminosity ( E ˙ = c ​ o ​ n ​ s ​ t . ˙ 𝐸 𝑐 𝑜 𝑛 𝑠 𝑡 \dot{E}=const. ) instead of a constant mass outflow (labeled Reference E) and a modified version of PLUTO model in which the PLUTO simulations don’t account for energy dissipation. In both cases, all of the internal energy is radiated.

	$ϵ$ [%]	$t_{\nu }$ [ms]	$N_{\mathrm{coll}}$	$E_{\mathrm{eng}}$ [${10}^{54}$ erg]	$E_{p,\mathrm{tot}}^{\mathrm{iso}}/E_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}$	$E_{\nu ,\mathrm{tot}}^{\mathrm{iso}}/E_{\gamma ,\mathrm{tot}}^{\mathrm{iso}}$	$E_{p,\max }$ [${10}^{12}$ GeV]
Reference	35.8 $\pm$ 1.4	55.2 $\pm$ 1.3	970.1 $\pm$ 3.3	1.75 $\pm$ 0.07	0.42 $\pm$ 0.03	0.29 $\pm$ 0.05	1.2 $\pm$ 0.4
Reference E	25.8 $\pm$ 1.0	54.5 $\pm$ 0.5	987.5 $\pm$ 2.8	2.4 $\pm$ 0.1	0.27 $\pm$ 0.01	0.080 $\pm$ 0.006	0.26 $\pm$ 0.02
PLUTO no diss.	45.3 $\pm$ 2.9	41.2 $\pm$ 1.1	1270 $\pm$ 19	1.38 $\pm$ 0.09	0.54 $\pm$ 0.16	0.26 $\pm$ 0.04	1.29 $\pm$ 0.30

Equations

$$\ln \left(\frac{{\Gamma }_k-1}{{\Gamma }_0-1}\right)=A\cdot x,$$

$${\Gamma }_m=\sqrt{\frac{m_r{\Gamma }_r+m_s{\Gamma }_s}{m_r/{\Gamma }_r+m_s/{\Gamma }_s}}$$

$$E_{\text{int},m}=m_r{\Gamma }_r+m_s{\Gamma }_s-(m_r+m_s){\Gamma }_m.$$

$$\eta =\frac{E_{\text{diss},m}}{E_{\text{int},m}}$$

$${\Gamma }_{\mathrm{fs}(\mathrm{rs})}={\Gamma }_m\cdot \sqrt{\frac{1+2{\Gamma }_m/{\Gamma }_{s(r)}}{2+{\Gamma }_m/{\Gamma }_{s(r)}}}.$$

$$\delta t_{\mathrm{em}}=\frac{l_r}{{\beta }_r-{\beta }_{\mathrm{rs}}}.$$

$$l_m=l_s\frac{{\beta }_{\mathrm{fs}}-{\beta }_m}{{\beta }_{\mathrm{fs}}-{\beta }_s}+l_r\frac{{\beta }_m-{\beta }_{\mathrm{rs}}}{{\beta }_r-{\beta }_{\mathrm{rs}}}.$$

$$ϵ\equiv E_{\mathrm{diss}}^{\mathrm{tot}}/E_{\mathrm{kin}}^{\mathrm{tot}}$$

$$d=\frac{d^{\ast }}{{\rho }_{\mathrm{Max}}}=\frac{{\rho }_{\mathrm{Max}}-{\rho }_{\mathrm{Min}}}{{\rho }_{\mathrm{Max}}},$$

$$\begin{array}{rl}{\bar{\Gamma }}_{\mathrm{r}}& =(M^2+(m_{\mathrm{r}}{)}^2-(m_{\mathrm{s}}{)}^2)/2m_{\mathrm{r}}M,\\ {\bar{\Gamma }}_{\mathrm{s}}& =(M^2-(m_{\mathrm{s}}{)}^2+(m_{\mathrm{r}}{)}^2)/2m_{\mathrm{s}}M\end{array}$$

$$\begin{array}{rl}{\Gamma }_{\mathrm{r}}& ={\bar{\Gamma }}_{\mathrm{r}}{\gamma }_{\mathrm{m}}-\sqrt{({\bar{\Gamma }}_{\mathrm{r}}^2-1)({\gamma }_{\mathrm{m}}^2-1)}\text{and}\\ {\Gamma }_{\mathrm{s}}& ={\bar{\Gamma }}_{\mathrm{s}}{\gamma }_{\mathrm{m}}-\sqrt{({\bar{\Gamma }}_{\mathrm{s}}^2-1)({\gamma }_{\mathrm{m}}^2-1)}.\end{array}$$

$$\frac{\partial \rho e}{\partial t}=\{\begin{array}{ll}-\chi \cdot \rho e& E_{\mathrm{diss},\mathrm{total}}\le {\eta }^{\ast }{E}_{\mathrm{int},\mathrm{th}}\\ 0& \mathrm{else}\end{array}$$

Full text

††software: PLUTO v4.2 (Mignone et al., 2007), NeuCosmA (Biehl et al., 2018)

1 Introduction

Gamma-Ray Bursts (GRBs) have been proposed as plausible candidates for the origin of the Ultra-High Energy Cosmic Rays (UHECRs) and neutrinos, invoking photohadronic interactions in the fireball scenario (Waxman & Bahcall, 1997). It is however evident from GRB stacking searches that GRBs cannot be the dominant source of the observed diffuse astrophysical neutrino flux (Abbasi et al., 2012; Aartsen et al., 2017). In radiation models of high-luminosity GRBs where all emission regions look alike (one-zone models), the nominal predictions for the neutrino flux are in tension with these stacking searches in spite of recently improved neutrino flux estimates (Hummer et al., 2012; Li, 2012; He et al., 2012). Dedicated scans of the source parameters (including the baryonic loading, the luminosity in protons versus gamma-rays) and fits to the UHECR spectrum and composition data confirm that the simple one zone emission picture is in tension with neutrino data for most of the parameter space (Baerwald et al., 2015; Biehl et al., 2018).

As a possible solution multi-collision models have been proposed, which tend to predict a lower neutrino flux (Globus et al., 2015; Bustamante et al., 2014, 2017) for the same baryonic loading. These models also demonstrate that the different messengers (neutrinos, cosmic rays, gamma-rays) are predominantly emitted at different radii within the GRB jet (Bustamante et al., 2014). Also the stochasticity pattern and the time delays between different energy bands in the light curves can contain information on the neutrino production efficiency (Bustamante et al., 2017).

As specific implementations of the fireball model (Rees & Meszaros, 1992, 1994), most GRB multi-collision models invoke emission from internal shocks  (Kobayashi et al., 1997; Daigne & Mochkovitch, 1998), in which a set of plasma shells is emitted from an intermittent central engine with a specific distribution of Lorentz factors. The emitted shells can have, for instance, equal masses or equal energies, which are related through the relation $E_{\mathrm{kin}}=\Gamma\,m$. Because of the different velocities, the shells will eventually catch up with each other and collide inelastically, forming a merged shell that continues to propagate through the jet. In each shell collision, shocks form that accelerate particles, converting internal energy (from the inelastic collision) into non-thermal radiation.

In principle, multi-collision models provide a natural explanation for the the fast time variability and variety of light curve shapes observed in GRBs; see e.g. Kobayashi et al. (1997); Daigne & Mochkovitch (1998); Spada et al. (2000); Beloborodov (2000); Kobayashi & Sari (2001); Daigne & Mochkovitch (2003). A fundamental problem is, however, the moderate dissipation efficiency of kinetic energy into non-thermal particles. On one hand, the afterglow emission caused by the shells running into the circumburst medium after the prompt phase might be too high if a large amount of energy remains in the jet. One the other hand extremely powerful engines are required to reach the observed gamma-ray luminosities. The Ultra Efficient Shock scenario has been proposed as a potential solution to this problem (Kobayashi & Sari, 2001), in which a fraction of the internal energy is re-converted into kinetic energy. Consequently, the two shells survive and remain in the system after the collision, resulting in swift thermalization. In this scenario a high overall dissipation efficiency can be achieved, even if the efficiency of each individual two-shell collision is low. A key ingredient of the scenario is the assumption that two shells remain after each collision, whereas hydrodynamical studies and analytical estimates (Kino et al., 2004) demonstrate that this assumption is depending on the collision parameters and that most collisions likely result in only one merged shell.

In this work, we revisit the multi-collision multi-messenger models in Bustamante et al. (2014, 2017) from the point of view of the collision dynamics and study the impact on the production of multiple messengers. We employ the methods from Baerwald et al. (2012); Hummer et al. (2012); Biehl et al. (2018) for the radiation calculations and use broken power-law target photon spectra that resemble observations. After recapping the previous approach, i.e., the fully inelastic case with equal mass injection in Section 2, we scan the parameter space for the two-shell collision for configurations with two post-collision shells, and scrutinize the plausability of this assumption in the fireball evolution. The multi-messenger production is studied in Section 4 for alternative collision models: A) a version of the reference model, in which a fraction of the internal energy is converted into adiabatic expansion of the merged shell; B) the ideal Ultra Efficient scenario according to Kobayashi & Sari (2001); and C) a “hybrid” model, where PLUTO is used to determine the fate of each two-shell collision individually. We discuss the impact on the observables (light curves and neutrino fluxes) in Section 5 and conclude in Section 7.

2 Methods and reference model

As the reference case, we pick the collision model from Bustamante et al. (2017), based on Kobayashi et al. (1997) and recapitulate the main ideas and mathematical expressions before entering the discussion of alternatives in the next sections.

The relativistic outflow of the jet is discretized as a one-dimensional sequence of spherical plasma shells with different velocities and masses. When a faster shell catches up with a slower one, they interact and collisionless shocks convert kinetic energy into internal energy (which will be radiated as non-thermal particles). During the interaction, the two shells merge into a single new shell which will continue to propagate as a part of the jet (see Bustamante et al. (2017) and first subsection for a detailed description). The internal energy remaining from the collision is converted into radiation through the interaction of accelerated particles, based on the model from an updated version of the NeuCosmA software (Biehl et al., 2018).111The update impacts the description of the optically thick (to photo-hadronic interactions) case, see App. C of Biehl et al. (2018) for details. This leads to slightly lower predicted cosmic ray fluxes at the highest energies. The baryonic loading is defined with respect to the injection luminosity (previously, the steady state density), resulting in a more transparent evaluation of the energy budget.

In this section, the “GRB 1” in Bustamante et al. (2017) acts as a reference parameters set, but with an equal-mass instead of equal-energy setup, resulting in a few quantitative differences, which are discussed in Appendix C.1. It will become clear in the next sections why this choice allows for a better comparison with the alternative models. For the sake of simplicity, we focus on the proton-only case in this study. Changing the composition would not affect the level of the neutrino flux, as it mainly depends on the energy injected per nucleon. However the maximal energy of neutrino flux is lower for heavier injection, as the maximal rigidity is reduced by $Z/A$. See Biehl et al. (2018) for a detailed discussion of nuclear injection in GRBs.

Collision model

A plasma shell is characterized by its mass ($m$), width ($l$) and the Lorentz factor ($\Gamma$) with respect to the engine frame. When a fast (“rapid”, index $r$) shell collides with a slow (index $s$) shell, forward and reverse shocks develop. They propagate from the contact discontinuity (”CD”, which is the surface separating the initial slow from the fast shell) through the joint density profile initially compressing the matter (see Fig. 2 for an illustration of the setup). After the shock waves reach the edges of the density profile, two rarefaction waves develop that propagate back towards the CD. They decompress the matter and reconvert internal energy into kinetic energy. The kinetic properties of the $k$-th shell are determined by the Lorentz factor $\Gamma_{k}$, the mass $m_{k}$ or the kinetic energy $E_{\text{kin},k}$, the radius $R_{k}$ and the width $l_{k}$. The other properties, such as the speed $\beta_{k}$, the volume $V_{k}=4\,\pi\,R_{k}^{2}\,l_{k}$ and the density $\rho_{k}$, can be derived. In the following, all unprimed quantities (unless explicitly given in the observer’s frame) are given in the engine frame, and all primed quantities refer to the shock rest (or merged shell) frame.

The simulation of the reference model initially contains 1000 shells with their Lorentz factors $\Gamma_{k}$’s sampled from a log-normal distribution

[TABLE]

where $x$ is distributed as a Gaussian $P(x)dx=\exp(-x^{2})/\sqrt{2\pi}dx$ with $\Gamma_{0}=500$ and $A=1.0$. We further assume that the spatial separation for all shells is equal to their width ($l_{k}=L_{k}=0.01\ \mathrm{c\ s}$, equivalent to the temporal separation $\delta t_{k}=0.01\ \mathrm{s}$ for relativistic shells) and that the innermost shell starts at a radius $R_{\text{min}}$. The initial distribution of kinetic energies $E_{\text{kin},k}$ is a free model choice, and typical assumptions are equal shell masses $m_{k}$, energies $E_{\text{kin},k}$ or densities $\rho_{k}$. As previously mentioned, the reference model uses the equal-mass case. From this initial setup the system evolves until all shells have reached the circumburst medium at $R_{\text{max}}$.

For the collision between a rapid and a slow shell the properties for the merged shell (index $m$) are determined from simple analytical expressions following Kobayashi et al. (1997). The Lorentz factor $\Gamma_{m}$ follows from the conservation of momentum as

[TABLE]

(using the approximation $\beta=1-1/(2\Gamma^{2})$ that holds for $\Gamma\gg 1$). The internal energy available for radiation is then the difference of kinetic energies before and after the (inelastic) collision:

[TABLE]

We define the efficiency of each collision as the ratio between the dissipated energy and the produced internal energy:

[TABLE]

The idealized assumption of $\eta=1$ in the Reference model will be modified for the alternative models in Section (4).

The Lorentz factors of the forward (fs) and the reverse shock (rs) are determined from

[TABLE]

The timescale of emission $t_{\mathrm{em}}$ is estimated from the time taken by the reverse shock to cross the rapid shell

[TABLE]

This timescale is observed Doppler-boosted in the observer’s frame. The width of the compressed merged shell is computed from

[TABLE]

For Eq. (5) to Eq. (7) we follow the description given in Kobayashi et al. (1997). This formulation might not be valid in the non-relativistic regime, however collisions in the non-relativistic regime dissipate little energy and therefore contribute only marginally to the multi-messenger emission. Even though $t_{\mathrm{em}}$ and $l_{m}$ might be slightly miss-estimated for those collisions, the impact on our main results is therefore small.

Radiation model

The collision parameters derived in the previous section are now used for the computation of the energy density and the emission spectra in the merged shell (primed frame). It is assumed that fractions of the injected luminosity $L^{\prime}_{\text{diss},m}\simeq E^{\prime}_{\text{diss}}\,c/l^{\prime}_{m}$ are distributed into protons ($\epsilon_{p}$), electrons ($\epsilon_{e}$), and the magnetic field ($\epsilon_{B}$), such that $\epsilon_{p}+\epsilon_{e}+\epsilon_{B}=1$. Normally it is assumed that the non-thermal electrons loose energy quickly via synchrotron radiation, implying that gamma-rays will carry a comparable amount of energy $\epsilon_{\gamma}\simeq\epsilon_{e}$. We assume a baryonic loading of $\epsilon_{p}/\epsilon_{\gamma}=10$. It is convenient to relate these quantities to the (observed gamma rays), consequently $\epsilon_{\gamma}=\epsilon_{B}=1/12$ and $\epsilon_{p}=10/12$.

We assume a power-law spectrum motivated by Fermi shock acceleration ($\propto(E^{\prime}_{p})^{-2}\exp(-E^{\prime}_{p}/E^{\prime}_{p,\text{max}})$) for the injected proton component. The maximal energy $E^{\prime}_{p,\text{max}}$ is found by comparing the timescales of efficient acceleration $t^{\prime-1}\simeq c/R^{\prime}_{L}$ ($R^{\prime}_{L}$ is the Larmor radius), photohadronic and adiabatic energy losses (for a detailed discussion of the maximum energies see also Samuelsson et al. (2018)). The proton spectrum is normalized to the available luminosity fraction ($\epsilon_{p}$) as outlined above.

For the target photon spectrum, we do not perform explicit self-consistent radiation calculations, but instead assume a shape motivated by observations– a broken power law with spectral indices $-1$ and $-2$ below and above the break energy $\varepsilon^{\prime}_{\text{break}}\simeq 1\ \mathrm{keV}$, respectively Gruber et al. (2014). The maximal photon energy is assumed to be limited to $\varepsilon^{\prime}_{\text{max}}=10^{6}\ \mathrm{GeV}$ and reduced in case $\gamma\gamma$-annihilation sets in at lower energies. The normalization is performed equivalent to the proton case with the fraction $\epsilon_{\gamma}$. With these power-law indices, the energy densities depend (at most) logarithmically on the maximal and minimal energies, reducing the impact of our choice of $\varepsilon^{\prime}_{\text{max}}$. Successful modeling of GRB prompt spectra within the internal shock model has been performed in e.g. Bosnjak et al. (2009); Daigne et al. (2011); Bošnjak & Daigne (2014). Realistic values for the synchrotron peak energy and spectral indices can be achieved, for instance, by allowing a small fraction of electrons to be accelerated to high energies (as in Eichler & Waxman (2005); Spitkovsky (2008)). However, an explicit modeling of photon synchrotron spectra is not within the scope of this work.

The coupled proton-neutron system is evolved to the steady state using NeuCosmA (Biehl et al., 2018), which takes photohadronic, pairproduction, adiabatic and synchrotron losses into account. The neutrons escape freely, whereas the protons are assumed to be magnetically confined and to escape only from the boundaries (within their Larmor radius), referred to as “direct escape” (Baerwald et al., 2013). This assumption implies that close to the photosphere the cosmic ray emission is dominated by neutron escape, and for large collision radii by direct proton escape (Bustamante et al., 2014). The mechanism for cosmic ray escape is currently discussed, and it is likely that several competing components contribute; see discussion in Zhang et al. (2018). The implementation for secondary particle emission is described in great detail in Biehl et al. (2018).

The model does not account for the emission from sub-photospheric collisions, for which the optical thickness to Thomson scattering is larger than one – resulting in a different, thermalized shape of the target photon spectra. Even if cosmic-ray acceleration could take place at such low radii, the high radiation densities will prevent the particles from reaching the UHE range. Since the pion production efficiency scales with the density similar to the Thomson optical depth, the neutrino production is most efficient close to the photosphere. There could be a significant contribution from sub-photospheric collisions; hence, our neutrino flux estimate shall be considered as minimal prediction. We chose examples in which the fraction of sub-photospheric collisions is small to reduce the impact of this effect.

Discussion of Reference model

We will use the format of Fig. 1 to characterize the behavior of different models, and hence explain it here in greater detail. Since the initial configuration of the 1000 shells is drawn from the distribution in Eq. (1), the model naturally produces stochastic fluctuations in the variables. We compute 100 representations if not otherwise noted. In the figures, we show the average as a solid curve and the region between the edge cases with a shaded band.

The left panel shows the differential energy dissipation in the engine frame for the different messengers: neutrinos, UHECRS above $10^{10}\,\mathrm{GeV}$, gamma-rays, and gamma-rays above $1\,\mathrm{GeV}$ only. In the dark-shaded region, all collisions occur below the photosphere, in the light-shaded region, the collision may be above or below the photosphere (depending on the other shell parameters). The result is consistent with earlier works (Bustamante et al., 2014, 2017): neutrinos originate near the photosphere due to the high density; UHECRs prefer intermediate collision radii since high magnetic fields are required for efficient acceleration – but not as high as that radiation losses (such as synchrotron radiation) limit the maximal energy; gamma-rays trace the region where most energy is dissipated (see middle panel and discussion below); HE gamma-rays above 1 GeV in the source frame prefer larger collision radii where maximal energies are not dominated by losses. The different astrophysical messengers originate from different radii of the same jet, thus drawing attention to the collision model that defines the connections among these different regions. The statistical fluctuations from the initial shell setup (shaded areas) imply some variability in this picture, but they do not change it qualitatively.

The middle panel of Fig. 1 shows the distribution of radius and dissipated energy of the individual collisions in arbitrary units. This figure demonstrates that most of the relevant collisions happen above the photosphere and there is a not very noticable anti-correlation between $R_{C}$ and $E_{\mathrm{diss}}$, confirming the relation to the gamma-ray output in the left panel.

The maximal cosmic ray energy in the upper right panel of Fig. 1 exhibits a similar shape as the UHECR curve in the left panel, that only accounts for cosmic rays with $E_{\text{p,max}}$ within the red-shaded area. The lower panel shows the magnetic field $B^{\prime}$ (left axis, black curve) and the optical thickness to photohadronic interactions evaluated at $E_{\text{p,max}}(R_{C})$ (right axis, green curve). Both scale with the collision radius, $\tau_{p\gamma}\sim R_{C}^{-2}$ and $B^{\prime}\sim R_{C}^{-1}$. The collisions close to the photosphere are optically thick to photohadronic interactions, and are responsible for most of the neutrino production (see left panel).

We define the overall energy dissipation efficiency from kinetic energy to radiation as

[TABLE]

In contrast to $\eta$ it describes the efficiency of the whole system and is an output instead of an input parameter. For the Reference model we find $\epsilon\approx 36\%$, i.e., somewhat higher than for the equal-energy setup in Bustamante et al. (2014, 2017). The reason for this is that the efficiency for individual collision $\eta$ is higher in the equal mass case (Kino et al., 2004). Compared to the equal-energy setup, the Reference model discussed here relatively efficiently converts the kinetic energy drawn from the GRB engine into secondaries. We discuss this issue in greater detail in Section (4) for alternative model assumptions, and we compare to the literature in Section (6.1). The kinetic energies assumed in our models (see Tab. 6a) are at the higher end of observations (Gruber et al., 2014), motivated by the higher baryonic loading required to reach the UHECR energy density observed at Earth.

3 Probability for two-shell final states

The key ingredient of the the Ultra Efficient model by Kobayashi & Sari (2001) is the emergence of two shells after a collision. Here, we study the probability of such events using ranges of collision parameters from stochastic GRB multi-collision models with hydrodynamical simulations, see Appendix (A) for details. We use the PLUTO code (Mignone et al., 2007) with a one-dimensional setup. A study for a two-dimensional setup demonstrated qualitatively similar behavior to one dimension (Mimica et al., 2004). Magnetic fields are neglected. The colliding shells are assumed to be cold. For a treatment of arbitrarily hot plasma shells see e.g. Pe’er et al. (2017), who find a possible suppression of the shock formation and a dependence of the energy per particle after the collision on the pre-collision plasma temperature.

The collision process and the post-collision shell configurations are illustrated in Fig. 2. Panel (a) demonstrates the equal energy case (in the source frame) and panel (b) the equal masses case. In both panels the shells are shown at $t^{\prime}=0$, the slow shell is depicted in blue, the fast one in red. The mass density profiles at $t^{\prime}_{\mathrm{shock}}$, when both shocks have crossed the respective shells, are shown in purple. In Fig. 2 (a), the resulting mass density profile is clearly single-peaked, while in Fig. 2 (b) two distinct peaks moving at different speeds can be seen. In order to identify two-shell post-collision configurations, we evaluate the mass density profile at the time when both shocks have crossed the shells and the internal energy is maximal. For a given snapshot in time, the mass density can exhibit multiple peaks; Kino et al. (2004) predict up to three peaks. Theoretical estimates comparing the time scales of the wave propagation (when the shock and rarefaction waves cross the two shells) may result in unrealistic approximations, since double-peaked profiles can rapidly evolve into single-peaked ones. After $t^{\prime}_{\mathrm{shock}}$, the density profile continues to evolve since the velocity profile is not uniform across the shell(s) (cf., lower panels of Fig. 2). Instead, the shell edges move in opposite directions in the CD frame. This leads to a dilution of the density profile as time evolves, invalidating the assumption of a constant shell width for the entire duration of a simulation (see also Pe’er et al. (2017)). As discussed below, the dissipation of internal energy reduces this effect by slowing down the thermal expansion that ultimately leads to a washed-out single-peaked profile.

For the classification of a two-shell collision, we define the relative depth of the dip between two mass peaks

[TABLE]

where $\rho_{\mathrm{Min}}$ is the density at the dip between the two maxima and $\rho_{\mathrm{Max}}$ the density at the lower peak (illustration see Fig. 2 (b)).

The main impact on the post-collision mass density profile comes from the pre-collision mass ratio ($m_{\mathrm{r}}/m_{\mathrm{s}}$), the Lorentz factors ($\Gamma_{\mathrm{r}}/\Gamma_{\mathrm{s}}$) and the shell widths ($l_{\mathrm{r}}/l_{\mathrm{s}}$). We fix the width ratio $l_{\mathrm{r}}/l_{\mathrm{s}}=1$ since in GRB internal shock models (see e.g. Kobayashi & Sari (2001); Kobayashi et al. (1997); Daigne & Mochkovitch (1998) and Globus et al. (2015); Bustamante et al. (2017); Bosnjak et al. (2009)), plasma shells are ejected at constant time intervals resulting in equal widths. An alternative choice $l_{\mathrm{r}}=0.1\ l_{\mathrm{s}}$ is discussed in the Appendix.

Fig. 3 shows the depth d as a function of $m_{\mathrm{r}}/m_{\mathrm{s}}$ and $\Gamma_{\mathrm{r}}/\Gamma_{\mathrm{s}}$ for the parameter ranges that enclose most of the collisions in our models (see Section (4)). As already noticed by Kino et al. (2004), pronounced dips, and hence double-peaked density profiles, occur for almost equal shell rest masses and for high Lorentz factor ratios $\Gamma_{\mathrm{r}}/\Gamma_{\mathrm{s}}$. The latter can be understood from the shock-/ rarefaction-wave timescales: If $\Gamma_{\mathrm{r}}/\Gamma_{\mathrm{s}}$ is high, then $l_{s}^{\prime}\ll l_{r}^{\prime}$. Therefore, the reverse shock takes substantially longer to cross the fast shell than the forward shock to cross the slow one. As a result, the rarefaction wave from the direction of the slow shell has enough time to create a pronounced dip separating the two shells.

The radiation model assumes that a fraction of the internal energy is converted into non-thermal electrons and/or ions which leave the system, which effectively cools the system. To study the impact on the collision dynamics we introduce a simplified energy dissipation term in our PLUTO simulations that removes internal energy from the system until a certain threshold has been reached (for more details see Appendix (A)). The rate of this process is assumed to be proportional to the available internal energy. The impact on the occurrence of two-shell configurations is demonstrated by the contours in Fig. 3 (b) and (c) for two efficiency choices. The dissipation of internal energy reduces the velocity and amplitude of the rarefaction waves, since those are powered from the available internal energy. Too slow rarefaction waves result in shallower dips in the mass density profile (see also Figure Fig. 7), decreasing $d$ and moving up the contours in Fig. 3.

For the Ultra Efficient shock scenario, this result means that the dissipation of internal energy into non-thermal particles significantly reduces the probability for two-shell final states. At the same time, the reduction of internal energy available to drive the rarefaction waves increases the lifetime of a two-shell configuration since the thermal expansion of a the double-peaked profile slows down. Hence, a higher energy dissipation rate reduces the emergence probability of two-shell configurations, but increases their lifetime. The model “PLUTO” described in the next section takes the first effect more rigorously into account.

4 Alternative collision models

In this section, we define three alternative models that use different assumptions for the collision of two shells, and study the impact on the multi-messenger production. The format of the discussion follows Section (2) where the details on the Reference model can be found, which is used as benchmark case.

The Reduced Efficiency model assumes that in each collision a fraction $\eta=0.5$ of the internal energy in Eq. (3) is dissipated as radiation, whereas the remaining kinetic energy goes into the adiabatic expansion of a single merged shell; the Ultra Efficient model assumes partially inelastic collisions in which the remaining internal energy is re-converted into kinetic energy, always resulting in two shells after the collision (Kobayashi & Sari, 2001); in the PLUTO model each two-shell collision is simulated individually with the PLUTO code.

Since each model would produce a different result for the same initial shell setup, we modify the setup of each model to re-produce comparative burst durations, variability times derived from the light curve, and total gamma-ray luminosities. Here, the variability timescale refers to the fastest time variability observed on top of longer-lasting light curve pulses. By construction, all GRBs are normalized to release the same amount of energy in photons in the optically thin regime. The burst duration ist determined by the initial size of the system (the sum of all initial shell widths and separations), which matches among the different models. The time variability primarily depends on the number of collisions for a constant burst duration, which scales with the number of initial shells. However, in the Ultra Efficient model shells do not merge when colliding, leading to substantially more collisions for the same number of initial shells. We compensate for this by reducing $N^{\mathrm{shells}}_{\mathrm{initial}}=1000$ to 125. While the resulting variability timescales are rather small, they do not necessarily translate into observed time variabilities, since these are limited by the instrument’s response (i.e., the actual variability could be smaller). We have verified that our results related to the multi-messenger production are not qualitatively affected by the choice of $1000$ initial of shells, see App. 11. In all cases we assume an engine with a constant mass outflow and constant up- and downtimes resulting in equal initial shell widths and separations. As previously discussed, a constant mass outflow is more likely to result in a splitting into two shells, which the Ultra Efficient model is based on. An overview of the model parameters is given in Tab. 4.