Multi-messenger model for the prompt emission from GRB 221009A

We present a multi-messenger model for the prompt emission from GRB 221009A within the internal shock scenario. We consider the time-dependent evolution of the outflow with its impact on the observed light curve from multiple collisions, and the self-consistent generation of the electromagnetic spectrum in synchrotron and inverse Compton-dominated scenarios. Our leptohadronic model includes UHE protons potentially accelerated in the outflow, and their feedback on spectral energy distribution and on the neutrino emission. We find that we can roughly reproduce the observed light curves with an engine with varying ejection velocity of ultra-relativistic material, which has an intermediate quiescent period of about 200 seconds and a variability timescale of $\sim1$~s. We consider baryonic loadings of 3 and 30 that are compatible with the hypothesis that the highest-energetic LHAASO photons might come from UHECR interactions with the extragalactic background light, and the paradigm that energetic GRBs may power the UHECR flux. For these values and the high dissipation radii considered we find consistency with the non-observation of neutrinos and no significant signatures on the electromagnetic spectrum. Inverse Compton-dominated scenarios from the prompt emission are demonstrated to lead to about an order of magnitude higher fluxes in the HE-range; this enhancement is testable by its spectral impact in the Fermi-GBM and LAT ranges.


INTRODUCTION
Gamma-ray bursts (GRBs) are extremely energetic explosions involving the collapse of a massive star or the merger of two compact objects. As their name suggests, GRBs release most of their electromagnetic output in γ rays within a short period of time (ranging from tens of milliseconds to hundreds of seconds). This brief and variable emission, known as the prompt GRB phase, is followed by the afterglow, a long-lasting multi-wavelength emission (for a review see, e.g. Kumar & Zhang 2015). The prompt emission is thought to be produced in a relativistic collimated plasma out-neutrinos (e.g. Waxman 1995;Vietri 1995). While the contribution of typical luminosity GRBs during their prompt phase to the diffuse neutrino flux measured by IceCube has been constrained to 1% of the diffuse astrophysical neutrino flux (Abbasi et al. 2022), the hypothesis that GRBs are UHECR sources cannot be ruled out yet (for a recent review see Kimura 2022). In fact, GRBs populating the high-end of the isotropic γ-ray energy distribution (E γ,iso > 10 54 erg) may provide the necessary energy output per event for powering UHE-CRs, requiring only a moderate baryonic loading (defined as the energy injected into non-thermal protons versus electrons) and without violating existing neutrino limits (Rudolph et al. 2022). In addition, energetic bursts may be detected as a single source (e.g. Gao et al. 2013).
On October 9th 2022, a very bright GRB was observed at redshift z = 0.151 (de Ugarte Postigo et al. 2022). The burst triggered the Gamma-Ray Burst Monitor (GBM) on board Fermi at 2022-10-09 13:16:59.000 UT (Veres et al. 2022), about an hour before the detection of a hard X-ray transient by the Burst Alert Telescope (BAT) of the Neil Gehrels Swift satellite (S. Lesage et al. 2022;Dichiara et al. 2022). The prompt phase of GRB 221009A consisted of a precursor (at about 10 s), followed by an extremely bright emission period about 200 s post GBM trigger. Overall, the prompt emission period lasted roughly 327 s and was composed of several peaks (S. Lesage et al. 2022). The preliminary Konus-Wind light curve showed several peaks of roughly 40 s duration (D. Frederiks et al. 2022); the initial peaks were separated from the latetime peak by a quiescent period of about 220 s. Some short-timescale variability on the order of seconds might be also visible in the 0.4-100 MeV AGILE MCAL light curve (A. Ursi et al. 2022) and INTEGRAL SPI-ACS light curve 1 . The burst was also observed in high energies (HE, 0.1 GeV) by the Fermi Large Area Telescope (LAT), starting about 200 s after the GBM trigger (i.e. during the prompt phase) and extending up to ∼ 25 ks into the afterglow phase (R. Pillera et al. 2022). In addition, very high energy (VHE, > 100 GeV) photons were detected by LHAASO, but their association to the prompt phase is not clear given that these photons were observed within a period of up to 2000 s after trigger (Yong Huang et al. 2022). It has been speculated that the highest energy photons (up to 18 TeV) might come from UHECR interactions with the extra-galactic background light (EBL), since such energetic photons escaping the source would be otherwise attenuated by the EBL (Das & Razzaque 2022;Alves Batista 2022;Mirabal 2022). These scenarios require a significant amount of energy carried by UHE protons, which might also leave signatures in the electromagnetic spectrum.
The extreme brightness of this burst caused pile-up in almost all GRB detectors, namely Fermi-GBM and LAT, KONUS-Wind, and AGILE. For this reason, detailed spectral analysis was unavailable at the time of writing. Preliminary analysis of LAT data (100 MeV -1 GeV) for 200 − 800 s after the GBM trigger provided an estimate of the photon index (−1.87 ± 0.04) and the photon flux (6.2±0.4)×10 −3 ph cm −2 s −1 ), but thet time interval excluded from the analysis was recently extended to 300 s . Hence, the relevance of these results to the main GRB episode detected by GBM is not yet clear. Moreover, preliminary analysis of Konus data during the brightest phase of the event (i.e. ∼ 180 − 200 s after the Konus trigger) produced a rest-frame peak energy of E peak = 1150 keV and E γ,iso 3 × 10 54 erg. No associated muon-neutrino track was detected by IceCube in a time range of [-1 hour, +2 hours] from the initial GBM trigger, which resulted in an upper limit on the muon neutrino fluence of 3.9×10 −2 GeV cm −2 assuming an E −2 neutrino spectrum (The IceCube Collaboration 2022). The inferred γ-ray isotropic energy, the proximity of this event to Earth, and the lack of prompt neutrino detection make GRB 221009A a unique case for the study of multi-messenger signatures from GRBs.
In this Letter, we present a multi-messenger model for the prompt emission of GRB 221009A. Under the assumption that protons and electrons are accelerated at the fastest possible rate in internal shocks occurring at different radii within the jet, we compute the multimessenger emission from each collision while taking into account the varying physical conditions in the outflow and the UHECR feedback on both the photon and neutrino emissions. Our goal is to test the hypothesis of UHECR acceleration in the GRB jet by comparing the self-consistently computed broadband photon spectrum and the accompanying neutrino flux with available observational information. Since details on the prompt spectrum are not yet available, and the mentioned pileup effects likely introduce some degeneracy, our model also has some predictive power. In this work we indicatively use the following observables: (1) the peak energy E peak = 1060 keV (as fitted for the onset of the bright emission period by Konus), (2) the estimated fluence in the GBM band F GBM = 2.91 · 10 −2 erg cm −2 (1−1000 keV) and (3) the approximate light curve structure observed by Konusand INTEGRAL SPI-ACS. Peak energy and fluence are reproduced within ±5 %.

