On the hadronic origin of the TeV radiation from GRB 190114C

The recently discovered TeV emission from Gamma-Ray Bursts (GRBs) has renewed the long-standing discussion about the hadronic versus leptonic origin of the observed GRB radiation. In this work, we investigate the possibility that the very high energy gamma rays observed by MAGIC from GRB 190114C (with energy from ∼0.1 to ∼0.8 TeV) are originated in a hadronic model. We developed a Monte Carlo (MC) simulation of the source internal state dynamics and of the photo-hadronic interactions at internal shock. We included in the simulation also the pair production process that the secondary gamma rays undergo in the GRB jet. We find upper limits on the internal shock model parameters by comparing our simulations to the sub-TeV observations of GRB 190114C. Neutrino flux predictions by the model are found to be consistent with experimental upper limits set by ANTARES and IceCube.


Introduction
GRB 190114C is a long-duration GRB observed to emit gamma rays in the TeV band.MAGIC detected, from about one minute after the burst, high-energy gamma rays (∼0.2 to ∼0.8 TeV) with high statistical significance, at the transition between the prompt and afterglow phases of the GRB emission (MAGIC Collab. 2019a).Several models have been proposed in order to reproduce the broadband emission observed in GRB 190114C, both hadronic (Melandri at al. 2022;Sahu & López Fortín 2020) and leptonic (Chand et al. 2020;Derishev & Piran 2019;Fraija et al. 2019b;Ravasio et al. 2019;Rueda et al. 2020;Wang et al. 2019).Although, firm conclusions on the production mechanisms of GeV-TeV emission have not been reached so far, being also limited by the large number of parameters involved in GRB modelling.
In this paper, we consider the hypothesis that part of the emission at high energy from GRB 190114C could be caused by the presence of a hadronic component, with the subsequent production of high-energy neutrinos (Eichler, Guetta & Pohl 2020;Guetta, Spada & Waxman 2001;Guetta 2015;Yacobi et al. 2014).
The GRB 190114C light curve exhibits at the beginning irregular multi-peaks due to bunches of γ rays superimposed on a smoothly varying emission component that extends beyond the highly variable emission period.Therefore observations do not exclude that part of the high energy emission may be released during the prompt phase.We consider the possibility that a fraction of the Very High-Energy (VHE) emission of GRB 190114C may be due to photomeson interactions within the internal shocks (IS) region.We consider the standard fireball model in which energy dissipation occurs at IS between shells in the relativistic outflow of the jet or through interactions with ambient matter (Eichler, Guetta & Pohl 2020;Guetta, Spada & Waxman 2001).Thus, a substantial part of the bulk kinetic energy is converted into internal energy, which is then distributed between electrons, protons, and magnetic field.The internally accelerated electrons are presumably responsible for the keV-GeV photons observed in the GRB, which are emitted through synchrotron or inverse Compton processes.Accelerated protons may interact with these (∼ MeV) photons and produce both neutral and charged pions, which in turn decay into high-energy photons and neutrinos, respectively.However, MAGIC observations did not show any sign of significant spectral change or sharp flux variations from prompt phase to prolonged emission that are characteristic features of the prompt phase.Furthermore, recent LHAASO observations (LHAASO Collab.2022) provide a compelling evidence that the TeV emission can come from the afterglow, even if it overlaps in time with the prompt emission.MAGIC observations started in the early afterglow phase, therefore several authors have interpreted this emission as due to the external shock model Derishev & Piran (2019); Sahu & López Fortín (2020).In this paper we consider the possibility that part of the TeV emission is due to the internal shocks.The comparison between MAGIC data and the results of the Monte Carlo (MC) simulation developed for this work can set an upper limit on the contribution to very high energy photons due to the prompt phase.The results described in this paper have been obtained by a full Monte Carlo simulation that described in detail in a previous paper (Fasano et al. 2021).For this simulation each accelerated proton is generated according to a power law (E −2 ) and is tracked into a radiation field until it interacts or escapes.The radiation field is generated according to the Fermi measurements (Fermi-GBM Collab.2019) as described in Sec. 3. Photo-mesons interactions are simulated for the production of charged and neutral pions.Photons and neutrinos are so obtained from mesons decays.For each photon from pion decay we evaluate the probability to escape from the source, if the photon interacts we evaluate the probability that a secondary photon can emerge from the source.This makes our simulation different from previous works i.e. (Sahu & López Fortín 2020).The paper is structured as follows: in Sec. 2, we describe spectral and temporal properties of GRB 190114C, assumed for the simulation.In Sec. 3, we present our MC program, developed to simulate the photo-hadronic interactions occurring during the prompt phase.We also describe the result of the simulation of the electromagnetic cascades initiated by the interaction of high energy photons produced in π 0 decays.In Sec. 4, we compare the photon flux resulting from our simulation to MAGIC data, deconvolved for the Extragalactic Background Light (EBL) interactions.From this comparison we obtain a set of best fit values for our model.We then predict a flux of high energy neutrinos and evaluate the expected number of events in present and future experiments.A discussion of our results is provided in Sec. 5.

