Neutrino and cosmic-ray emission from multiple internal shocks in gamma-ray bursts

Gamma-ray bursts are short-lived, luminous explosions at cosmological distances, thought to originate from relativistic jets launched at the deaths of massive stars. They are among the prime candidates to produce the observed cosmic rays at the highest energies. Recent neutrino data have, however, started to constrain this possibility in the simplest models with only one emission zone. In the classical theory of gamma-ray bursts, it is expected that particles are accelerated at mildly relativistic shocks generated by the collisions of material ejected from a central engine. We consider neutrino and cosmic-ray emission from multiple emission regions since these internal collisions must occur at very different radii, from below the photosphere all the way out to the circumburst medium, as a consequence of the efficient dissipation of kinetic energy. We demonstrate that the different messengers originate from different collision radii, which means that multi-messenger observations open windows for revealing the evolving GRB outflows.


INTRODUCTION
Gamma-ray bursts (GRBs) are violent outbreaks of energy, most of that being detected as gamma-rays during the so-called prompt phase (see Refs. [1][2][3] for reviews). The common view is that relativistic jets are ejected from a central engine, triggered by a collapsing star or a neutron star merger, in the direction of the observer. The inhomogeneity in jets naturally leads to internal shocks, at which charged particles can be accelerated. In the classical GRB scenario, the observed gamma-ray emission is attributed to synchrotron radiation from nonthermal electrons, which we however do not need to assume in our model. It is natural that protons are accelerated as well, and GRBs have also been considered as a possible candidate class for the origin of the ultra-high-energy cosmic rays (UHECRs) [4][5][6]. Since the charged cosmic rays cannot be traced back to their origin because of magnetic field deflections, neutrinos from GRBs, which are generated via proton-gas or proton-radiation interactions, can provide crucial clues to the UHECR mystery [7].
While PeV neutrinos from probably extragalactic sources have now been detected in the IceCube neutrino telescope [8], even approximately having the predicted spectral shape for GRBs [9][10][11][12], stacking searches for prompt GRB neutrinos using the timing and directional information have been so far unsuccessful [13,14]. Because some of the early analytical predictions of the GRB neutrino fluxes [14,15] have shortcomings which are independent of astrophysical uncertainty (although these do not exist in some numerical works such as Ref. [9]), the model used in Ref. [14] has been revised by about one order of magnitude [16][17][18]. The current data are even pushing into the expected regime of the latest predictions, enabling us to address whether GRBs can be the sources of the UHECRs, and what the neutrinos can tell us about that.
In most of the earlier discussions, the simple one-zone model is assumed in the sense that GRB parameters are fixed during the GRB duration. In particular, the internal shock radius R C , which is crucial for neutrino and UHECR production, is often estimated from geometric arguments [19] using the Lorentz factor Γ, variability timescale t v , and redshift z as The variability timescale can be obtained by inspection of the light curve; the Lorentz factor can only be estimated using various approaches [20,21]; and the redshift can be estimated via the observation of the host galaxy of the GRB. In the internal shock model, using the typical variability time (about three times shorter than the pulse width of ∼ 1 s [22,23]), R C ∼ 10 8 -10 10.5 km is expected [24] and neutrino predictions correspondingly vary [9,18,25]. In particular, in dissipative photospheric scenarios [10,[26][27][28], internal shocks may occur under or around the photosphere, at which the Thomson optical depth for eγ scattering [29] is unity. Since the photospheric radius R ph ∼ 10 8 -10 8.5 km is small, neutrino production is expected to be highly efficient. UHECR production is also sensitive to R C . There has been some evidence that the composition of UHE-CRs is heavy [30]. Heavy nuclei can be largely loaded in GRB jets [31,32], and acceleration to ultra-high energies and survival against photodisintegration are possible only if R C is large enough [25,33]. The UHECR escape also depends on GRB parameters [34]. Although it is often assumed that UHECRs can escape after the dynamical time, this is not the case if magnetic fields do not decay. Especially strong constraints on the UHECRneutrino connection can be obtained if cosmic rays escape only as neutrons which are produced in the same interactions as the neutrinos ("neutron-escape model") [35,36]. While this specific model is essentially excluded [37,38], a hard flux of protons leaking from the sources (further on called "direct escape" around the maximum energy and/or "diffusion escape" at lower energies) can dominate the UHECR emission, which is largely allowed by neutrino observations [34,37]. As demonstrated in Ref. [34], the dominant UHECR escape mechanism is in fact a function of the shell parameters.
Since the one-zone model is not realistic in the internal shock picture, R C and Γ should evolve even within one GRB. The R C -dependence of neutrino production efficiency has been discussed in the internal shock model [12,18,25,39], but its integrated effects have not been studied in detail. Hence, it is conceivable that not all collisions occur at the same radius, which has significant consequences for the neutrino and cosmic-ray production, as we show. For example, different UHECR escape mechanisms will dominate in different phases of the evolving GRB.

MODEL AND METHODS
In order to demonstrate neutrino and cosmic-ray emission from various R C , we follow the internal collision model of Ref. [40]; see Fig. 1 for illustration. We set out l l l l l The shells propagate, collide, and merge as soon as they meet other shells (multiple collisions are allowed), whereupon their widths and speeds change. The energy dissipated in the collision is assumed to be radiated away immediately. a number N sh of shells from a central engine with equal initial kinetic energies but a spread in the bulk Lorentz factor of where Γ k,0 is the initial Lorentz factor of the k-th shell and x follows a Gaussian distribution P (x)dx = e −x 2 /2 / √ 2πdx. Note that large A Γ 0.1 are required in order to achieve high efficiencies [41], as we have explicitly tested. The shells are assumed to be emitted with an uptime of the emitter δt, followed by an equally-long downtime, which is similar to the value of t v obtained from the simulation later. While the shells evolve, their widths, masses, and speeds (i.e., their Lorentz factors) are assumed to be constant, and their mass density decreases ∝ r −2 , with r the radial distance to the emitter. Because of the different Γ of the different shells, a shell will collide and merge as soon as it meets another one; see Fig. 1. We assume that the new shell immediately cools by prompt radiation of the internal energy into gamma-rays, cosmic rays, and neutrinos. Derivations of the properties of the newly formed shell are given in Refs. [40,42] and are maintained in the present system. Our results match the analytical predictions for the dissipation of modest-amplitude fluctuations from Refs. [43,44]. Note that we simplify the evolution of the internal shocks in several points, although our approach is enough for the purpose of this work. First, we assume situations where most of the dissipation occurs in the optically-thin regime. If significant dissipation occurs in the optically-thick regime, the internal energy scales adiabatically ∝ r −2/3 , which is spent to accelerate the outflow. Second, since we do not consider cases where only a fraction of the internal energy made available after a collision is released as radiation [41], this means the efficiency issue of the internal shock model may remain unresolved [45]. Third, observed light curves may have slow variability components as well as fast variability components [46], which cannot be easily explained by a discrete number of shells from a continuous emitter, whereas continuous outflow models give better agreements [47][48][49].
In this study, we choose Γ 0 = 500, N sh = 1000, δt = 0.01 s, and A Γ = 1 as well as a perfect acceleration efficiency of η = 1 (defined by t −1 acc ≡ η c 2 e B /E p , see Ref. [34]). The simulation yields 990 collisions, t v 0.06 s from the average obtained rise time of the light curve pulses (see Ref. [40]), a burst duration T N coll t v ≈ 59 s, and an average Γ ≈ 370 (Lorentz factor of the merged shells, corresponding to the observable Γ), i.e., the GRB is sufficiently close to conventional assumptions in neutrino production models. We normalize the total isotropic photon energy of all collisions in the source frame to E iso = 10 53 erg. Note that the fraction of photon energy dissipated in subphotospheric collisions is only about 9%, which means that a renormalization of the gamma-ray energy output to collisions above the photosphere would hardly affects our result. For the cosmicray and neutrino production, we follow Refs. [34,50] to compute the spectra for each collision individually, choosing equal energies in electrons (i.e., photons) and magnetic field, and a baryonic loading of ten (i.e., ten times more dissipated energy in protons than in photons). The target photon spectrum is assumed to be a broken power law with spectral indices α γ = 1 and β γ = 2, respectively, with a fixed break energy of γ,break = 1 keV in the merged-shell rest frame (primed quantities are in the merged-shell rest frame). That is, it is implied that the target photon spectrum corresponds to conventional GRB observations regardless of the underlying radiation processes leading to this spectral shape.
The light curve of the simulated burst is shown in Fig. 2, black curve. Although we show only one representative simulation in this study, we will present a more detailed parameter space study in a future work [51]. We do not investigate effects of the spectral evolution during the dynamical time for one collision [52], as we imply that taking into account contributions from multiple shells is more relevant as in the case of gamma-rays [42]. Note that, although we do not calculate hadronic cascades, their feedback on neutrino spectra is unimportant for the baryonic loading factors used in this work.

RESULTS
We show in Fig. 3 the neutrino fluence, maximal proton energy, and maximal gamma-ray escape energy for each collision (dot) as a function of R C . The maximal proton energies are obtained from comparing acceleration, dynamical, synchrotron loss, and photohadronic (for protons) timescales. As a result, we find that the collisions are spread between about 10 6 km and our choice of 5.5 · 10 11 km for the deceleration radius [53], where outflows terminate by the external shock, with most collisions around 10 10 km -slightly above the estimate from the geometry, R C ≈ 1.6 · 10 9 km. Red dots mark collisions in the neutron-escape model regime (optically thick to pγ interactions) and blue empty circles, collisions in the direct cosmic-ray escape regime.
Black squares mark subphotospheric collisions, i.e., those for which the Thomson optical depth is larger than unity. The optical depth is obtained by calculating the proton density from the masses of the shells and assuming that the electron density is as high as the proton density, which is expected for an electrically neutral plasma. In fact, the electron and positron densities may be somewhat higher if there is a significant non-thermal contribution from electron-positron pair production. The obtained photospheric radius R ph ≈ 2 · 10 8 km is somewhat larger than the conventional expectation calculated using the dissipated energy in gamma-rays (R ph ≈ 3 · 10 7 km). This estimate is affected by the efficiency of the conver- , maximal proton energy (in the source frame, for ideal (η = 1) acceleration), and maximal allowed gamma-ray energy (in the source frame, where τγγ(Eγ,max) = 1) as a function of the collision radius. Each dot represents one collision: red filled dots represent collisions where cosmic rays mainly escape as neutrons (optical thickness to pγ interactions larger than unity), blue empty circles represent collisions where cosmic-ray leakage dominates over the neutron-escape model, and black boxes denote subphotospheric collisions or collisions where this picture cannot be maintained (i.e., where the Thomson optical depth is large). In the right panel, the energy ranges which can be reached by the Fermi-GBM, Fermi-LAT, and Cerenkov Telecope Array (CTA) instruments are illustrated as colored bands. Collisions in which photons with energies above 10 6 GeV are able to escape are marked as upward-pointing arrows.
sion from kinetic energy into dissipated energy, which is roughly 25% in our cases. However, the more important reason is the significant baryonic loading: since most of the energy is dissipated into protons, the masses of the shells have to be upscaled to match the required energy output in gamma-rays (10 53 erg), which leads to larger radii of the photosphere because of higher electron densities. It can therefore be expected that the large baryonic loadings that are needed to describe the UHECR observations [25,37] will lead to larger fractions of subphotospheric collisions.
We find that the obtained range of collision radii is large, from below the photosphere out into the deceleration radius, since dissipation occurs for a wide range of R C especially when the spread of the Lorentz factor A Γ is large [43,44]. Note that 12% of collisions occurs below the photosphere for the chosen parameter set, altogether 118 out of the total 990 collisions, but most of the energy dissipation occurs at large radii > 10 10 km.
In the internal shock model, gamma-ray emission should be produced beyond the photosphere, so we only consider collisions beyond the photosphere in the following -unless noted otherwise. This is conservative, since the baryonic loading may be smaller below the photosphere [10] and particle acceleration becomes inefficient for radiation-mediated shocks [54].
Most importantly, Fig. 3, left panel, demonstrates that neutrinos are dominantly produced at small collision radii R C 10 9 km, close to the photospheric radius R ph ≈ 2 · 10 8 km. This result can be understood as follows. In each collision, the emitted gamma-ray energy, E iso γ−sh , is a fraction of the total dissipated energy. The pion production efficiency (fraction of energy going into pion production) at the spectral break depends on R C as where Γ m is the Lorentz factor of a merged shell. Since the internal shock model predicts [43] Hence, since A Γ has to be sufficiently large for efficient energy dissipation, neutrino production should be dominated by collisions at radii around the photosphere. UHECR protons come from collisions in the range 10 8.5 km R C 10 10 km; see middle panel of Fig. 3. First the maximal proton energy increases with collision radius as (close to the peak) synchrotron losses limit the maximal proton energy, and the magnetic field drops with R C . The peak occurs where adiabatic losses take over, and the decline comes from a decrease of the acceleration efficiency for dropping magnetic fields; the expressions for the different energy-gain and energy-loss timescales can be found, e.g., in Refs. [9,34]. Note that the UHECRs come from two different components dominating at different collision radii: for R C 10 8.5 km, neutron escape dominates, and for R C 10 8.5 km, protons directly escaping from the source dominate -which are obviously not related to strong neutrino production; see left panel. In the chosen example, the main contribution to cosmic rays actually comes from direct escape. Finally, the right panel of Fig. 3 illustrates that high gamma-ray energies, which can only be observed in CTA or other next-generation imaging atmospheric Cherenkov telescopes, come from large collision radii R C 10 9 km, since for lower radii the optical depth for γγ interactions is too high. As a consequence, neutrinos, cosmic rays, and Fermi-LAT/CTA gamma-rays probe different emission radii. Neutrinos are useful to probe dissipation at small radii, including subphotospheric dissipation. For dissipation at large radii, where heavy nuclei survive, the TeV gamma-ray diagnostics of a GRB would be useful [25]. In order to obtain an even more quantitative statement of how much energy is released as a function of collision radius, we show in Fig. 4 fractional binned distributions for the prompt gamma-rays, neutrinos, and cosmicray protons, which are all directly calculable within our model. We note that the energy per messenger per bin is obtained as a product of energy released per collision, and the number of collisions occuring per R C -bin; especially the latter number is important to get the proper weighing of R C . The result confirms the above observations: the neutrino production is dominated by small vales of R C just beyond the photosphere from within a relatively narrow region R C ≈ 10 8.5 to 10 9 km, the cosmic-ray production by intermediate R C ≈ 10 9 to 10 10 km, and the prompt gamma-ray emission is, in fact, dominated by large R C , at around 10 10 to 10 11.5 km -compatible with what is typically expected in the literature [24]. These results have significant implications: our knowledge of the prompt phase of GRBs is obtained from gamma-rays, of course, and, consequently, R C is derived from gamma-ray observations. This collision radius is, however, not the one to be used for neutrino or cosmic-ray calculations. It is therefore conceivable that multi-zone predictions are different from the naive one-zone expectation based on the gamma-ray emission radius. One can also read off from Fig. 4 that a significant amount of energy in UHE-CRs is transported away by direct escape, unrelated to neutrino production, which may affect the predicted neutrino flux if normalized to the observed UHECRs, as in, e.g., Ref. [37].
We show in Fig. 5 the predicted quasi-diffuse neutrino spectra from collisions beyond the photosphere as thick orange curves for three different values of Γ 0 , where the middle panel correspond to our standard assumptions. Note that the neutrino fluence per burst has been rescaled to a quasi-diffuse flux prediction by assuming 667 (identical) bursts per year and is significantly below the current diffuse neutrino signal reported by IceCube at the level of 10 −8 GeV cm −2 s −1 sr −1 flux [55]. The dashed curves correspond to the standard assumption that all collisions occur at the same radius, derived from gammaray observations. To generate these curves, we use the parameters N coll , t v , Γ , and T obtained from the simulation assuming identical shells with a collision radius obtained from Eq. (1) (R C ≈ 10 9.2 km in the middle panel). Note that the middle reference flux is significantly lower than the prediction in Ref. [16]. In that reference, the same parameters as in the IceCube analysis [56] were used for comparison, implying that R C ≈ 1.9 · 10 8 km. That is about one order of magnitude smaller than the R C used here; cf., Eq. (3) for its impact on the neutrino flux. The reference flux in the left panel is comparable to Ref. [16].
We first of all find that the neutrino spectra from collisions beyond the photosphere (thick orange curves) all exhibit the same flux level quite independently of Γ 0 (and even of A Γ , as we have explicitly tested). The expected neutrino flux per flavor is at the level of E 2 J ∼ 10 −11 GeV cm −2 sr −1 s −1 , peaking between 10 5 and 10 7 GeV. This contribution can be regarded as a minimal prediction for the neutrino flux, as it can be inferred from gamma-ray observations and hardly depends on the parameters. Note that this flux is probably outside the sensitivity of the existing IceCube experiment, but it should be well within the reach of the planned highenergy volume upgrade. There is a significant qualitative difference to conventional models such as Refs. [7,15], for which the pion production efficiency contains a factor Γ −4 coming from the collision radius estimate in Eq. (1) applied to Eq. (3). However, the optical thicknesses to Thomson scattering and photohadronic interactions both scale ∝ R −2 C , which leads to the following estimate for the pion production efficiency at the photosphere independent of Γ [10]: Here e is the fraction of the dissipated energy going into photons, and ε is the dissipation efficiency (ratio between dissipated and kinetic energies). Especially Γ drops out of f ph pγ -unless the break energy is fixed in the observer's frame, in which case there is a single factor of Γ retained. When the innermost collisions give the dominant contributions, the time-integrated neutrino fluence roughly scales as where N coll (f pγ 1) is the number of collisions with efficient neutrino production close to the photosphere, N tot coll is the total number of collisions, and p is the fraction of energy going into protons. Since the number of dominant collisions contributing to the neutrino flux is of order ten almost independently of the model parameters (see thin solid curves in Fig. 5), the neutrino flux prediction is relatively robust. The neutrino prediction above the photosphere hardly depends on the baryonic loading ( p / e ) as well, as long as most energy is dissipated into protons.
Increasing the baryonic loading in Eq. (5) is compensated by a correspondingly smaller e in Eq. (4). As a result, the neutrino flux is rougly independent of p / e -which we have explicitly tested numerically. This implies that the predicted flux E 2 J ∼ 10 −11 GeV cm −2 sr −1 s −1 can be hardly exceeded. The only exception may be increasing E iso , see Eq. (5), which may in fact be required to match the injected energy needed to describe UHECR observations; see Section 2 in Ref. [37].
The photon spectra can still be approximated by the Band function up to a Thomson optical thickness of ten or so [57,58], which occurs under the photosphere. This means that we can extrapolate our assumptions to under the photosphere to some degree. In the most extreme case, all energy may be dissipated into neutrinos, whereas the energy of neutrons, protons, and gammarays is reconverted -which is, however, very speculative, as nonthermal particle acceleration may not occur efficiently [10]. We show the corresponding subphotospheric extrapolations for the neutrino spectra as highly uncertain shaded regions in Fig. 5, corresponding to the black boxes in Fig. 3. Since the photospheric radius increases with decreasing Γ, the number of subphotospheric collisions increases with it, and their contribution in the left panel can be much higher than in the middle panel (and in the right panel much lower). As a consequence, the subphotospheric extrapolation may even reach the current sensitivity limit, and can be already constrained with current data. However, note again that this extrapolation is highly uncertain, as gamma-ray data cannot be used to obtain information from below the photosphere.
Finally, we show the "neutrino light curve" for our standard parameter set in Fig. 2 as a dotted (red) curve. It can be clearly seen that the neutrino flux is typically much lower than the gamma-ray flux except from some rare cases where the collision occurs close to the photosphere. Furthermore, the variation of the neutrino flux is larger due to the strong dependence of the pion production efficiency on R C .

SUMMARY AND CONCLUSIONS
In summary, we have studied neutrino, gamma-ray (at different energies), and cosmic-ray production in an evolving GRB outflow. We have demonstrated that they are produced at different collision radii. Consequently, the typical emission radius derived from prompt gammarays cannot be directly applied to neutrino and UHECR production, and the GRB will look very different from the point of view of different messengers. This concept is well known from conventional astronomical observations, where astrophysical objects look very different in different wavelength bands.
The neutrino spectra derived from gamma-ray observations are dominated by the emission close to the photosphere at R C ≈ 10 8.5 to 10 9 km, as the pion production efficiency depends on the collision radius in a non-linear way. UHECRs have been shown to be produced at intermediate collision radii R C ≈ 10 8.5 to 10 10 km, where the magnetic fields are high enough for efficient acceleration, but not so high that synchrotron losses limit the maximal proton energies. We have taken into account two possibilities for UHECR escape: emission as neutrons, which are not magnetically confined, and emission as protons from the edges of the shells -the dominant mechanism depends on the parameters of the colliding shells. Since the neutrons come from photohadronic interactions, their production dominates at smaller collision radii, whereas protons tend to be directly emitted at higher radii. The main energy in gamma-rays is deposited between around 10 10 to 10 11.5 km, compatible with earlier estimates. In particular, gamma-rays at the highest energies, such as in the energy range only accessible to CTA, cannot come from collision radii 10 9 km as the photon densities are too high there to let them escape.
For the quasi-diffuse neutrino flux prediction, we have identified two distinctive contributions. Above the photosphere, gamma-ray observations can be used to infer the pion production efficiency, which leads to a neutrino flux per flavor E 2 J ∼ 10 −11 GeV cm −2 sr −1 s −1 for an assumed isotropic energy of 10 53 erg emitted in gammarays. Especially, there is no strong dependence on the Lorentz boost Γ, in contrast to conventional one-zone models, as both the photosphere and the pion production efficiency scale with the collision radius in the same way. This is the minimal neutrino flux which one would at least expect in stacking analyses based on the actual gamma-ray observations, such as Ref. [14]. There is also a significantly milder dependence on the baryonic loading, as this parameter changes the photosphere of the model at the same time that it rescales the neutrino flux. The flux is significantly lower than earlier predictions [16] because a) we have explicitly excluded subphotospheric contributions, b) large photospheric radii have been obtained as a consequence of significant baryonic loadings (ten) and the moderate energy dissipation efficiency of the fireball (25%), and c) only a small number of collisions beyond the photosphere occurs at radii where the neutrino production efficiency is high. This expected "minimal" flux is beyond the sensitivity of the current IceCube experiment, but could be reached in future high-energy extensions. No gamma-ray information from deep under the photosphere can be directly obtained, and the neutrino production in that regime is more speculative [54]. In principle, a neutrino flux at the current sensitivity limit is conceivable, which means that current neutrino observations can already constrain the subphotospheric neutrino production for sufficiently large baryonic loading factors.
Our results imply that model-dependent studies of the multi-messenger connection, such as a GRB stacking analysis of neutrino fluences, can be improved, and give a stronger case for testing the hypothesis that UHECRs originate from GRBs. Compared to the one-zone model, some additional assumptions need to be made for the distribution of the collision radii. In particular, the width of the initial distribution of the bulk Lorentz factor A Γ , with which the shells are set out by the central engine, turns out to be the key additional parameter. However, it can in principle be obtained from comparing the light curves between simulation and observation. On the other hand, we have the advantage that the uncertainty in R C , which is the key issue in the standard model, disappears, as a collision radius distribution is predicted by the theory. There should also be new opportunities from our observations: different messengers can be used to study different regions of an evolving GRB outflow. For instance, direct neutrino and gamma-ray observations in CTA of a single GRB would open windows to very different regions of the GRB.
Note added. During completion of this work, Ref. [59] appeared. Although it shares some common aspects, we explicitly study the emission from all zones, compared to snapshots of the GRB evolution in Ref. [59].
We thank V. Mangano, E. Waxman, and B. Zhang for discussion and comments. This work is supported by NASA through Hubble Fellowship Grant No. 51310.01 awarded by the STScI, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under Contract No. NAS 5-26555 (K.M.). M.B. and W.W. would also like to acknowledge sup-port from DFG grants WI 2639/3-1 and WI 2639/4-1, and the "Helmholtz Alliance for Astroparticle Physics HAP", funded by the Initiative and Networking Fund of the Helmholtz Association. P.B. acknowledges support from NASA grant NNX13AH50G.