MODEL AND IMPLEMENTATION
Our model is fully described in an accompanying paper (Rudolph et al. 2022) to which we refer the interested readers for details. Here we summarize the key points of the model and present its main parameters in Tab. 1. In short, a relativistic outflow is discretized as shells ejected with different Lorentz factors from a central engine operating over a time t eng (as measured in the engine rest frame). Shells catch up with each other and merge at a radius R Coll from the central emitter, where energy is dissipated and particles are accelerated. Shells continue to propagate, merge and dissipate energy until their velocity distribution is such that no more collision occur. The dissipated energy in each inelastic collision, which can be obtained from energy and momentum conservation, is distributed into non-thermal electrons, protons and magnetic fields (for simplicity we assume that no energy remains in thermal plasma). We introduce the partition parameters f p/e (baryonic loading) and f B/e (magnetic loading) describing the ratio between proton/electron and magnetic/electron energy density, respectively (here only referring to accelerated, non-thermal particles). Assuming that the sub-MeV prompt spectrum is predominately produced by synchrotron radiation of accelerated electrons, we normalize the fireball kinetic energy to the total energy transferred to non-thermal electrons E tot e,NT that is needed to produce a given E γ,iso . Based on our simulations the energy dissipation efficiency 2 of the fireball is ε diss 0.04, for the specific Lorentz factor distribution assumed here. The required fireball kinetic energy (isotropic equivalent) is then indirectly related to the energy partition parameters and, in particular, increases with increasing baryonic loading as E kin,ini = ε −1 diss E tot e,NT (1+f p/e +f B/e ). In contrast to Rudolph et al. (2022), we normalize all models to a similar E γ,iso , leading to higher E kin,ini for cases with a low electron synchrotron efficiency.
To reproduce the structure of the observed light curve we use a varying Lorentz factor profile, see Fig. 1 (left panel), here without addressing the question which engine properties would lead to such a profile. The initial Lorentz factor distribution is obtained in two steps: First, we reproduce the broad structure of the observed light curve (a dip and two bright peaks, followed by a quiescent period and late-time emission) through sine waves with relative amplitudes matching the observed flux variations and durations of these engine activity intervals. The engine activity is correspondingly characterised by three time intervals: The activity time of the main emission period t main , the engine quiescent time t quiet and activity time for the late-time emission t late . Second, after inferring a short-time variability timescale of δt var 1s from observations we add modulations on this timescale, with amplitude drawn randomly from a normal distribution with standard deviation 0.08 · Γ . From 20 realisations of this random process, we select the initial configuration that best matches the observed light curve by eye.
We present two different scenarios for the initial Lorentz factor distribution -denoted as "R 16 " and "R 17 " -that produce collisions at an average radius R Coll ∼ 10 16 cm and R Coll ∼ 2 · 10 17 cm respectively, and probe higher and lower typical plasma densities. We verified that in both scenarios the bulk of energy dissipation occurs below the estimated deceleration radius for typical parameters of the circumburst medium. These scenarios are motivated by estimates for the optical thickness of the emitting plasma to γγ pair production ) and the requirement that the bulk of energy is dissipated below the approximate deceleration radius. In principle, high Γ factors are also expected from the empirical E γ,iso − Γ relationships (Ghirlanda et al. 2018). As an example, we show for "R 16 " the two-dimensional distributions of the non-thermal electron energy and the Lorentz factor of the emitting shell with collision radius in the middle and right panels of Fig. 1.
We self-consistently evaluate the emission from nonthermal electrons and protons, accelerated in shocks from the collision of shells, and injected with power-law distributions into the radiation zone, which is the hot plasma of the merged shell. The full-burst emission is then integrated over all collisions, taking into account the curvature of the emitting surface. The maximal electron and proton energies are limited by the dominating energy loss processes assuming efficient acceleration (operating at the Bohm limit). For each scenario, the minimum electron Lorentz factor γ e,min is set such the peak energy of 1060 keV is reproduced; The minimal proton energies are generally fixed to γ p,min = 10. We assume a power-law slope of primary protons of p p = 2.0 as in Rudolph et al. (2022). In our model the primary electron power-law slope would usually be determined by the high-energy photon index. As there was no reliable measurement for the high-energy slope at the time of writing, we adapt p e = 2.2 as suggested for mildly Power-law index of accelerated protons pp 2.0 Minimum Lorentz factor of accelerated electrons γ p,min 10 Notes. -All listed energies refer to the isotropic equivalent values. Energies and times are reported in the source (engine) frame. Averaged quantities are computed weighing with the shell mass ( Γini and Γfin ) or with the dissipated energy of a collision ( Γem , RColl , Ep,max ). The maximal proton energies are computed considering synchrotron and adiabatic losses and equally refer to the engine frame.
relativistic shocks (Crumley et al. 2019). We then numerically solve 3 the coupled system of integrodifferential equations describing the temporal evolution of particle distributions (electrons/positrons, photons, protons, neutrons, pions, muons, neutrinos) including the relevant processes, such as synchrotron emission and absorption, inverse Compton scattering (ICS), photo-pair and photo-pion production, γγ pair production, adiabatic cooling, and escape. Our approach fully captures the electromagnetic cascade induced by photo-hadronic and γγ pair production in each merged shell 4 . As discussed in detail in Rudolph et al. (2022), the electromagnetic spectrum will be, even in the lepto- We therefore discuss a scenarios synchrotron ("SYN")dominated and inverse Compton ("IC")-dominated one with different values of f B/e (see Tab. 1), characterized by different photon spectra and light curves. We also include the effects of EBL attenuation using the model of Dominguez et al. (2011), calculated with the opensource gammapy-package (Deil et al. 2018;Nigro et al. 2019). A crucial parameter in every multi-messenger model is the baryonic loading. Rudolph et al. (2022) showed that f p/e 3 is required for E γ,iso 3 · 10 54 erg, if such energetic GRBs ought to power the UHECRs. Much higher baryonic loadings in combination with low R Coll lead to spectral distortions of the photon spectrum (even in the GBM and LAT bands), and efficient neutrino production, which might be in tension with multi-messenger limits (for details, see Sec. 6 in Rudolph et al. 2022). A different argument comes from the observed VHE photons in LHAASO, if these were produced by interactions of UHECRs with the EBL (Das & Razzaque 2022;Alves Batista 2022;Mirabal 2022). Das & Razzaque (2022) derive an isotropic-equivalent E p,iso 3.9 · 10 54 erg for escaping cosmic rays in the range 0.1-100 EeV, which yields E p,iso 2·10 55 erg bolometrically corrected for our considered energy range, if all UHECRs free-streamingly escape. In our model, this energy budget is met for f p/e 3. For Alves Batista (2022), the given fraction into UHECRs is lower, but the bolometric correction factor is much higher, as only UHECRs at the highest energies were considered in that work. In what follows, we choose f p/e = 3 for "R 16 " and f p/e = 30 for "R 17 ". The former is compatible with the above estimates, while the latter is a more aggressive version (still compatible with neutrino bounds) that will challenge the GRB energetics in terms of the kinetic energy required. Purely leptonic results will not be shown, as we found no substantial modifications due to hadronic contributions in the observable photon spectra for the baryonic loadings considered.