GRB 190114C: spectral and temporal properties
The first detection of GRB 190114C is due to Swift (Swift Collab.2019), further observations have followed by GBM (Fermi-GBM Collab. 2019) Collab. 2019a).Within this temporal window, the time-dependent analysis of VHE data showed a systematic decrease in flux normalization, as well as a steepening trend over time.The high-energy photon spectrum reported by MAGIC, for the time interval 68-110 s, entirely overlaps with the T 90 estimated by Fermi -GBM for the prompt emission.During this time interval, the intrinsic burst spectrum in the 0.2-1 TeV band is characterised by a pure power-law (∝ E −ξ ) with ξ = 2.16 +0.29 −0.31 (MAGIC Collab.2019b).The photon spectral slope at the source has been derived deconvolving, from the observed spectrum the severe attenuation of the gamma ray flux due to its propagation to the Earth within the EBL (MAGIC Collab.2019b), according to the Dominguez et al. model (Dominguez et al. 2011).We hence consider such an intrinsic spectrum and compare it to the prediction of the gamma-ray flux emerging from our simulation of phenomena happening in the IS region.Additionally, we investigated the effects of adopting a different EBL model, e.g. the one from Franceschini et al. (Franceschini & Rodighiero 2017, 2018), finding results consistent with the Dominguez et al. model within the statistical uncertainty of the MAGIC measurements.The temporal overlap between the MAGIC observations and T 90 is not sufficient to definitely attribute the high-energy photons measured by MAGIC to the prompt or to the afterglow phase of the GRB emission.In this paper we assume the hadronic scenario and we compare the MAGIC observation with the radiation emerging from photomeson interactions where accelerated protons interact with the Band-like target radiation field in the prompt phase.The two main features of the hadronic scenario are the bulk Lorentz factor Γ of the relativistic jet, and the amount of energy channeled into relativistic protons E iso,p .This last quantity can be expressed with the baryon loading f p = E iso,p /E iso .Both the values of Γ and f p can be constrained to reproduce the VHE MAGIC observed spectrum.

