Deciphering the ∼18 TeV Photons from GRB 221009A

On 2022 October 9, at T0 = 13:16:59.000 UT (Veres et al. 2022), a long-duration gamma-ray burst (GRB) identified as GRB 221009A (also known as Swift J1913.1+1946) was detected in the direction of the constellation Sagitta by the Gamma-ray Burst Monitor (GBM; Meegan et al. 2009) on board the Fermi Gamma-ray Space Telescope. The prompt emission was also detected by several other space observatories, such as the Fermi Large Area Telescope (LAT), Swift (Dichiara et al. 2022; Krimm et al. 2022), AGILE (Piano et al. 2022; Ursi et al. 2022), INTEGRAL (Gotz et al. 2022), Solar Orbiter (Xiao et al. 2022), SRG (Lapshov et al. 2022), Konus (Frederiks et al. 2022), GRBAlpha (Ripa et al. 2022), and STPSat-6 (Mitchell et al. 2022). The GRB 221009A is located at the coordinate R.A. = 288.282 and decl. = 19.495 (Pillera et al. 2022). Fermi-LAT detected the most energetic photon of energy, 99.3 GeV (at t₀ + 240 s). It is the highest-energy photon ever detected by Fermi-LAT from a GRB in the prompt phase (Bissaldi et al. 2022; Pillera et al. 2022). The afterglow emission was also observed at different wavelengths (Das & Razzaque 2022), and the optical follow-up observation estimated a very low redshift of z ≃ 0.151 (de Ugarte Postigo et al. 2022). The total emitted isotropic-equivalent gamma-ray energy from GRB 221009A is estimated to be (2–6) × 10⁵⁴ erg (de Ugarte Postigo et al. 2022; Kann & Agui 2022). This is the brightest long-duration GRB and arguably one of the nearest, and possibly the most energetic, GRBs ever observed. It has also been reported that GRB 221009A produced a significant ionization of the Earth's lower ionosphere (∼60–100 km; Hayes & Gallagher 2022) and is the strongest ionization effect ever recorded from a GRB.

The Large High Altitude Air Shower Observatory (LHAASO) with the water Cherenkov detector array (WCDA) and the larger air shower kilometer square area (KM2A) detector observed more than 5000 very high energy (VHE) photons within T₀ + 2000 s in the 500 GeV–18 TeV energy range, making them the most energetic photons ever observed from a GRB (Huang et al. 2022). Surprisingly, the ground-based Cherenkov detector Carpet-2 at Baksan Neutrino Observatory reported the detection of what is undoubtedly a very rare air shower originating from a 251 TeV photon 4536 s after the GBM trigger from the direction of GRB 221009A (Dzhappuev et al. 2022). Observations of these unusually VHE gamma rays by LHAASO and Carpet-2 from GRB 221009A are incomprehensible and led to speculation about nonstandard physics explanations of these observed events. However, there is a caveat concerning the observation of the 251 TeV gamma ray. The angular resolution of Carpet-2 is several degrees, and the two previously reported Galactic VHE sources, 3HWC J1928+178 and LHASSO J1929+1745, are located close to the position of GRB 221009A (Fraija & Gonzalez 2022). It remains uncertain whether the observed 251 TeV photon is from GRB 221009A or either of these Galactic sources. Nevertheless, the temporal and spatial coincidence of this event with GRB 221009A is worth exploring (Finke & Razzaque 2022; Mirabal 2023; Alves Batista 2022). In the present context, we will delve into the VHE emission observed by LHAASO.

The VHE gamma rays observed by the Cherenkov telescopes from the extragalactic sources undergo energy-dependent attenuation by interacting with the extragalactic background light (EBL) through electron–positron pair production (Stecker et al. 1992; Ackermann et al. 2012). As a result, the shape of the spectrum at VHEs changes significantly. Several well-known EBL models have been developed to study the attenuation at different redshifts.These models have been used successfully by highly sensitive imaging atmospheric Cherenkov telescopes such as VERITAS (Holder et al. 2009), HESS (Hinton 2004), and MAGIC (Cortina 2005) to analyze the observed VHE gamma rays from sources of different redshifts. The observed VHE gamma-ray flux from the source can be written in terms of the intrinsic flux F_in and the survival probability of the VHE photon as (Hauser & Dwek 2001)