RESULTS
Light curves. By construction of our model, the rough structure of the light curves shown in Fig. 2 reflects the initial Lorentz factor distribution in Fig. 1. Each pulse of the synthetic light curves is formed by shells colliding over a wide range of radii (as can be inferred from the colour coding in Fig. 1). Collisions within the same pulse may thus have different opacities, potentially introducing time-lags between different energy ranges.
In the left panel, we indicate the observed light curves of INTEGRAL SPI-ACS and Konus, shifted to match the beginning of our synthetic light curve that does not include the precursor emission. The general structure of the observed light curve is reproduced by our model, but our results should also be representative for light curves with a similar general structure. The relative intensity of the two bright pulses in the observed light curve may strongly be impacted by pile-up effects in the detector. Here, we assumed that the first bright peak is intrinsically much more energetic than the second one. In our model the relative height of the peak is controlled by the ratio of Lorentz factors of the colliding shells (e.g. a bright peak is produced when this ratio is large). Available γ-ray light curves at the time of writing indicate short-time variability on the second timescale, which was modeled by stochastic variations of Γ ini on a timescale δt var = 1.4 s (see Fig. 1). However, if refined analysis of the light curves revealed variability on a shorter timescale, this would shift the typical collision radius inwards where the densities are higher. As shown in Rudolph et al. (2022), this would lead to stronger signatures of secondary particles, such as secondary leptons and neutrinos, for the same baryonic loading.
The synthetic HE light curve (shown in the middle panel of Figure 2) is very similar to the keV-MeV light curve. Preliminary analysis of LAT data using gtburst 5 showed no evidence for a fourth peak in the 20 MeV − 300 GeV light curve. If this is later verified by detailed LAT analysis, then late-time collisions with low Lorentz factors will be needed in our model in order to suppress the late-time HE emission due to high internal opacity to γγ pair production. While the shape of the HE light curve is similar in the SYN-and Figure 2. Synthetic light curves for three energy ranges obtained in the "R16" scenario for the IC-and SYN-dominated cases (in all panels indicated as red and blue curves respectively, see legend in right panel). In the left panel we show only the ICdominated case as it resembles the SYN-dominated one. We also add observed light curves of Konus and INTEGRAL SPI-ACS that operate in a comparable energy range, re-normalised to match the same scale and shifted in time to match the same peak times (-175 s for Konus, and -225 s for INTEGRAL SPI-ACS). We set T obs = 0 for the synthetic light curve as the minimum collision time in the observer's frame.
IC-dominated scenarios, the HE flux is higher in the latter case, suggesting that primary electrons are cooling more efficiently via ICS (see also next paragraph). Nonnegligible VHE emission between 1-10 TeV is also expected in the IC-dominated scenario (right-hand panel), with similar light curve as the in the lower energy bands.
Photon spectra. Fig. 3 shows the predicted timeintegrated photon fluences as a function of observed energy; the inset is a zoom-in on the spectra around the peak including the photon index (that is defined as the spectral slope of dN ph /dE obs ).
We first discuss the results for the SYN-dominated cases (blue lines in both panels). The spectrum around the peak is by construction independent of the typical collision radius and given by the standard synchrotron fast-cooling predictions: approaching a photon index −3/2 below the peak and −(p e + 2)/2 above the peak. Differences between "R 17 " and "R 16 " are visible at low and high energies with respect to the MeV peak. For "R 16 ", where the densities are higher, emission of secondaries increases the fluence in the eV range, leading to a low-energy spectral break. However, the spectrum in the LAT band is still dominated by the synchrotron emission of fast-cooling primaries. For "R 17 ", the powerlaw spectrum with photon index −1.5 extends to the lowest energies (till synchrotron self-absorption becomes important), while higher maximal synchrotron energies yield a harder high-energy spectrum.
In the "R 16 " IC-dominated case (red curves in left panel), the spectral slope below the peak is still −1.5, but the HE emission (in the LAT range) is modified by the ICS emission of primary and secondary leptons that creates an almost flat spectrum (i.e. photon index ∼ −2). In this case, the peak of the broadband spec-trum is shifted to ∼ 1 GeV. In the "R 17 " IC-dominated case (red curves in right panel) the photon index below the MeV peak (but not within the GBM band) becomes asymptotically ∼ −1.25, which is the highest value that can be produced in our model. This is indicative of electron cooling via ICS in the Klein-Nishina regime (Daigne et al. 2011). In the LAT band the spectrum is a powerlaw that extends to about 10 TeV (where the Klein-Nishina ICS emission of primary electrons with Lorentz factor γ e,min dominates). However, due to EBL attenuation this spectral feature is washed out. We note that in the "R 17 " IC-dominated case the contribution of secondary leptons in the GBM and LAT bands is negligible. For a more detailed discussion on spectra and their decomposition, we point to Sec. 4 and 5 in Rudolph et al. (2022).
In addition to the full time-integrated spectra, we computed the spectra for the three main emission periods, namely the first two bright pulses and the late-time pulse (not explicitly shown). Finding little difference in the spectra for these three pulses we predict no significant spectral evolution during these three emission periods.
Neutrinos. We include in Fig. 3 the predicted neutrino spectra (per flavor) 6 and the corresponding Ice-Cube limits (The IceCube Collaboration 2022) for the different scenarios as dashed curves. We also compute the number of expected neutrino events in IceCube with the appropriate point-source effective area for the declination range 18-21 degrees (Abbasi et al. 2021), and find n νµ = 0.012 (0.006) for the "R 17 " SYN-(IC-) dom- Figure 3. Modelled spectra E obs FE obs for the (left) "R16" and "R17" (right) scenarios, in both cases including a SYN-and a ICdominated scenario (for parameters see Table 1). Transparent curves correspond to the fluence without taking EBL attenuation into account, dashed curves to the per-flavour neutrino fluences. The latter are compared with the IceCube upper limit reported in The IceCube Collaboration (2022). The inset shows a zoom-in to the spectrum around the peak in the Fermi GBM/LAT range (indicated as a shaded region in all plots), with a dotted vertical line at 1060 keV to indicate the position of the observed peak. For this energy range we further indicate the photon indices, that are defined of the slope of dN ph /dE obs . inated case, and n νµ = 0.17 (0.29) for the "R 16 " SYN-(IC-) dominated case. The predicted neutrino fluences are thus below the IceCube limits in all cases and consistent with non-detection, as a result of the relatively large R Coll paired with the chosen baryonic loadings in consistency with the findings in one zone models (Ai & Gao 2022;Murase et al. 2022). For "R 16 ", however, a baryonic loading of f p/e 3 is expected to be in tension with neutrino limits, and if observations eventually favor a variability timescale δt var 1 s the baryonic loading would be limited to an even smaller value.
Contrary to Rudolph et al. (2022), we normalise the initial engine kinetic energy to achieve the same γ-ray isotropic energy, which for the IC-dominated scenarios increases the required energy (see Tab. 1). This increases the energy transferred to non-thermal protons, and subsequently the neutrino fluences. The effect can be noticed clearly for the "R 16 "-scenario. For the "R 17 "scenario the neutrino production efficiency is limited by the low(er) maximal proton energies, due to the lower magnetic fields obtained in the IC-dominated case. It is interesting that the peak neutrino energies of 10 17 to 10 19 eV exceed the expectation of the standard neutrino model for GRBs (10 15 eV, see e.g. Hummer et al. (2012)) by at least two orders of magnitude in energy. This is a result of the synchrotron-cooling dominated spectral indices below the peak (the photon number density peaks at lower energies) paired with weak magnetic field effects on the secondaries as a consequence of large R Coll . Therefore, energetic GRBs may be a target for future radio detection experiments.