Monte Carlo simulation
Modeling the physical processes occurring inside the IS region of the expanding GRB fireball requires characterization of the site where particles propagation and interactions take place.
Here we consider a simplified stationary one-zone scenario (Asano et al. 2009;Hummer et al. 2012;Murase & Nagatki 2006) in which mildly relativistic shells of plasma collide at a typical radius (Bearwald et al. 2015) where all the GRB energy is released.As variability timescale, we assume the value t var = 6 ms as suggested by observations during the prompt phase of GRB 190114C (Ajello et al. 2020).The bulk Lorentz factor Γ is treated as a free parameter of the model, and it will be fixed as the one that best reproduces MAGIC data.
The MC calculation is performed in the IS frame, assuming a spherical geometry (Baerwald, Hummer & Winter 2012) and a shell width ∆R IS = Γct var /(1 + z).We simulate a flux of accelerated protons with energies ranging from 1 GeV to 10 9 GeV, according to dN p /dE p ∝ E −2 p , consistently with a Fermi I-order acceleration process.Such a large energy range has been selected in order to avoid biases in the results.The target photon energy distribution dn γ /dϵ γ in the IS frame is assumed to reproduce the Band function observed by Fermi -GBM (Fermi-GBM Collab.2019) in the prompt phase.In order to allow the particle ∆ + resonance1 production with the entire spectrum of accelerated protons, the Band function has been extended to high energy following the trend of the higher part of the spectrum.Each generated proton can either interact with ambient photons, if the center of mass energy is above the interaction threshold condition, or it propagates further in the shell.The average interaction length λ pγ (s) = [n γ σ pγ (s)] −1 is evaluated, depending on the pγ center of mass energy and on the density of photons in the IS frame n γ (ϵ IS th,∆ ), requiring a center of mass energy above the ∆ + production threshold: , m ∆ and m p being respectively the ∆ + and proton masses.
According to the average interaction lenght λ pγ (s) we extracted the proton path before the interaction x p .If x p < R IS and if the threshold condition for the ∆ + production is satisfied, the photo-meson interaction occurs.The energies of protons (E IS p ) and target photons (ϵ IS ph ) that satisfy the photo-meson production condition are shown in Fig. 1.
Target photon energy ϵ IS ph vs proton energy E IS p for events above the ∆ + production threshold.Energies are given in the IS frame (Γ = 800 and t var = 6 ms).The colour code indicates the number of simulated events.
Multi-pion generation beyond resonant ∆ + production is also simulated and secondary particles like photons, muons, and neutrinos originated by decays are followed.Secondary protons, from ∆ + decay, are tracked until they leave the IS region to account for possible re-interactions along their path.
To account for proton energy losses during the propagation inside the IS region, we compare the proton acceleration timescale t acc (E IS p ) = r L (E IS p )/c (r L being the particle Larmor radius), to the average ∆ + production collision timescale t pγ (E IS p ) = [n γ (ϵ IS th,∆ )cσ ∆ pγ K pγ ] −1 .In the previous expression, K pγ = 0.13 and σ ∆ pγ = 5 × 10 −28 cm 2 are respectively the inelasticity coefficient and the cross section for the interaction at the threshold (Atoyan & Dermer 2001).
Charged particles may be affected by synchrotron losses, therefore, in order to evaluate these effects, we calculate the magnetic field value at IS by equating the magnetic energy density U B = B 2 /(8π) to the kinetic energy density of the accelerated electrons.Naming ϵ B and ϵ e respectively the fraction of jet energy converted into magnetic field and carried by electrons, we assume ϵ e = ϵ B = 0.1 (Waxman 2000).Proton synchrotron losses are relevant for energy above E IS p,cut ≃ 45, 65, 100, 400 PeV assuming Γ = 300, 500, 800, 1000 and t var = 6 ms.These evaluations confirm our choice to simulate the proton energy up to 10 9 GeV.Pion and muon propagations into magnetic field are also affected by synchrotron energy losses.We have compared the particles lifetime with the synchrotron timescale, t syn (E IS ) = 3m 4 x c 3 /4σ T m 2 e ϵ π U B , (Guetta & Granot 2003), where m x is the mass of the particle which we are taking into account.We found E IS µ ± ,cut ≃ 0.5, 2.6, 10, 20 TeV and E IS π ± ,cut ≃ 10, 50, 200, 400 TeV (Γ = 300, 500, 800, 1000, t var = 6 ms).The MC simulation of pγ interactions follows the methods described in (Fasano et al. 2021).The energy particle distribution for interacting protons, charged pions, muons and neutrinos assuming Γ = 800 and t var = 6 ms is shown in Fig. 2.
Concerning the neutral pion production and decay, we consider the possibility that each originated gamma ray might escape from the IS region, or interact with target photons, via e ± pair production.The average radiation length for this process is evaluated following (Waxman 2003) as with the condition ϵ th E IS γ (1 − cos θ) > 2(m e c 2 ) 2 , being 0 < θ < π the angle between the photons in the center of mass reference frame, and σ γγ the e ± pair production cross section (Vernetto & Lipari 2016).For the target photons we assumed the spectral energy distribution dn γ /dϵ γ approximated by a Band function as described in Sec. 2. For each π 0 -originated photon with energy E IS γ , we extract its free path x γ according to the average radiation length λ γγ in Eq. (3.2).If x γ > R IS the photon escapes the IS region and will eventually be observed in the laboratory frame with energy E obs = ΓE IS /(1 + z).Most of the photons originated by π 0 decays are absorbed and initiate an electromagnetic cascade.We track the products of the interactions, until their energy can be included within the MAGIC observed energy range.
The contribution of photons from electromagnetic cascades to the MAGIC observed spectrum depends on the bulk Lorentz factor value used in the simulation: for a given π 0 energy this contribution from cascades is larger for lower values of Γ.The flux of escaping photons can be directly compared with the intrinsic spectral energy density evaluated on the basis of MAGIC published results (MAGIC Collab. 2019b).The contribution of photons from electromagnetic cascades has been considered previously for GRBs (Wang 2018) and other astrophysical sources (i.e.blazars (Bottcher 2013)).For our work, instead of the semi-analytical treatment described in (Kelner and Aharonian 2008), we have developed a full MC simulation and applied it to GRB190114C.