$$\begin{eqnarray}&&{F}_{\gamma }({E}_{\gamma })={F}_{\mathrm{in}}({E}_{\gamma }){e}^{-{\tau }_{\gamma \gamma }({E}_{\gamma })},\end{eqnarray} \tag{ 1 }$$

where E_γ is the observed VHE photon energy, and τ_{γ γ} is the optical depth for the pair-production process. The optical depth for an 18 TeV photon at a redshift of z = 0.151 is 18.3 in the EBL model of Franceschini et al. (2008) and 19.4 in the EBL model of Dominguez et al. (2011), which corresponds to the survival probability of the VHE photon, ${e}^{-{\tau }_{\gamma \gamma }}\sim 1.1\times {10}^{-8}$ and 3.6 × 10⁻⁹, respectively, in both of these models. Thus, for an 18 TeV photon energy, the observed flux will be suppressed by a factor of ∼10⁻⁹–10⁻⁸. As the observation of an 18 TeV photon from a source at redshift z = 0.151 is difficult to comprehend, it is viewed as a signature of new physics, such as Lorentz invariance violation (Li & Ma 2022; Baktash et al. 2022; Finke & Razzaque 2022; Zhu & Ma 2022), oscillation of a photon to a pseudoscalar particle (axion-like particle, ALP; Galanti et al. 2022; Lin & Yanagida 2022; Troitsky 2022), ALP abundance enhanced with its mass caused by a first-order phase transition in a hidden sector (Nakagawa et al. 2022), heavy neutrinos as the means of propagation to avoid the energy attenuation (Cheung 2022), and sterile neutrinos produced via mixing with active neutrinos (Brdar & Li 2022). On the other hand, from the standard physics point of view, these gamma rays are argued to be the secondaries arising from the interactions between the ultrahigh-energy cosmic rays emanating from GRB 221009A and the cosmological photon background on their way to the Earth (Alves Batista 2022). Also, observation of neutrinos from such a bright GRB is discussed (Murase et al. 2022).

Since the VHE spectra of most of the flaring high-energy peaked BL Lac blazars (HBLs) of different redshifts are explained very well using the EBL models of Franceschini et al. (2008) and Dominguez et al. (2011), the recent observation of an ∼18 TeV photon from GRB 221009A falls short of this expectation. The obvious question is, can it be due to the intrinsic flux from the source? If we look into Equation (1), the depletion in the flux due to ${e}^{-{\tau }_{\gamma \gamma }}$ can, in principle, be compensated for by increasing the intrinsic flux. However, this may not be possible in most of the situations. As noted previously, GRB 221009A is very special, as its afterglow outburst outshone all other GRBs seen before, despite the fact that GRB 221009A is obstructed by the Milky Way. Furthermore, the burst was so powerful that it ionized Earth's atmosphere and disrupted longwave radio communications. It is estimated that, at low redshifts, such energetic GRBs are extremely rare events and may occur once in a century (Atteia 2022). Thus, it is possible that the intrinsic VHE flux from the source might have increased manifolds, which could compensate for the depletion from the EBL effect. In this letter, we would like to pursue such a scenario and its impact on the observation of ∼18 TeV photons by LHAASO.

The emission mechanisms in blazars (a subclass of active galactic nuclei, AGNs) and GRBs have many features in common (Urry & Padovani 1995; Gehrels & Razzaque 2013). Such common features are found to prevail in the synchrotron luminosity and Doppler factor between GRBs and AGNs (Wu et al. 2011). In several studies, it was observed that the jets in blazars and GRBs share common features despite large differences in their masses and bulk Lorentz factors (Wang & Wei 2011; Nemmen et al. 2012; Wu et al. 2016). It is natural to use such mechanisms and processes to study the multi-TeV flaring of HBLs to study the afterglow phases of GRBs.