DISCUSSION
There are a number of effects which can be included in order to enrich model. For example, the different pulses may be produced by collisions of shells with very different Lorentz factor ranges or even microphysics parameters. This can cause a combination of SYN-and ICdominated scenarios in different peaks (see also Zhang et al. 2022, for the reverse shock), or suppression of VHE emission in others (by low Lorentz factors enhancing the γγ optical thickness, and also the neutrino production). One may speculate that a non-observation of the late-term peak by Fermi-LAT or a contribution to the LHAASO signal could be produced by such effects. Within each peak, low Lorentz factors and a strong correlation between collision radius and observation time can cause an early suppression VHE photons, which may be interpreted as a delay (Bustamante et al. 2017). While our spectral index above the peak can be adjusted by the electron injection index and the efficiency of ICS, our model may not accommodate very hard low-energy photon indices of ∼ −1, unless additional components, such as photospheric thermal emission, are considered. In lepto-hadronic models these components would increase the number of target photons available for photohadronic interactions. We note that the highest photon index of −1.25 was obtained in the "R 17 " IC-dominated scenario, but well below the GBM band. Whether such a high dissipation radius is indeed realistic should be confirmed by observations of the variability timescale and estimates of the Lorentz factor.
The GRB central engine may be a newly formed accreting black hole or magnetar. For energetic bursts, such as GRB 221009A, the latter scenario can be excluded because the available energy is limited by the magnetar's rotational energy to ∼ 2 · 10 52 erg (Usov 1992;Thompson et al. 2004). The rotational energy, E rot , of an accreting black hole can be extracted via electromagnetic fields, provided there is a strong largescale magnetic field threading the black hole horizon (Blandford-Znajek mechanism (BZ) Blandford & Znajek 1977). The isotropic equivalent energy of the jet can be written as E jet,iso = η j f −1 b E rot , where η j < 1 is the fraction of rotational energy ending up in the jet, f b = 1 − cos(θ j ) ≈ θ 2 j /2 is the beaming factor, θ j is the jet half-opening angle, E rot = f (a)M BH c 2 , f (a * ) = 1− (1 + 1 − a 2 * )/2, and a * is the dimensionless black hole spin. Adopting θ j = 3.5 deg (D'Avanzo et al. 2022) 7 and η j = 0.5 we find that a maximally spinning (a * = 1) black hole with M BH = 10 M can produce E jet,iso 1.4 · 10 57 erg, which is comfortably larger than E kin,ini required for the "R 16 " scenario. Moderate spins (a * ∼ 0.5) would require M BH ∼ 40 M to meet the model's energetic requirements. The formation of such massive black holes in the collapsar scenario for long-duration GRBs (with stellar masses < 40 M in the zero-age main sequence) is not expected (Woosley 1993) -see, however, Siegel et al. (2021) for very massive collapsars. The "R 17 "-scenario with baryonic loading 30 corresponds to E kin,ini ≥ (3 − 5) · 10 57 erg which would need a * = 1 and M BH ≥ 20 M . Hence, this scenario is unlikely on energetic grounds. We point out that although a BZ-powered jet would be initially Poynting-flux dominated, our internal shock model implies a matter-dominated jet at the dissipation region (see e.g. Granot et al. 2011;Giannios & Uzdensky 2019;Gottlieb et al. 2022, for possible scenarios for energy conversion in jets). For a non-negligible magnetisation at the dissipation radius, a different mechanism for particle energization may instead be invoked, such as magnetic reconnection (e.g. Sironi & Spitkovsky 2014;Guo et al. 2016;Werner et al. 2018).
The energy budget is further typically constrained from the afterglow brightness. Here it is noteworthy to mention that a part of the energy in the afterglow being dissipated in VHE may relax the requirements on the prompt-phase efficiency, especially for protonsynchrotron models for the VHE emission (Isravel et al. 2022). Note that a part of the afterglow energy going into thermal particles may equally increase the allowed kinetic energy of the blastwave after the prompt phase. Afterglow observations are also used to infer typical parameters of the GRB, for example Ren et al. (2022) find a Lorentz factor Γ 0 = 190 (with room for slightly higher Γ 0 , see their Fig. 3) for the afterglow. This would favour our "R 16 "-scenario, that has a typical Lorentz factor Γ fin ∼ 230 after the prompt phase. Our obtained (averaged) maximal proton energies are in the range between 10 20 and 2 · 10 21 eV under the assumption of efficient particle acceleration, see Table 1. Higher proton energies are expected in the SYNdominated scenarios than the IC-dominated ones because the acceleration rate is higher due to higher magnetic fields. Our maximal energies are compatible with the values used e.g. in Das & Razzaque (2022) (100 EeV) to describe the LHAASO VHE photons from EBL interactions. We note that the time delay induced by the extragalactic magnetic fields (EGMFs) requires extremely low field values paired with large proton energies, which means that this challenge can be somewhat mitigated in our SYN-dominated scenario by the high proton energies. We also note that the EGMF induced delay is very large (assumed to be limited to the observed LHAASO window of 2000s) even under aggressive assumptions for the magnetic field, which means that the protons must be accelerated in the prompt phase of the GRB and further delays induced by the afterglow cannot be accommodated. We find that most of the UHECR protons are emitted within the first 100s in our model, which is compatible with this picture. Our obtained maximal proton energies are also higher than the maximum corresponding rigidity R max 1 − 3 EV required to describe UHECR data (Heinze et al. 2019) -where details depend on the assumed cutoff shape.