Results: gamma rays and neutrinos from photo-hadronic interactions
We perform MC simulation for different values of t var and Γ, obtaining spectra of highenergy gamma rays and neutrinos emerging from the interaction region.We then convert these spectra into expected fluxes on Earth by taking into account the cosmological nature of GRBs.We evaluate the quantity E 2 γ ϕ(E γ ) on Earth rescaling the IS energy spectrum according to: i) the luminosity distance d L , ii) the first time interval of MAGIC observations ∆t = 42 s, and iii) the energy conversion factor from the IS to the observed frame, as: For each simulated Γ and t var , the baryon loading value f p has been assumed in the interval 1 − 30 to find the best agreement between the EBL-deconvolved MAGIC data, in the 68-110 s time interval (MAGIC Collab. 2019b) and the simulated fluxes, through a χ 2 statistical test.For the assumed values of Γ = 300, 500, 800, 1000 we obtain f p = 23 ± 6, 13 ± 3, 11 ± 2, 7 ± 2 as best-fit parameters.These values are well consistent with the expectations (Murase 2007).
For each Γ value, we evaluated the expected flux on Earth E 2 γ Φ(E γ ).This result is shown in Fig. 3 where the shaded area represents the 1σ statistical uncertainty on f p .For the assumed values of Γ, we found that the majority of photons emerging from the source is originated in the electromagnetic cascades due to interactions of photons from π 0 .A direct proof of the hadronic origin of the observed TeV radiation might come from coincident neutrino observations.Because of the transient nature of GRBs, the detection of a single neutrino event would allow identifying these sources as extreme hadronic accelerators.
To evaluate the neutrino energy flux on Earth for muon neutrinos and antineutrinos we assume Γ = 800, we use the energy conversion factor described in Eq.( 4.1), and we take into account neutrino oscillations, assuming normal ordering and the standard three-flavor scenario (Mascaretti & Vissani 2019).To check whether present and future neutrino detectors might be able to find a signal from an astrophysical source similar to a GRB 190114C, we have evaluated the expected neutrino, plus antineutrino, fluence dN νµ+ νµ /dE νµ+ νµ dS.We can then compute the expected number of track-like events with the where A νµ eff (E νµ , δ) is the Neutrino Telescope effective area, given as a function of the neutrino energy and of the source declination δ.Public effective areas of ANTARES (ANTARES Collab.2012) and IceCube (IceCube Collab.2014) are evaluated at analysis level and for the declination band that includes the position of GRB 190114C.For KM3NeT/ARCA the effective area is evaluated at trigger level and averaged throughout the sky (KM3NeT Collab.2016).In Fig. 4 we report the number of muon neutrino and antineutrino expected events for the three neutrino detectors as a function of the neutrino energy.In Table 1, we provide the total number of neutrino-induced events expected in the three detectors, obtained by the simulation with t var = 6 ms, Γ = 800 and f p = 11.Both the ANTARES and IceCube Collaborations have unsuccessfully searched for coincident neutrino-induced signals from the direction of GRB 190114C.In the case of ANTARES, the derived 90% confidence level integrated upper limit amounts to 1.6 GeV/cm 2 .For IceCube the same limit amount to 0.44 GeV/cm 2 (ANTARES Collab.2021; IceCube Collab.2019).In both cases, the constraints are limited to the muon neutrino component reaching Earth.In fact, because angular precision is a crucial feature in reducing the atmospheric background entering the search cone angle, muon neutrino charged-current interactions constitute the better astronomical channel, as muons emerging from these can be identified in neutrino telescopes as long tracks.The non-detection of current instruments is compatible with the hadronic model expectations for the values of f p and Γ that better reproduce the sub-TeV gamma-ray MAGIC data (IceCube Collab.2021).

