Subphotospheric Emission from Short Gamma-Ray Bursts: Protons Mold the Multimessenger Signals

The origin of the observed Band-like photon spectrum in short gamma-ray bursts (sGRBs) is a long-standing mystery. We carry out the first general relativistic magnetohydrodynamic simulation of an sGRB jet with initial magnetization σ 0 = 150 in dynamical ejecta from a binary merger. From this simulation, we identify regions along the jet of efficient energy dissipation due to magnetic reconnection and collisionless subshocks. Taking into account electron and proton acceleration processes, we solve for the first time the coupled transport equations for photons, electrons, protons, neutrinos, and intermediate particle species up to close to the photosphere (i.e., up to 1 × 1012 cm), accounting for all relevant radiative and cooling processes. We find that the subphotospheric multimessenger signals carry strong signatures of the hadronic interactions and their resulting particle cascades. Importantly, the spectral energy distribution of photons is significantly distorted with respect to the Wien one, commonly assumed below the photosphere. Our findings suggest that the bulk of the nonthermal photon spectrum observed in sGRBs can stem from hadronic processes occurring below the photosphere and previously neglected, with an accompanying energy flux of neutrinos peaking in the GeV energy range.

The observed non-thermal nature of sGRB electromagnetic spectra (Ghirlanda et al. 2004;Burgess et al. 2019b;Poolakkil et al. 2021) hints towards efficient particle acceleration taking place in these jets (see e.g.Rees & Mészáros 1994;Spitkovsky 2008, for the related theoretical background).However, the physical processes driving particle acceleration in the sGRB prompt phase and the related gamma-ray emission are unset-tled.Specifically, sGRB observations are difficult to reconcile with optically thin synchrotron models (Mészáros et al. 2001;Ajello et al. 2019;Li 2020;Ryde et al. 2022, see however also Burgess et al. (2019a)).An alternative is provided by photospheric models of (modified) thermal spectra released as the outflow becomes optically thin (see e.g.Ramirez-Ruiz 2005;Pe'er & Ryde 2017;Beloborodov & Mészáros 2017).Such models may be especially interesting for sGRBs, which generally have harder low-energy spectral index compared to long GRBs (Poolakkil et al. 2021;Acuner & Ryde 2018;Dereli-Bégué et al. 2020).
In order to obtain gamma-ray spectra compatible with observations, dissipative effects below the photosphere are required.Among the typically invoked processes are synchrotron and inverse Compton radiation from a non-thermal lepton population (Pe'er et al. 2006;Giannios 2008;Vurm et al. 2013;Vurm & Beloborodov 2016;Bhattacharya & Kumar 2020), radiation-mediated shocks (Bromberg et al. 2011;Lundman & Beloborodov 2021;Samuelsson & Ryde 2022) and proton-neutron (pn) collisions induced as neutrons drift with respect to the proton flow in the jet (Beloborodov 2010;Vurm et al. 2011).While much work on these scenarios has been carried out relying on analytical outflow models, the effects of dissipative processes on the photon spectra remain to be demonstrated in consistent numerical simulations (e.g. as performed with a Monte-Carlo approach in Ito et al. 2015Ito et al. , 2021Ito et al. , 2023;;Parsotan et al. 2018).In addition, below the photosphere, the photon spectrum could be further modified by lepto-hadronic and hadronic processes (Kelner et al. 2006;Kelner & Aharonian 2008), but these effects have been largely neglected.
If either dissipative processes co-accelerate baryons contained in the jet or heating occurs due to pn collisions, production of high-energy neutrinos is also expected.Such sub-photospheric neutrino signals were previously calculated based on simple models (Razzaque et al. 2003;Ando & Beacom 2005;Tamborra & Ando 2016;Murase 2008;Wang & Dai 2009;Gao et al. 2012;Asano & Mészáros 2013;Murase et al. 2013;Xiao et al. 2017).More recently, relying on GR-MHD simulations of the jet, neutrino signals at lower energies (≲ 10 5 GeV) were predicted both for collapsars (Guarini et al. 2023) and sGRBs (Gottlieb & Globus 2021).In fact, because of the strong mixing with the cocoon (obviously not taken into account in naïve models), the jet is loaded with baryons, and neutrino production is hindered in the deepest jet regions (Guarini et al. 2023;Gottlieb & Globus 2021).Such recent developments call for advanced modeling of the particle production jointly with the jet evolution in order to forecast the expected multimessenger signals reliably.Gottlieb & Globus (2021); Guarini et al. (2023), however, focused on the jet region up to several 10 10 cm, thus well below the photosphere (and before breakout from the dynamical ejecta in Gottlieb & Globus (2021)).At these radii, a thermal (Planck or Wien) spectrum was assumed (Bégué & Pe'er 2015;Chhotray & Lazzati 2015;Gottlieb et al. 2019).However, as the outflow approaches the photosphere, non-thermal processes may leave imprints both on the photon and neutrino distributions.Hence, it becomes necessary to self-consistently calculate the coupled evolution of photons, protons, leptons and intermediate species.Such non-thermal processes may also have a non-negligible (yet unknown) impact on the observable signals.
This Letter explores the multi-messenger signals produced in sGRBs up to close to the photosphere, including leptonic, hadronic as well as lepto-hadronic processes, relying on the benchmark GR-MHD sGRB jet simulation introduced in § 2. The regions of efficient energy dissipation along the jet are outlined in § 3.After introducing the coupled transport equations for all particles in § 4, we present our results on the temporal evolution of the spectral energy distributions of photons and neutrinos.A discussion on our findings follows in § 5, and we conclude in § 6.Additional characteristic jet quantities are provided in Appendix A, while details on the numerical evolution of the particle distributions are outlined in Appendix B.

SHORT GAMMA-RAY BURST JET MODEL
We perform a GR-MHD simulation of a black hole (BH)-powered jet in homologously expanding merger ejecta using the GPU-accelerated code h-amr (Liska et al. 2022).We reproduce the fiducial simulation of Gottlieb et al. (2022), α3d5, with the only differences between our and their simulation being that here the jet is launched with a high magnetization of σ 0 = 150 (leading to an asymptotic Lorentz factor of O(10-20) compatible with those expected in sGRB jets (e.g.Zou et al. 2018)) and evolves for longer times and larger distances.The setup assumes that the binary merger remnant is a meta-stable neutron star that lasts for t d = 0.5 s before collapsing into a BH with mass M BH = 3M ⊙ and dimensionless spin a = 0.9375.During the metastable NS lifetime, a torus with mass M t = 0.2M ⊙ and a characteristic gas-to-magnetic pressure ratio of 1000 forms around the compact object.Additionally, the merger debris reaches homologous expansion between velocities v min = 0.05c and v max = 0.25c.Thus, the simulation is initiated with such torus around the BH, Figure 1.Top: Charateristic properties of our benchmark jet simulation at 7 s and in the x-z plane.From left to right we show the radial component of the Lorentz factor Γ rad , the logarithmic comoving mass density log 10 (ρ ′ ), and the logarithmic magnetization log 10 (σ).In order to highlight the location of the relativistic jet, we plot here the jet region with viewing angle −0.4 rad < θ < 0.4 rad; the dotted white line marks θ = 0 to guide the eye.The blue, purple, red, and yellow isocontour lines correspond to the radial Lorentz factor Γ rad equal to 1.5, 3, 10, and 25.At 7 s, the relativistic jet sits around 1.4-2 × 10 11 cm and it is surrounded by a mildly relativistic cocoon, whose comoving mass density (magnetization) is larger (smaller) than that of the jet.The jet simulation inputs from the shaded region at r ≳ 1.8 × 10 11 cm for the 7 s snapshot are not considered in our investigation of the particle acceleration sites.Instead, an extrapolation procedure based on a comoving shell located between 1.48-1.8× 10 11 cm for the 7 s snapshot is adopted (see main text for details and region delimited by the white solid lines in the top panels); we then extend such extrapolation up to ≃ 10 12 cm, which is slightly below the photosphere, i.e. beyond the jet evolution computed through the GR-MHD simulation.We model the acceleration and particle production following the evolution of the comoving shell, moving out from 1 × 10 11 cm, as sketched in the left panel.Bottom: Evolution of Γ rad , ρ ′ and σ for θ = ϕ = 0 (thus, along the dotted line in the upper panel) for 5.0, 5.5, 6.0, 6.5, and 7.0 s.In the left plot, we indicate the average Lorentz factor of the last snapshot ⟨Γ rad ⟩7s with a dashed purple line.The shaded band marks the jet region considered in our multi-messenger emission modeling, which in the two right plots corresponds to the thick, non-transparent lines.In the middle bottom plot of ρ ′ , we show the power-law extrapolation of the density as a dashed black line, while in the right plot of σ we indicate the average magnetization at 7 s within the jet region as dashed purple line (see main text for details).
with homologous ejecta with mass M ej = 0.05 M ⊙ (inspired by GW 170817, see e.g.Kasen et al. (2017)), and azimuthally-symmetric baryon density profile: here θ is the polar angle and ρ 0 is determined by requiring that M ej = 0.05 M ⊙ .
We employ an ideal equation of state with a relativistic adiabatic index of 4/3, within the full simulation setup provided in Gottlieb et al. (2022).We follow the jet evolution from its launch by the post-merger BH for 7 s until its head reaches ∼ 2 × 10 11 cm.
The upper panels of Fig. 1 display isocontours of the radial Lorentz factor Γ rad , the comoving mass density ρ ′ , and the magnetization σ, from left to right respectively, extracted at 7 s for representative purposes, in the xz-plane (thus ϕ = 0; ϕ being the azimuthal angle) for 7 × 10 10 cm < r < 2.1 × 10 11 cm, and viewing angle −0.4 rad < θ < 0.4 rad1 .From the profile of Γ rad (left panel), one can identify the relativistic jet between 1.4-2.0× 10 11 cm, roughly propagating along the z-axis, and a mildly relativistic cocoon surrounding the jet.The jet extension corresponds to an engine activity time of about 2 s, roughly consistent with the typical duration of sGRBs.It is important to note that the outflow is still optically thick at 7 s with optical depth τ ≳ 10 4 below 2.1×10 11 cm along the z-axis (see Fig. 5 in Appendix A).The high-density structures surrounding the jet stem from the cocoon (light pink regions in the middle panel of Fig. 1).In contrast, the relativistic jet has a lower density.The magnetization, σ = B ′2 /4πρ ′ c 2 (where B ′ is the comoving magnetic field), in the right panel of Fig. 1 follows the inverse pattern: the largest σ is visible in regions of low ρ ′ .
In order to investigate the efficiency of particle acceleration up to the photosphere (located at r ≳ 10 12 cm) and related multi-messenger emission, we focus on the region with 10 11 cm ≲ r ≲ 10 12 cm.Within this radial range, we follow the comoving evolution of a shell representing the relativistic jet, that e.g. at 7 s is located between 1.48 × 10 11 cm ≲ r ≲ 1.8 × 10 11 cm (see the top panel of Fig. 1), where the radial Lorentz factor Γ rad ≳ 10.At the viewing direction θ = ϕ = 0, the jet extension in radius is maximal, while the internal energy reaches its highest values.In this sense, it may dominate the overall observed emission and be representative of the jet region.The jet region with r ≳ 1.8 × 10 11 cm is shaded in the top panel of Fig. 1, as the simulation here is subject to numerical artifacts that plague the evolution of the magnetization.Consequently, this region is ignored in our post-processing of the jet simulation to compute the multi-messenger emission.Since the full evolution up to the photosphere is not captured by the GR-MHD simulation, we rely on extrapolations obtained as follows.
At each radius, the shell is characterized by its Lorentz factor Γ, magnetization σ, internal energy density u ′ , and rest mass density ρ ′ .In order to extrapolate the behavior of these quantities, we select the snapshots taken at [5.0, 5.5, 6.0, 6.5, 7.0] s (see the bottom panels of Fig. 1 and Appendix A); note that we choose to begin the computation of particle emission at 5 s when the jet characteristic quantities almost reached asymptotic values, approaching the end of acceleration and adiabatic expansion, the photon spectral energy distribution being thermal earlier.Both the radial Lorentz factor as well as the magnetization change marginally between the last snapshots; thus, we assume constant Γ ≡ ⟨Γ rad ⟩ 7s = 12.5 and σ ≡ ⟨σ⟩ 7s = 0.18 (here, the label 7 s refers to the considered jet simulation snapshot).For u ′ and ρ ′ , which due to expansion of the system clearly decay with time, we instead extrapolate the time evolution from the jet simulation.In order to obtain representative values of these quantities up to the photosphere, we first perform a fit for each quantity of interest, relying on the shaded cells displayed in the bottom left panel of Fig. 1 and then calculate the average of each quantity.This procedure yields ρ ′ (r) = 2.3 × 10 28 (r/cm) −2.96 g/cm 3 and u ′ (r) = 5.6 × 10 55 (r/cm) −3.56 erg/cm 3 .Note that this deviates from the behavior expected for the simple, self-similar expansion of the shell (that would yield ρ ′ ∝ r −2 ).This hints that the acceleration has not fully ceased (in fact, we find that the average specific enthalpy ⟨h ′ ⟩ = ⟨1 + 4p ′ /(ρ ′ c 2 )⟩ ≳ 2 along θ = ϕ = 0 in the considered region).Future GR-MHD simulations up to even larger radii are imperative to tackle this problem fully.These fits are somewhat dependent on θ, yet for all viewing angles we find that the decay of ρ ′ is steeper than r −2 and the difference in the power-law indexes for u ′ and ρ ′ is between 0.2 and 0.4, with u ′ always decaying faster than ρ ′ .The latter potentially affects the formation of the Wien spectrum stemming from the Planck spectrum and the evolution of the thermal photon peak (Vurm et al. 2013;Bégué & Pe'er 2015).

PARTICLE ACCELERATION SITES
We now identify the regions of energy dissipation in the jet and subsequent particle acceleration, focusing on the jet region at r ≳ 10 11 cm and assuming the representative shell evolution introduced above.Analyzing our benchmark simulation, we identify the following mechanisms of energy dissipation: 1. Magnetic reconnection.The low magnetization σ = 0.18 in principle inhibits magnetic reconnection (Blandford & Znajek 1977;Drenkhahn & Spruit 2002;Drenkhahn 2002;Gill et al. 2020).However, the magnetization is locally enhanced up to σ = 7.5 in the jet region between 1-1.6×10 11 cm and for the narrow range of θ covering the relativistic jet, with slight variations in the exact positions/distributions of the peaks.Since it is too expensive to solve the transport equations (see § 4) while considering the small fluctuations in σ, we take into account these local maxima of σ by assuming |dσ/dr| = const.,where the magnetization builds up between 1-1.3 × 10 11 cm before decaying between 1.3-1.6 × 10 11 cm (see Fig. 2).Although this simplified profile may not be representative of the complex reality, we do not expect this choice to impact our results as long as the average |dσ/dr| and σ are similar in both cases.Following the GR-MHD simulation, we do not account for magnetic dissipation above 1.6 × 10 11 cm.The induced rate of energy dissipation during the decay of the magnetic field is 2. Internal shocks and collisionless sub-shocks.If the jet Lorentz factor varies on a timescale t v (length scale ct v ), internal shocks take place at the radius 2t v Γ 2 c (Rees & Mészáros 1994).We identify variability in the jet Lorentz factor of our benchmark simulation on timescales t v ∼ 50-500 ms.Internal shocks are thus expected at 4.7 × 10 11 cm ≲ r ≲ 4.7 × 10 12 cm.These shocks convert kinetic into internal energy, which causes a decrease in the Lorentz factor of the plasma.Hence, the dissipated energy can be estimated through the expected slow-down of the plasma.To compute the dissipated energy, we follow a procedure similar to the simple internal shock model (e.g.Barraud et al. 2005) based on the interaction between a fast (f) and a slow (s) plasma shell.Within the considered region, the fastest plasma moves with Γ f ≃ 16 at r f = 1.5 × 10 11 cm, whereas the slowest plasma has Γ s ≃ 9 at r s = 1.5 × 10 11 + t v c cm.For the maximal variability timescale of t v = 500 ms, this simple model yields ∆Γ ≃ 0.5.Although single shocks may dissipate energy at specific radii, the continuous velocity profile can be expected to yield a smooth dissipation of energy.We emulate this by assuming dΓ/dr = ∆Γ/∆r IS =const.We thus compute the comoving rate of energy dissipation via internal shocks as When internal shocks are radiation mediated (as expected in our case for the optically thick plasma below the photosphere, τ ≳ 1 in the region of 10 11 cm ≲ r ≲ 10 12 cm; see Appendix A), particle acceleration is inhibited (Levinson & Nakar 2020).However, collisionless sub-shocks that form within the radiation-mediated shock may enable particle acceleration if the plasma is mildly magnetized (0.1 ≲ σ ≲ 1) and holds (Beloborodov 2017).While, for simplicity, we do not follow numerically the evolution of the shock structure consistently with the spectra of the accelerated particles (Budnik et al. 2010; Beloborodov 2017), at 4.7 × 10 11 cm ≲ r ≲ 1 × 10 12 cm, we find that collisionless sub-shocks meet the criteria introduced above and are therefore an efficient mean of accelerating particles in our jet model; the energy released in sub-shocks is obtained through Eq. 3 (which is effectively an upper limit, for the sub-shocks embedded in the radiation mediated shocks).
3. MHD turbulence.Acceleration at MHD turbulences (possibly in connection with shocks) has been proposed as an alternative particle acceleration mechanism in GRBs (Lemoine 2021;Bresci et al. 2023).Given that it is not expected to operate efficiently in the collisional plasma below the photosphere, it is omitted in our modeling.
In synthesis, as sketched Fig. 2, in our benchmark simulation, the jet crosses three different regions as it propagates between 10 11 -10 12 cm: 1.The reconnection region in the range [1, 1.6] × 10 11 cm, where a spike in σ induces magnetic reconnection; 2. the expansion region, in the range [1.6, 4.7] × 10 11 cm, where no additional energy dissipation occurs; 3. the sub-shock region above 4.7 × 10 11 cm, where particles are accelerated at collisionless sub-shocks.In Fig. 2, we further schematically show how σ and Γ are affected in the dissipation regions.Note that the photosphere is located still at slightly larger radii, as can also be inferred from the Compton-Y -parameter that equals ∼ 10 at 10 12 cm (see Fig. 6 in AppendixA).

MULTI-MESSENGER EMISSION
We now introduce the transport equations of electrons, positrons, protons, neutrons, pions, muons, neutrinos, and photons, as well as their related particle distributions.We also present our findings on the temporal evolution of the particle spectral distributions.

Particle transport equations
In a purely thermal plasma, electrons and photons may interact through Compton scatterings, as well as through thermal bremsstrahlung, double Compton, and cyclotron emission or absorption.In the case of efficient particle acceleration, numerous leptonic and hadronic emission processes may shape the multi-wavelength spectra: First, all charged particles in a magnetic field emit synchrotron radiation and participate in inverse Compton scatterings.The plasma may be enriched by energetic electron-positron pairs by means of γγ-annihilation or photo-hadronic pair production, where the secondary pairs may contribute significantly to the photon spectra through the synchrotron and inverse Compton processes.Finally, hadronic processes such as proton-proton interactions and photo-pion production yield energetic pions that decay into photons, muons, electrons, and neutrinos.Here, the secondary neutrinos serve as unique probes of hadronic interactions.
Using the inputs at each r from our jet model, as illustrated in § 2, we solve the time-dependent coupled partial differential equations of (non-thermal) electrons, positrons, protons, neutrons, pions, muons, neutrinos and photons: with δ/δ x ≡ δ x , n being the particle number (or occupation number, for momentum space), and x representing the dimensionless particle energy/momentum.In our treatment for photons, x represents the momentum and the pre-factor A(x) = 4πx 2 .For all other (non-thermal) particle species, x represents the energy and the pre- The energy distributions of non-thermal electrons, positrons, protons, neutrons, pions, muons and neutrinos are evolved relying on a modified version of the AM 3 code (Gao et al. 2017).In addition, we extend the original code that includes synchrotron, inverse Compton, photo-pion production, Bethe-Heitler pair production, γγ-annihilation, and adiabatic cooling; our new version of AM 3 also accounts for proton-proton interactions (see Appendix B).Secondary particles are treated selfconsistently, undergoing the same interactions as primary particles.
The photon evolution is instead computed through a newly developed code.The latter treats the evolution in momentum space, accounting for Comptonization by a thermal electron population through the Kompaneets equation (Kompaneets 1957).We extend the Kompaneets equation to include cooling, absorption/emission terms due to thermal emission and absorption processes (bremsstrahlung, cyclotron, and double Compton) as well as the effects of adiabatic expansion.Non-thermal processes are included by offering an interface to AM 3 (γγ-annihilation, synchrotron, inverse Compton, and pion decay).The initial spectrum at r ≃ 10 11 cm is a Planck black-body spectrum.For the thermal protons, neutrons, and electrons, no separate treatment is necessary.For these instead, at each time step, the tem-perature and number density are updated, relying on the evolution of the internal and rest mass density of the plasma shell.The assumption of Maxwellian distributions for the thermal populations is justified since the thermal relaxation time due to Coulomb collisions is shorter than the dynamical time, for the radial range under consideration.For details on the considered processes and the numerical treatment, we refer the interested reader to Appendix B. The interface between the GR-MHD simulation and the radiative calculations is thus implemented on the one hand through the thermal populations (following the evolution of ρ and u), and on the other hand through the injection terms for the non-thermal electrons and protons that are a result of the inferred particle acceleration regions and processes.

Accelerated particle distributions
We describe the accelerated spectrum for the particle species i (i.e., protons or electrons; see injection term in Eq. 5) as a power-law distribution with exponential cutoffs both at the minimum Lorentz factor γ ′ min and the maximal Lorentz factor γ ′ max : The factor N i is computed by normalizing with respect to the energy dissipation rate du ′ /dt ′ .The splitting of du ′ /dt ′ between the different particle populations is parameterized through ϵ th (that defines the fraction of energy transferred to thermal protons, neutrons, and electrons) and ϵ i (the fraction of non-thermal energy transferred to particles of species i).Then, N i can be calculated through Similarly, the thermal plasma is heated by ϵ th du ′ /dt ′ .The latter increases the temperature of the thermal proton, neutron and electron populations.The minimum and maximum Lorentz factors, as well as the power-law slopes p i , ϵ th , ϵ i depend on the acceleration process at work.As these parameters specify the properties of the primary injected distributions, they are not altered by the radiative processes.In contrast, the cooled comoving distributions may differ significantly; for example, secondary lepton pairs can dominate the electron distribution.
We summarize our assumptions on the factors entering the accelerated particle distributions for magnetic reconnection and collisionless sub-shocks in Table 1.Here, we assume that magnetic reconnection operates similar to electron-ion reconnection in the optically thin regime, as explored e.g. in Sironi et al. (2015); Werner et al. (2018); Zhang et al. (2023).However, even in a collisionless plasma, magnetic reconnection may be slowed down due to the injection of electron-positron pairs (Hakobyan et al. 2019)-further work is required regarding reconnection in collisional plasmas (Uzdensky 2011).
The maximum particle energy attainable for acceleration at sub-shocks can be computed by balancing the acceleration characteristic timescale with the energy losses (and their corresponding loss rates t ′ loss ).For a particle of Lorentz factor γ i and mass m i , the acceleration rate is (Globus et al. 2015) with q i being the particle charge.As for magnetic reconnection, the synchrotron burnoff may be exceeded for acceleration along the electric field, where the maximum energy scales with the magnetization (see e.g.Cerutti et al. 2013;Kagan et al. 2016Kagan et al. , 2018)), thus γ ′ e,max = 5σm p /m e and γ ′ p,max = (γ ′ e,max −1)m e /m p +1.
Notes: 1.We neglect non-thermal proton acceleration in magnetic reconnection for σ ≲ 10.In this case, ϵpdu ′ /dt ′ is instead a heating term for the thermal populations.11.9 log10(r/cm) Figure 3. Evolution of the comoving spectral energy distribution of photons with radius, split in three regions (see also Fig. 2).In the magnetic reconnection region (first panel from the left), the local enhancement of the magnetization enables magnetic reconnection and the subsequent acceleration of electrons.In the expansion region (second panel), no energy dissipation takes place and therefore the photon distribution peak shifts to lower energies, with overall lower photon density as a result of the plasma expansion.Finally, in the sub-shock region (third panel), both protons and electrons are accelerated at collisionless sub-shocks, leading to the appearance of a non-thermal high-energy tail and a softening of the spectrum below the distribution peak.In order to highlight the impact of hadronic processes, the photon distribution obtained without non-thermal protons is also shown in the sub-shock region (fourth panel).To highlight the evolution of the photon distribution, a power-law fit at the final snapshot of each region is provided.
Figures 3 and 4 show our findings on the time evolution of the photon and neutrino spectral energy distributions, respectively, between 10 11 cm and 10 12 cm obtained by solving Eq. 5 and relying on the evolution of the shell introduced above.In the following, we outline the main features of the particle distributions considering the three jet regions introduced before (see Fig. 2): 1. Magnetic reconnection region.During the rampup of the magnetization σ at 1.0 × 10 11 cm < r < 1.3 × 10 11 cm, the photon distribution evolves as a thermal one; adiabatic expansion and the decreasing thermal electron temperature are responsible for the shift of the distribution peak at lower energies and a decrease in number density.The slope of such a spectrum is As magnetic reconnection becomes active for r ≳ 1.3 × 10 11 cm, two effects manifest: 1.A nonthermal photon population is injected.The latter appears as a non-thermal tail at E ′ γ ≳ 10 2 keV in the initial power-law distribution.At lower energies, efficient Comptonization of the synchrotron seed photons induces a softening of the spectrum between E ′ γ ∼ 10 −1 -1 keV.At the lowest energies, the plasma is still in thermal equilibrium, its thermal shape induces a "bump"-feature at the turnover energy.2. Due to the low σ ≲ 10, protons are not accelerated above thermal energy.Sub-sequently, roughly half of the dissipated energy heats the thermal population.Since photons and electrons are still coupled, the peak of the photon thermal spectrum is consequently shifted to higher energies.As the energy dissipation ceases and the electron acceleration stops, the non-thermal signatures directly disappear and a narrow Wien spectrum approximately scaling as E ′2.43 γ is evident (light orange line in the leftmost panel of Fig. 3).Due to proton acceleration being inefficient, no neutrino production is expected.

Expansion region. As the jet expands without
further energy dissipation, the plasma cools and dilutes.Hence, the peak of the photon distribution moves to lower energies, and lower photon densities are achieved.The redistribution of photons abundantly present due to the previous dissipation phase yields a spectrum scaling approximately as E ′1.12 γ (in the region where the spectrum is not in thermal equilibrium for E ′ γ ≳ 10 −2 keVsee the second panel to the left of Fig. 3).As the plasma moves outwards, the turnover energy at which photons are still in equilibrium shifts to lower energies, broadening the the soft part of the spectrum.
3. Sub-shock region.Protons and electrons accelerated at collisionless sub-shocks lead to the formation of a non-thermal high-energy tail, as well as a strong softening of the spectrum below the distribution peak (with a spectrum scaling approximately as E ′0.7 γ , as displayed in the third panel of Fig. 3; the thermal distribution peak is approximately preserved despite the dissipation processes. We highlight the important role of hadronic processes in shaping the photon distribution.In fact, if only electron acceleration is accounted for (rightmost panel of Fig. 3), the photon distribution basically maintains its former shape.This strong effect of hadronic dissipation is linked to the opacity of the plasma, where a secondary cascade of lepton pairs is injected due to γγ-annihilation of the very-high-energy pion decay photons as well as through Bethe-Heitler pair production.In the sub-shock region, proton-proton interactions operate efficiently, being the dominant energy loss process for low-energy protons.
The evolution of the neutrino energy distribution produced by the subsequent decay of energetic pions is depicted in Fig. 4. The relatively low maximum energy of neutrinos is attributed to efficient cooling of protons due to photo-pion and Bethe-Heitler processes.We also find that the neutrino energy distribution is mainly shaped by proton-proton interactions, the proton-photon channel being sub-leading and leading to the production of neutrinos with slightly higher energies on average.Note that we neglect the contribution of kaons in shaping the neutrino energy distribution.Kaon decay is expected to affect the neutrino spectrum at high energies (Ando & Beacom 2005;Asano & Nagataki 2006).Below the photosphere and for a collapsar jet with σ 0 = 200, Guarini et al. (2023) showed that kaon decay leads to a bump above O(10 3 ) GeV in the neutrino distribution due to the large magnetic field.However, in our case, the low maximum energy of neutrinos and the dominance of proton-proton interactions suggest that the kaon contribution would lead to negligible modifications.Note that our simulation includes no escape, but an adiabatic expansion term for neutrinos.In this sense, the plotted energy distributions are representative of time-integrated ones, although corrected for the increasing comoving volume.Let us highlight two aspects characterizing the temporal evolution of the neutrino distribution.After the first (green) line, the later neutrino spectra are almost identical with their corresponding lines lying on top of each other.This is indicative of an almost steady-state proton population.Second, the peak of the neutrino distribution lies around 10 5 keV.Note that, if represented as energy flux, this peak is shifted to an energy of 1 GeV, which for a typical redshift of z = 2 would yield an observed peak energy of E ν = 4.1 GeV.This rather low peak energy may be understood by the low maximum proton energies, where photo-pion reactions potentially producing neutrinos of higher energies are inhibited, similar to the findings of Guarini et al. (2023).
Overall, the period of magnetic energy dissipation at small radii tends to increase the number of available photons and yields a softened spectrum in the Wien zone.The cascades initiated by hadronic processes in the sub-shock region soften these spectra even further, in addition to a high-energy non-thermal tail.
The magnitude of the reported signatures naturally scales with the energy dissipation rate.Further, we focus on specific regions where magnetic reconnection and internal shocks are expected.However, magnetic dissipation and sub-shocks due to variability on shorter timescales, may also affect the spectrum in the expansion region, but are not addressed in this work.  .Evolution of the comoving spectral energy distribution of neutrinos (summed over all flavors).Neutrino production is dominated by proton-proton interactions.We only focus on the third region (sub-shock region), since neutrino production is inefficient in the region of magnetic reconnection.Note that no escape term is included in the differential equation for the neutrinos; and the shown spectra rather represent time-integrated results (corrected for adiabatic expansion).The same color scale as in Fig. 3 is adopted.

DISCUSSION
The origin of the spectral shape of the gamma-ray spectrum observed for (s)GRBs remains unclear (Nakar 2007;Baiotti & Rezzolla 2017).Optically thin synchrotron spectra tend to yield spectra that are too soft/wide, while (re-processed) thermal spectra are rather too hard or narrow.On the other hand, the photon spectral energy distribution, shown in Fig. 3 at r ∼ 10 12 cm, scales as dn ′ /dE ′ γ ∝ E ′−0.3 γ below the peak.This behavior is compatible with the power-law slopes observed for a fraction of sGRBs, although in most cases the spectral index at low energies is still lower (Burgess et al. 2019b;Poolakkil et al. 2021).This could be resolved through enhanced energy dissipation, such as stronger shocks due to larger differences in Lorentz factors, considering more shocks or an extended magnetic reconnection region.Further, softer spectra can be obtained if the evolution of u and ρ is such that the energy released by sub-shocks in comparison to the internal energy density is larger than in the current scenario.
It is important to stress that, if hadronic processes are not accounted for, the photon distribution is significantly harder with dn ′ /dE ′ γ ∝ E ′0.25 γ (see Fig. 3), which is potentially difficult to reconcile with sGRB gammaray spectra.Note that similar to the neutrino spectra, photon spectra are mostly impacted by pp-interactions instead of photo-pion production.At high energies, we find that synchrotron emission of secondary non-thermal electrons results in a high-energy tail of the spectrum scaling as dn , which is steeper than the average observed slope-although for sGRBs a low highenergy spectral index has also been reported (Acuner & Ryde 2018;Poolakkil et al. 2021).Again, this signature is only visible if hadronic processes are included.Being dominated by secondary pairs, the high-energy slope differs from the one expected for an electron population with p e ∼ 2.1-2.2.Finally, the comoving peak energy of photons at O(1 keV), obtained in Fig. 3 and a relic of the peak of the thermal distribution at 10 11 cm, would require a high(er) jet Lorentz factor in order to be compatible with typical sGRB peak energies of O(10 2 -10 3 keV) (Burgess et al. 2019b;Poolakkil et al. 2021).
Here we point out that noticeable relics of the thermal spectrum are still present around the peak in our final photon distribution; stronger energy dissipation, potentially also taking place beyond the photosphere, may be required to wash these features out (see also Vurm et al. 2011;Vurm & Beloborodov 2016).
Dominated by the products of proton-proton interactions, the neutrino distribution [E ′ ν dn/dE ν ] obtained in our simulation and shown in Fig. 4, peaks at E ′ ν ≃ 0.02 GeV.The corresponding energy flux [E ′2 ν dn/dE ν ] peaks at E ′ ν ≃ 1 GeV, equivalent to an observed energy of E ν = 4.1 GeV (for a typical redshift of z = 2) as a consequence of the low Lorentz factor Γ = 12.5.This finding is compatible with the results presented for collapsars in Guarini et al. (2023) and with the ones of Gottlieb & Globus (2021) obtained for sGRBs, where the production of neutrinos at collisionless sub-shocks at r ∼ 10 9 cm was considered and a maximal neutrino energy of 10 3 GeV obtained.Relying on the findings of Guarini et al. (2023); Gottlieb & Globus (2021), we expect that the low-energy tail of the neutrino distribution might be boosted by less efficient neutrino production occurring at r ≲ 10 11 cm.However, our results are in contrast with the ones previously reported in e.g.Razzaque et al. (2004); Ando & Beacom (2005); Tamborra & Ando (2016) and Murase & Ioka (2013) for lowluminosity GRBs, forecasting that neutrinos with TeV-PeV energy could be produced in choked jets harbored in collapsars because of collisionless internal shocks or collimation shocks (the latter are unlikely to be efficient in our case because of the large optical depth of the jet).The neutrino spectra from the jet could be complemented by neutrino emission from other components, see e.g.(Decoene et al. 2020;Kimura et al. 2018).However, we expect these to be sub-dominant along the line of sight of the jet.
The detection of this low-energy neutrino flux may be challenged by the atmospheric background in this energy range, however timing information could favor its detection in Hyper-Kamiokande and IceCube Deep Core (Guarini et al. 2023).
Future work should be addressed to simulate the longterm evolution of sGRB jets and the related neutrino and radiation transport, reaching up to optical depth τ ≃ O(0.1).In this context, the still ongoing acceleration of the plasma should be accounted for.Note, however, that as long as the magnetization in the subshock region supports collisionless subshocks, we do not expect a major impact on our main conclusions.The evolution of the shock structure should be also modeled consistently with the spectra of the accelerated particles.Moreover, growing evidence points towards a major impact of the magnetic field profile, its amplification mechanism, as well as a non-trivial coupling between energy deposition by neutrino annihilation and the magnetic field in the jet launching and its dynamics (Christie et al. 2019;Ciolfi 2020;Mösta et al. 2020;Gottlieb et al. 2022Gottlieb et al. , 2023a,b,c;,b,c;Hayashi et al. 2022Hayashi et al. , 2023;;Combi & Siegel 2023;Kiuchi et al. 2023a,b); all these factors may further affect the sites of energy dissipation as well as the associated multi-messenger emission.
Observations of sGRBs suggest a large variety in the temporal evolution of their signals.This may be due to a different distribution of the energy dissipation regions, as well as a different jet magnetization, and magnetic field configurations.We intend to investigate the impact of the jet properties on the multi-messenger observables in a companion paper (Rudolph et al. 2023b).Along the same vein, one may also consider different shells propagating outwards and representing a non-homogeneous jet region, in addition to taking into account the dependence on the viewing angle/angular jet structure.

CONCLUSIONS
In this Letter, we explore the multi-messenger signatures of sub-photospheric energy dissipation in sGRBs.We carry out a state-of-the art GR-MHD simulation of a binary neutron star merger aftermath, where we follow the propagation of a relativistic jet with initial magnetization σ 0 = 150 in massive ejecta for 7 s.For the first time, we simulate the coupled evolution of non-thermal and thermal particles in a plasma shell representative of the jet region that propagates between 10 11 and 10 12 cm (the photosphere being located at a slightly larger radius), including all relevant thermal and non-thermal processes.The evolution of the characteristic parameters of the plasma shell (such as the internal energy, rest mass density, Lorentz factor, and magnetization) have been extrapolated from our GR-MHD simulation up to 10 12 cm, the Lorentz factor and magnetization having already reached their asymptotic values at 10 11 cm.
In our model, we identify a region of magnetic energy dissipation at [1, 1.6] × 10 11 cm, where σ is locally enhanced by 7.5 enabling magnetic reconnection.While acceleration of protons above non-thermal energies is inhibited (hence, neutrinos cannot be produced), radiation of non-thermal leptons introduces a weak nonthermal high-energy tail and a low-energy feature in the spectral distribution of photons.In the region between [1.6, 4.7]×10 11 cm, the plasma expansion cools the thermal populations, and the photons are redistributed at lower energies, with a consequent softening of the photon distribution below the thermal peak.
We identify variability in the Lorentz factor, leading to internal shocks between [4.7 × 10 11 , 10 12 ] cm.Although particle acceleration is inefficient at internal shocks in a radiation dominated environment, protons and electrons can be efficiently accelerated at collisionless subshocks because of the mild magnetization in the region [4.7 × 10 11 , 10 12 ] cm considered in our model.The resulting non-thermal radiation from the primary accelerated particles and the induced secondary cascade significantly alters the spectrum with low-and high-energy asymptotic distributions scaling as n , while the distribution peak derives from the thermal one at 10 11 cm.The main modifications of the spectral energy distribution of photons in this region are due to hadronic processes, highlighting the need to include such processes and the resulting particle cascades when modeling sub-photospheric GRB spectral energy distributions.The obtained asymptotic slopes are compatible with the one observed in a fraction of sGRBs.Further energy dissipation not accounted for in the present scenario (both below and above the photosphere) may yield photon distributions closer to the typical observed spectra.The spectral peak of the photon energy distribution is softened during the expansion to ∼ 1-10 keV, which may require a larger Lorentz factor in order to reach the typical observed peak energy of sGRBs.Neutrinos are also produced in the subshock region, where efficient proton-proton interactions, Bethe-Heitler induced cooling as well as sub-dominant photo-pion production lead to a distribution peaking in the sub-GeV to GeV regime.
We conclude that photon spectra compatible with sGRB observations may be produced from radiation of protons and leptons accelerated below or close to the photosphere.Initiating an effective cascade, hadronic processes are crucial for transitioning from a Wien distribution to a true non-thermal one.Signatures of hadronic processes may be naturally carried by neutrinos as well.Future work extending our analysis to and above the photosphere, as well as including contributions from different viewing angles would be desirable to further confirm our findings.

ACKNOWLEDGMENTS
We are grateful to Eli Waxman for insightful discussions and to the anonymous referee for constructive feedback.This project has received funding from the Villum Foundation (Project No.  .Evolution of the jet characteristic quantities within the shell considered for particle acceleration, as it propagates outwards from 10 11 cm to 10 12 cm.The upper panels show the dimensionless electron thermal electron temperature (calculated from the electron temperature Te as θe = kBTe/mec 2 ) on the left and the Compton-y parameter with a [black, red] dotted horizontal line indicating where [yCompton = 10, yCompton = 1] on the right.The lower panels show the comoving magnetic field B ′ (left) as well as the energy density and rest mass density of the thermal population (composed of protons, neutrons and electrons, on the right).In all plots we shade the magnetic reconnection region in grey, and the sub-shock region in light blue.
and number densities of the thermal populations.We have verified that, for the parameters and jet radial range of interest, neutrons, protons, and electrons remain in equilibrium (i.e., the Coulomb relaxation time is much smaller than the dynamical time of the system).The density of each particle species is proportional to ρ ′ (r) of the plasma shell (with the electron-to-baryon ratio Y e ≃ 0.5, see e.g.Just et al. (2022)).The evolution of the thermal temperature θ e can be obtained by equating here du ′ γ /dt ′ exp represents the photon cooling term due to the expansion of the plasma.Second, we calculate the evolution of non-thermal leptons, muons, pions, protons and neutrinos from Eq. 5 relying on a modified version of the AM 3 code (Gao et al. 2017).The latter has been altered with respect to its original version in the following aspects: 4. Quantum synchrotron radiation.Owing to the large magnetic field [O(10 8 G)], quantum synchrotron radiation introduces a cut-off in the electron synchrotron emission that may be expressed as a modified synchrotron emission kernel (Brainerd 1987).For a more efficient computational implementation, we follow Imamura et al. (1985) and model this effect as a simple upper limit on the emitted photon energy (or a lower limit on the energy of the radiating electron).We have however verified that the results are in good agreement with respect to the ones obtained through the full treatment proposed in Brainerd (1987).Finally, photons are evolved with a separate, novel code that contains both emission/absorption from non-thermal and thermal particles, accounting for Comptonization due to the thermal electron population.The emission and absorption terms due to non-thermal processes are calculated with AM 3 and passed to the photon simulation code.In addition, we account for emission and absorption due to thermal bremsstrahlung, cyclotron and double Compton emission (Rybicki & Lightman 1986;Vurm et al. 2013;Bégué & Pe'er 2015).As the electrons are in thermal equilibrium, Kirchhoff's law can be applied.As a consequence, the emission ϵ(x) and absorption rate α(x) are linked through the following relation Here n is the photon occupation number, n e the total number density of thermal electrons and σ T is the Thomson cross-section.Adiabatic expansion leads to the same cooling and escape terms introduced above.The evolution of the photon population in energy and time is calculated relying on the Chang and Cooper scheme (Chang & Cooper 1970;Park & Petrosian 1996), where we choose an equally spaced grid in ln(x).Note that for all particle species we omit escape terms, since we are still in the optically thick regime even at the largest radius of 10 12 cm.

Figure 2 .
Figure2.Sketch of the particle acceleration sites in our sGRB jet simulation (not to scale).The central BH (dark red) launches a jet (cyan) which is surrounded by a mildlyrelativistic cocoon (dark purple).The jet propagates outwards, breaking out from the merger ejecta (pink); note that the jet is not present anymore at smaller radii since it was active for ≃ 2 s-the sub-relativistic cocoon expands in the region previously occupied by the jet.At 10 11 cm, it has reached a magnetization σ = 0.18 and Lorentz factor Γ = 12.5; we neglect the jet region below this radius, since we are interested in processes close to the photosphere potentially shaping the observed spectrum.A local increase in σ occurs between 1-1.3 × 10 11 cm, while the decay of σ enables magnetic reconnection (equivalent to 1 s duration).The induced change in magnetization is ∆σ = 7.5.Variability in the Lorentz factor of the jet leads to collisionless sub-shocks between 4.7 × 10 11 cm and 1 × 10 12 cm; as a consequence, a decrease in the Lorentz factor of ∆Γ = 0.5 occurs.The photosphere at which the jet becomes optically thin is represented through the solid ochre line.

Figure 4
Figure4.Evolution of the comoving spectral energy distribution of neutrinos (summed over all flavors).Neutrino production is dominated by proton-proton interactions.We only focus on the third region (sub-shock region), since neutrino production is inefficient in the region of magnetic reconnection.Note that no escape term is included in the differential equation for the neutrinos; and the shown spectra rather represent time-integrated results (corrected for adiabatic expansion).The same color scale as in Fig.3is adopted.
37358), the Carlsberg Foundation (CF18-0183), and the Deutsche Forschungsgemeinschaft through Sonderforschungsbereich SFB 1258 "Neutrinos and Dark Matter in Astro-and Particle Physics" (NDM).IT also thanks the Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy for partial support during the final stages of this work.OG is supported by Flatiron Research and CIERA Fellowships.OG also acknowledges support by Fermi Cycle 14 Guest Investigator program 80NSSC22K0031, and NSF grant AST-2107839.An award of computer time was provided by the ASCR Leadership Computing Challenge (ALCC), Innovative and Novel Computational Impact on Theory and Experiment (INCITE), and OLCF Di-

Figure 6
Figure6.Evolution of the jet characteristic quantities within the shell considered for particle acceleration, as it propagates outwards from 10 11 cm to 10 12 cm.The upper panels show the dimensionless electron thermal electron temperature (calculated from the electron temperature Te as θe = kBTe/mec 2 ) on the left and the Compton-y parameter with a [black, red] dotted horizontal line indicating where [yCompton = 10, yCompton = 1] on the right.The lower panels show the comoving magnetic field B ′ (left) as well as the energy density and rest mass density of the thermal population (composed of protons, neutrons and electrons, on the right).In all plots we shade the magnetic reconnection region in grey, and the sub-shock region in light blue.

5.
Cooling of intermediate particles.Relying on the AM 3 code version used inRudolph et al. (2023a), we include adiabatic, synchrotron, and inverse Compton cooling of intermediate pions and muons.Note that we neglect any contribution from kaons, following the procedure outlined inHümmer et al. (2010).
ϵ(x) = α(x) e x/θe − 1 , (B3)where x represents the momentum of the photon and θ e = k B T e /m e c 2 is the dimensionless electron temperature.Comptonization by a Maxwellian distribution of electrons is included through the Kompaneets equation:

Table 1 .
Summary of the quantities entering the accelerated particle distributions.