Previously, we have used the photohadronic process to study the multi-TeV flaring from HBLs (Sahu 2019; Sahu et al. 2019, 2020). In the photohadronic scenario, protons in the blazar jet are accelerated to VHEs and then collide with the background seed photons to produce Δ-resonance (p γ → Δ⁺) with the following kinematical condition (Sahu 2019):

$$\begin{eqnarray}&&{E}_{p}{\epsilon }_{\gamma }=0.32\,{\rm{\Gamma }}\,{ \mathcal D }{(1+z)}^{-2}\,{\mathrm{GeV}}^{2},\end{eqnarray} \tag{ 2 }$$

where E_p and _γ are the proton energy and background seed photon energy, respectively, in the observer's frame. In the process, the observed VHE photon carries about 10% of the proton energy, E_γ ≃ 0.1Ep. The bulk Lorentz factor and the Doppler factor, respectively, are given by Γ and ${ \mathcal D }$. As the jets of the observed HBLs and GRBs beam toward us, ${\rm{\Gamma }}\simeq { \mathcal D }$. The Δ-resonances decay to neutral pions that subsequently decay to VHE gamma rays. These are the blueshifted photons observed by the Cherenkov telescopes on Earth. This model is very successful in explaining the VHE gamma-ray spectra from several HBLs, and the intrinsic flux F_in is given by

$$\begin{eqnarray}&&{F}_{\mathrm{in}}={F}_{0}\,{E}_{\gamma ,\mathrm{TeV}}^{-\delta +3},\end{eqnarray} \tag{ 3 }$$

where E_γ,TeV is the photon energy in TeVs. The normalization constant F₀ can be fixed from the observed spectrum, and the spectral index δ = α + β is the free parameter in the model (Sahu 2019; Sahu et al. 2019). Note that F_in is independent of Γ and ${ \mathcal D }$. The high-energy protons in the jet have a power-law differential spectrum ${dN}/{{dE}}_{p}\propto {E}_{p}^{-\alpha }$, where E_p is the proton energy, and we take α = 2 (Dermer & Schlickeiser 1993), a generally accepted value. For HBLs, the seed photon flux also follows a power law ${{\rm{\Phi }}}_{\gamma }\propto {\epsilon }_{\gamma }^{\beta }\propto {E}_{\gamma }^{-\beta }$ (Sahu 2019; Sahu et al. 2019). For HBLs, the value of δ always lies in the range 2.5 ≤ δ ≤ 3.0, which corresponds to a β value in the range 0.5 ≤ β ≤ 1.0, indicating that the seed photons are in the low-energy tail region of the synchrotron self-Compton (SSC) spectrum (Sahu et al. 2019). Recently, it has been shown that for GRBs, the value of β can be positive or negative (Sahu & Fortin 2020). Here β > 0 implies that the seed photons are in the self-Compton regime, and β < 0 locates the seed photons in the synchrotron regime. It was previously shown that the VHE spectra of GRB 190114C and GRB 190829A are due to the interaction of the high-energy protons with the low-energy tail region of the background SSC photons in the jet with β > 0 (Sahu & Fortin 2020; Sahu et al. 2022). Also shown was that the VHE spectrum of GRB 180720B is from the interaction of high-energy protons with the synchrotron seed photons in the jet environment with β < 0 (Sahu & Fortin 2020). This negative value of β corresponds to the falling part of the synchrotron spectrum.

With its two detectors, WCDA and KM2A, LHAASO detected ≥5000 photons above 500 GeV from GRB 221009A within T ∼ 2000 s of the prompt emission. The number of photons N_γ detected at a time interval T by any of these detectors at zenith angle θ and effective area A(E_γ , θ) is (Zhao et al. 2022)