Detector
Declination band  Since the coincident background rate expected from atmospheric neutrinos is by orders of magnitude lower than the signal rate expected from GRB 190114C, a detection from a single GRB appears to be very difficult.This result suggests that stacking several GRBs of the same kind could be a viable solution for a significant detection.Despite the experimental efforts, no clear signal has emerged so far in data from the ANTARES and IceCube telescopes (IceCube Collab.2014; ANTARES Collab.2012) when searching for spatial and temporal coincidences with the prompt emission of a sample of stacked GRBs.This leads to a constraint on the possible contribution of standard GRBs to less than 10% of the diffuse cosmic neutrino flux (IceCube Collab. 2021;ANTARES Collab. 2013, 2017, 2021).Furthermore, it is still unclear whether TeV-emitting GRBs behave as the standard GRB population.

Discussion and Conclusions
The successful observation of radiation from GRBs extending up to the TeV domain has provided further evidence of the extreme nature of these sources.However, there is still poor understanding of the processes that characterize these TeV emissions that might witness the presence of effective hadronic acceleration in GRB jets, possibly already during the prompt phase of the emission or, more realistically, during the afterglow phase.This occurrence would establish the connection between GRBs and Ultra-High-Energy Cosmic Rays (UHECRs), which remains a long-standing paradigm still to be proven.
To test the hypothesis of the hadronic origin of TeV radiation, we developed a MC simulation of photo-meson interactions between high-energy protons and target photons distributed according to a Band-like spectrum, as indicated by Fermi-GBM observations, during the prompt phase.After computing the spectra of secondary particles emerging from these interactions, we additionally simulated the electromagnetic absorption that gamma rays undergo in the IS shell.The spectrum of escaping photons so obtained has been compared to the intrinsic source spectrum derived by the MAGIC observations of GRB 190114C.We obtained our results by developing a full MC simulation of the interactions of high energy accelerated protons with the IS radiation.High energy photons from π 0 , are followed until they interact in the IS region or escape.The total intrinsic flux of photons, directly from π 0 or from electromagnetic cascades, escaping from the interacting regions is then evaluated.This flux transported to Earth is compared with MAGIC observations deconvolved by EBL interaction.Several simulations have been performed, in the framework of both IS scenarios, assuming different values of the relevant model parameters, like the Lorentz factor Γ and the baryon loading f p .From the comparison of the simulation results and the MAGIC data, we extracted the parameter values that better reproduce the observations.In this paper we have considered the possibility that part of the TeV emission may be due to the IS model responsible of the prompt phase.We have found the best fit parameters for this model that can reproduce the maximum contribution due to the IS interactions.
The same MC code has been used to evaluate the flux of neutrino from charged pion decay arriving at Earth.Applying the same MC simulation to the entire sample of the observed TeV GRBs could provide insights on the physical mechanisms responsible for the TeV emission.A better knowledge of the high energy photon spectrum, that at present seems to be limited by the late response of imaging atmospheric Cherenkov telescopes in pointing, would help in the GRB characterization.Confirmation of the hadronic origin of sub-TeV radiation is expected from neutrino observations.However our simulation indicates that such a detection from individual astrophysical sources like GRB 190114C appears extremely unrealistic, as confirmed by the lack of spatial correlations in data reported by the ANTARES and IceCube neutrino telescopes.

Figure 2 .
Figure 2. Particle spectra resulting from pγ interactions, obtained in the IS frame from the simulation with Γ = 800 and t var = 6 ms.

Figure 3 .
Figure 3. IS simulation: comparison between the EBL-deconvolved flux of GRB 190114C, as provided by MAGIC in the temporal interval 68-110 s (MAGIC Collab.2019b) and the MC results.Photon fluxes are due to π 0 -decays and to electromagnetic showers following internal gamma-ray absorption on the IS for different parameters of the model.

Table 1 .
Number of neutrino induced events expected by the GRB Internal Shock MC simulation for different neutrino telescopes, as due to ν µ and νµ interactions during the [68;100] s time interval of the prompt emission of GRB 190114C, for the model with t var = 6 ms, Γ = 800 and f p = 11.