SUMMARY AND CONCLUSIONS
In this letter we have presented a state-of-the-art multi-messenger emission model for the prompt emission of GRB 221009A. In this model, plasma shells are ejected from a central engine with varying Lorentz factors; these eventually catch up and energy is dissipated in internal shocks. Our radiation model includes the effect of UHECR protons self-consistently, such as the electromagnetic cascade in each shell, generated from secondary electrons, positrons and photons produced in photo-hadronic interactions. Our assumptions for the baryonic loading (3 and 30) have been motivated by the paradigm that energetic GRBs, such as GRB 221009A, could be sources of the UHECRs; they are also consistent with the hypothesis that the highest-energy LHAASO photons come from interactions of UHECRs with the EBL.
We have demonstrated that an intermittent engine can reproduce the observed prompt light curves if a quiescent period of about 200 seconds is included and assuming a variability timescale of ∼ 1s. We have implemented relatively large Lorentz factors and therefore collision radii supported by various arguments (such as γγ optical thickness, neutrino non-observation). Our predicted electromagnetic spectra exhibit synchrotron fast cooling dominated spectral indices below the MeV peak, except for an ICS dominated scenario with energy dissipation at high radii that yielded a softer spectral index. Above the peak, ICS can affect the spectral index very strongly, an effect which can extend up to the highest energies and enhance the VHE fluxes. On the other hand, hadronic contributions did not significantly alter the photon spectra. Since pile-up effects affect the analysis for nearly all instruments due to the high brightness of the GRB, the predictive power of our model may be useful. The predicted neutrino emission was consistent with the non-observation of neutrinos by IceCube due to high predicted peak energies in the range interesting for radio neutrino telescopes. For lower emission radii (due to variability on shorter timescales and/or lower Lorentz factors) the baryonic loading would likely be constrained to lower values than the ones assumed here. Our findings are consistent with the available rotational energy which can be extracted from a maximally spinning black hole with a mass of the order of 10 M if the baryonic loading is not too far away from energy equipartition; our standard assumption of a baryonic loading of 3 is consistent with this picture.
We conclude that GRB 221009A is an interesting object to test the internal shock model and the paradigm that energetic GRBs could be the sources of UHECRs. While direct signatures of cosmic rays (such as neutrinos) have not been seen, the LHAASO observation of TeV photons could point towards UHE proton acceleration. Our model connects the different messengers for the prompt phase of this GRB.
Note. -During completion of this work Liu et al. (2022) appeared. In contrast to their paper, we selfconsistently model the photon spectra and account for several emission regions along the jet. Their low(er) emission radii are the consequence of a variability timescale of 82 ms that was derived from an analysis of the GBM light curve up to 219 s, whereas we inferred the variability timescale from the INTEGRAL SPI-ACS light curve also including the bright emission period most relevant for the spectra. Overall, their limits on the baryonic loading are compatible with our findings.