$$\begin{eqnarray}&&{N}_{\gamma }=T{\int }_{0.5\,{TeV}}\displaystyle \frac{{{dN}}_{\gamma }}{{{dE}}_{\gamma }}A({E}_{\gamma },\theta ){e}^{-{\tau }_{\gamma }({E}_{\gamma })}{{dE}}_{\gamma },\end{eqnarray} \tag{ 4 }$$

where the differential photon spectrum can be written as 

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{\gamma }}{{{dE}}_{\gamma }}={F}_{0}\,{E}_{\gamma ,\mathrm{TeV}}^{-\delta +1}\,{\mathrm{TeV}}^{-2}.\end{eqnarray} \tag{ 5 }$$

The source was observed at the zenith angle of 30° ≲ θ ≲ 35° that we adopt in Equation (4). Taking into account the areas of LHAASO-WCDA and LHAASO-KM2A (Cao et al. 2022), we evaluate the integral in Equation (4) for δ = 2.5, 1.7, and 1.2. For the present analysis, we consider the EBL model of Franceschini et al. (2008). We assume that these two detectors observe photons above 500 GeV in the range 5000 ≤ N_γ ≤ 6500. By fixing the value of N_γ , we calculate the value of F₀, which is then used to calculate the VHE photon flux and the integrated flux ${F}_{\gamma }^{\mathrm{int}}$ in the energy range 100 GeV ≤ E_γ ≤ 18 TeV.

In Figure 1, we show the predicted spectra for δ = 2.5, 1.7, and 1.2 by taking into account the effective area of the LHASSO-WCDA detector and fixing N_γ = 5500. The relative energy resolution of LHAASO-WCDA is ≃50% at energies around 18 TeV (Figure 26 of Chapter 1 of Cao et al. 2022). For δ = 2.5, the flux starts from a maximum value of F_γ ∼ 10⁻⁸ ergcm⁻² s⁻¹ at E_γ = 100 GeV and decreases slowly up to ∼4 TeV. Beyond ∼4 TeV, it falls faster due to the EBL effect. The spectrum intersects with the sensitivity curve of LHAASO with 2000 s exposure at E_cut = 9.94 TeV, which is at the lower edge of the energy resolution (at 9 TeV). The δ = 2.5 value implies β = 0.5 with the intrinsic flux ${F}_{\mathrm{in}}\propto {E}_{\gamma ,\mathrm{TeV}}^{0.5}$. This corresponds to seed photons in the lower tail region of the SSC spectrum in the GRB jet. The accelerated high-energy protons in the jet interact with these seed photons to produce VHE gamma rays, a situation very similar to the VHE flaring of HBLs.

Zoom In Zoom Out Reset image size

We repeat the calculation for δ = 1.7, which corresponds to β = −0.3. As discussed previously, the negative value of the seed photon spectral index β corresponds to photons in the descending part of the synchrotron spectrum toward higher _γ values and ${{\rm{\Phi }}}_{\gamma }\propto {\epsilon }_{\gamma }^{-0.3}$. Thus, in this case, the high-energy protons interact with the seed photons in the synchrotron regime of the external forward shock region to produce gamma rays. The spectrum starts with F_γ ∼ 10⁻⁹ erg cm⁻² s⁻¹ at E_γ = 100 GeV, increases very slowly up to ∼4 TeV, and then falls faster as the exponentially decaying term from the EBL dominates. The curve intersects with the LHAASO sensitivity curve at E_cut = 11.53 TeV. The intrinsic flux increases as ${F}_{\mathrm{in}}\propto {E}_{\gamma ,\mathrm{TeV}}^{1.3}$.

Finally, we consider a smaller value of δ = 1.2, which is shown in Figure 1. This value of δ gives β = −0.8. In the photohadronic context, this corresponds to ${{\rm{\Phi }}}_{\gamma }\propto {\epsilon }_{\gamma }^{-0.8}$, which is the descending part of the synchrotron spectrum toward higher _γ values like the ones for δ = 1.7. However, in this case, the seed synchrotron spectrum in the external forward shock region falls faster than the one for δ = 1.7. The spectrum increases and reaches a maximum flux at E_γ ∼ 4.5 TeV and then decreases exponentially for large values of E_γ intersecting the LHAASO curve at E_cut = 12.44 TeV. The intrinsic flux in this case behaves like ${E}_{\gamma ,\mathrm{TeV}}^{1.8}$.

We repeat the calculation by using the effective area of LHAASO-KM2A for δ = 2.5, 1.7, and 1.2 and N_γ = 5500. The results are shown in Figure 2. For 18 TeV photons, the relative energy resolution of LHAASO-KM2A is ≃36% (Figure 2 of Chapter 1 of Cao et al. 2022), which puts the observed photon energy in the range 11.52–24.48 TeV. For a given δ, the pattern of the spectra in LHAASO-WCDA and LHAASO-KM2A is similar in both cases, but the E_cut value for LHAASO-KM2A is smaller than that of LHAASO-WCDA. Also, the E_cut for LHAASO-KM2A is less than 11.52 TeV, which shows that LHAASO-KM2A may not be able to detect these photons.

Zoom In Zoom Out Reset image size

To account for the behavior of the VHE spectrum for different N_γ , we fix N_γ = 5500 and 6500 to calculate F₀, the integrated flux ${F}_{\gamma }^{\mathrm{int}}$, and the luminosity L_γ and E_cut values using the effective areas of the LHAASO-WCDA and LHAASO-KM2A detectors. These are given in Table 1. The results of LHAASO-KM2A are the bracketed values in Table 1. It can be seen that increasing N_γ from 5500 to 6500 leads to an increase in all quantities. This implies that by knowing N_γ and the maximum value of E_γ , we can predict the VHE gamma-ray spectrum, provided the EBL contribution is well understood.

Table 1. Estimate of Observables Using the Effective Areas of the Detectors WCDA and KM2A

δ	Nγ	F0	${F}_{\gamma }^{\mathrm{int}}$	Lγ,48	Ecut
2.5	5500	2.11 (0.54)	2.49 (0.63)	1.63 (0.41)	9.94 (8.17)
	6500	2.50 (0.63)	2.95 (0.75)	1.92 (0.49)	10.18 (8.41)
1.7	5500	1.41 (0.48)	1.22 (0.41)	0.80 (0.27)	11.53 (10.48)
	6500	1.67 (0.56)	1.44 (0.49)	0.94 (0.32)	11.70 (10.63)
1.2	5500	0.92 (0.36)	1.07 (0.42)	0.70 (0.27)	12.44 (11.32)
	6500	1.08 (0.42)	1.26 (0.50)	0.83 (0.32)	12.55 (11.50)

Note. Using the LHAASO-WCDA (30° ≤ θ ≤ 45°) and LHAASO-KM2A effective detector areas, as well as different values of δ and numbers of events N_γ , we have calculated the flux normalization factor F₀ in units of 10⁻⁸ erg cm⁻² s⁻¹ and the integrated flux ${F}_{\gamma }^{\mathrm{int}}$ in units of 10⁻⁸ erg cm⁻² s⁻¹ in the energy range 100 GeV–18 TeV, along with the corresponding luminosity L_γ,48 in units of 10⁴⁸ erg s⁻¹. Here E_cut is the value of E_γ in TeV where it intersects with the LHAASO sensitivity curve with 2000 s exposure time. The bracketed values are the results using the LHAASO-KM2A detector area.

Download table as:  ASCIITypeset image

For a given value of N_γ , the E_cut value increases and approaches ∼18 TeV as δ decreases from 2.5 to 1.2. Moreover, by further decreasing δ, one can reach E_cut ∼ 18 TeV, which corresponds to a very stiff synchrotron spectrum and may be problematic. Also, for a given δ, by increasing N_γ , the E_cut value increases. From our analysis, we observed that LHAASO-WCDA is more likely to observe photons of energy ∼18 TeV than LHAASO-KM2A. From the dependence of E_cut on δ, we infer that the interaction of high-energy protons with the descending part of the synchrotron seed photon spectrum is more likely to produce ∼18 TeV photons than the high-energy protons interacting with the low-energy tail region of the seed SSC photons in the GRB jet.

Additionally, we calculate the chance probability of N_γ ≥ 1 for E_γ ∼ 18 TeV by taking the normalization constant F₀ as a variable in the range 10⁻¹⁰ ≤ F₀ (erg cm⁻² s⁻¹) ≤ 2 × 10⁻⁷ for a fixed δ = 1.2. We fix the total number of observed photons above 500 GeV to be 5500 and T = 2000 s. We assume that the noise in the data has a Gaussian distribution with an unknown standard deviation of σ and variance σ² (Gregory 2010). Bayesian inference is implemented by using the Markov Chain Monte Carlo method to estimate the posterior probability distribution function (PDF) as a function of F₀ for a given σ value. This we have done for σ in the range 100 ≤ σ ≤ 1500. Using the PDFs for different σ values and E_γ ∼ 18 TeV, we evaluate

$$\begin{eqnarray}&&{N}_{\gamma }({E}_{\gamma })=T{\int }_{(1-{\delta }_{E}/2){E}_{\gamma }}^{(1+{\delta }_{E}/2){E}_{\gamma }}\frac{{{dN}}_{\gamma }}{{dE}{{\text{'}}}_{\gamma }}A(E{\text{'}}_{\gamma },\theta ){e}^{-{\tau }_{\gamma }(E{\text{'}}_{\gamma })}{dE}{{\text{'}}}_{\gamma },\end{eqnarray} \tag{ 6 }$$

using the Monte Carlo simulation. Here δ_E = 0.5 is the value for the 50% uncertainty in the energy around 18 TeV for the LHAASO-WCDA. For a given value of σ, we repeat the procedure 10⁶ times. This way, we obtain the percentage of the chance probability of N_γ ≥ 1 = 10⁻⁴ × no. of times N_γ (18 TeV) ≥ 1. The results are plotted in Figure 3. It is observed that if the posterior PDF is symmetric around the midpoint, then the percentage of chance probability of detection of N_γ ≥ 1 is very small. Similarly, for an asymmetric PDF with the weight factor leaning more toward smaller values of F₀, the percentage of chance probability is also very small. However, for the asymmetric PDFs with the weight factor leaning more toward larger values of F₀, the percentage of chance probability is large, with a maximum value of 40%.

Zoom In Zoom Out Reset image size

In conclusion, the recent observation of ∼18 TeV photons by LHAASO from GRB 221009A sheds doubt on the applicability of the well-known EBL models for photons of energy >10 TeV at a redshift of ≳0.151, even though these EBL models work very well in explaining the VHE spectra of so many other TeV sources. This incompatibility has led toward new physics solutions. However, there is still a conventional way to delve into the problem, which we propose here. We argue that high-energy protons interacting with the synchrotron photon background in the GRB jet will be able to produce photons of energy close to 18 TeV. Assuming that the error in the data has a Gaussian distribution and using the area of LHAASO-WCDA, we obtain a maximum 40% chance probability of observing N_γ ≥ 1 for E_γ ∼ 18 TeV. Our analysis shows that LHAASO-WCDA is more likely to observe photons of energy ∼18 TeV than LHAASO-KM2A. We anticipate that the publication of the GRB 221009A results should be able to either confirm or rule out either most or all of the explanations discussed here.

The authors are indebted to Miguel Enrique Iglesias Martínez and José Guerra Carmenate for helping in the Bayesian analysis. We are also thankful to the anonymous referee for making useful suggestions. S.S. is thankful to S. Salokya for reading the manuscript. The work of S.S. is partially supported by DGAPA-UNAM (Mexico) project No. IN103522. B.M.-C. and G.S.-C. would like to thank CONACyT (México) for partial support. Partial support from CSU-Long Beach is gratefully acknowledged.