Introduction

Gamma-ray bursts (GRBs) are the most luminous stellar explosions in the universe. These events generally appear as brief and intense γ-rays followed by a long-lived afterglow emission. The GRB prompt emission originates from relativistic jets that dissipate the energy and accelerate particles either via internal shocks or magnetic reconnection, with high variability and usually lasts from milliseconds to thousands of seconds1,2,3. Most of the observed spectrum of GRB prompt emission in the keV to MeV energy range usually can be described by a smoothly joint broken power-law function (called the Band function4). Despite decades of intensive investigation, our understanding of the physics behind the prompt emission of GRBs remains limited.

The existence of X-ray or gamma-ray emission and absorption lines in the GRB energy spectrum has been debated. For example, in prompt emission phase, the Konus instrument detected absorption lines at 30–70 keV and emission lines at 400–460 keV in the energy spectra of some GRBs5,6. HEAO-1 observed absorption-like features in the spectra of some GRBs7,8. The Japanese Ginga Gamma-Ray Burst Detector (GBD) observed two absorption-like features in three GRBs (GRB 870303, GRB 880205 and GRB 890929), which may be interpreted as the first and second cyclotron absorption lines9,10,11. Additionally, some 6.4 keV iron K-α spectral lines were claimed to have been found in some bursts. Within the energy range of 3.8 ± 0.3 keV, a possible transient Fe absorption feature was identified in the prompt X-ray spectrum of GRB 990705. This feature appeared during the initial rising phase of the burst profile and disappeared thereafter12. In another study, Frontera et al.13 analyzed the prompt emission spectrum of GRB 011211 and found potential indications of transient Fe absorption features around 6.9 ± 0.6 keV during the rise of the main pulse. However, the statistical significance of these features is found to be below the 5σ threshold. Even when extending the spectral lines search to the afterglow phase and conducting large-scale searches using X-ray detection satellites such as Chandra14, Swift X-ray Telescope15, X-ray Multi-Mirror Mission Newton (XMM-Newton)16, Advanced Satellite for Cosmology and Astrophysics (ASCA)17, Satellite per Astronomia X (BeppoSAX)18, no credible X-ray line feature has been detected in GRBs afterglow19,20,21.

More recently, a highly significant (>5σ) narrow emission feature around 10 MeV has been detected in the Fermi data of GRB 221009A22,23. These intriguing features appear during the decay phase of the brightest pulse, with the central energy of the Gaussian distribution gradually shifting towards lower energies over time (about 37 MeV to 6 MeV), while the ratio of the line width to the central energy is nearly constant (about 10%). At the same time, in the Konus-Wind data of GRB 221009A, a similar narrow emission feature has been found with a significance level below 2σ24. Two independent satellites simultaneously detected the narrow emission feature in GRB 221009A, further bolstering the credibility of the narrow emission feature observed in GRB 221009A.

In this work, we perform a spectral analysis of the prompt emission from GRB 221023A using Fermi Gamma-Ray Burst Monitor (GBM)25 data. We find a marginally significant narrow emission feature around 2.1 MeV. The highest statistical significance of this feature is observed in the time interval 8–30 s, with the chance probability value  <2.56  × 10−5 (after accounting for the look-elsewhere effect), corresponding to a Gaussian-equivalent significance >4.20σ. We find that the relativistic hydrogen-like high-atomic-number ions entrained in the GRB jet can generate such narrow MeV emission lines through the de-excitation of excited electrons.

Results

Light curve and spectral analysis

GRB 221023A triggered the GBM onboard Fermi at 20:41:34.92 UT on 23 October 202226. Simultaneously, this event was also detected by Konus-Wind27 and AGILE (Astrorivelatore Gamma ad Immagini LEggero)28. The GBM light curve shows one bright peak with a total duration time T90 about 39 seconds (s) in the 50–300 keV energy band26, and the fluence reported in the Fermi-GBM catalog is F = 3.41 × 10−4 erg cm−2 in the energy range 10–1000 keV29. The Fermi Large Area Telescope (LAT)30 instrument was triggered during this event, and the highest-energy photon detected is a 17 GeV event with a 99% probability which is observed 576 seconds after the GBM trigger31.

Panels a, b of Fig. 1 presents light curves for GRB 221023A at different energy bands. We analyzed the spectral evolution of the GRB prompt emission by spectra in 5 adjacent time intervals (labelled A (0–5 s), B (5–8 s), C (8–30 s), D (30–36 s) and E (36–60 s)). The spectra of time intervals 0 –5 s, 5–8 s, 30–36 s and 36 − 60 s can be fitted with a Band function (see methods subsection Spectral fitting), detailed analysis results are summarized in Table 1. Interestingly, when fitting the spectra in time interval 8–30 s using the Band function, as shown in the a and b panels of Fig. 2, revealing a distinct narrow and bright emission feature between 1 MeV and 3 MeV. This narrow emission feature can be well modeled by adding a Gaussian component on top of the Band function (see methods subsection Spectral fitting), the best-fit parameter values are \(\alpha=-0.9{3}_{-0.01}^{+0.01}\), \({E}_{p}=891.0{7}_{-33.19}^{+3.03}\,{{\rm{keV}}}\), \(\beta=-2.6{5}_{-0.02}^{+0.05}\), \({E}_{{{\rm{gauss}}}}=2154.6{0}_{-65.07}^{+53.37}\,{{\rm{keV}}}\), \({\sigma }_{{{\rm{gauss}}}}=229.3{6}_{-45.29}^{+93.57}\,{{\rm{keV}}}\), the c and d panels of Fig. 2 displays the corresponding fitted counts rate and νFν spectrum. We also present the best-fit νFν model spectra in 5 adjacent time intervals in Fig. 3. In time interval 8–30 s, comparing models with and without the Gaussian component, we obtained ΔAIC = 51.87, \({{\rm{\ln }}}({{\rm{BF}}})=9.99\), and Δχ2 = 40.14, which strongly supports the presence of an additional narrow emission feature (see methods subsection Model comparison). In order to exclude the impact of background subtraction on the extracted spectrum, we calculated the background spectrum by selecting several different time windows. The result of the narrow emission feature is substantially unaffected (see methods subsection Background).

Fig. 1: Multiwavelength light curves and temporal evolution of spectral parameters.
figure 1

Panels (a) and (b) display the multi-energy band light curves of GRB 221023A observed by Fermi-GBM, with a bin size of 64 ms for each band. The time intervals for spectral analysis are indicated by vertical red dashed lines, labeled as A (0–5 s), B (5–8 s), C (8–30 s), D (30–36 s) and E (36–60 s). Panels (c), (d), and (e) show the temporal evolution of the low-energy spectral index α, peak energy Ep, and high-energy spectral index β of the Band model (deep blue points), respectively. Panels (f), (g) and (h) show the temporal evolution of the central energy Egauss, width σgauss, and Fluxgauss of the Gaussian component (purple points), respectively. All error bars represent 1σ uncertainties. Source data are provided as a Source Data file.

Full size image
Table 1 The spectral fitting results for each time interval of GRB 221023A
Full size table
Fig. 2: Spectral fitting within the 8 − 30 s time interval.
figure 2

The counts rate spectrum in the panel a and the νFν spectrum in the panel b are obtained from fitting the Band function. Data are from GBM's two sodium iodide (NaI) detectors (n0: dark purple, n1: blue) and one BGO detector (b0: light cyan-green). The narrow feature appears as an excess around 1 MeV–3 MeV in the b0 detector data. Panels c and d show the same spectra fitted with the Band function plus a Gaussian component to model the observed excess. Error bars indicate the 1σ uncertainty on data points. Source data are provided as a Source Data file.

Full size image
Fig. 3: Energy spectrum evolution.
figure 3

Best-fit νFν model spectra for the time-resolved data in different time intervals, five time intervals are color-coded, with the corresponding shaded colors show the 68% confidence levels.

Full size image

We performed a time-resolved spectral analysis on time interval 8–30 s to further investigate the presence of the observed feature and to characterize its evolution. We used a fixed window size of 13 s, sliding it in steps of 3 s to divide the time intervals, resulting in four subintervals, referred to as 8–21 s (C.1), 11–24 s (C.2), 14–27 s (C.3), 17 –30 s (C.4). In four subintervals, we still extracted the spectra by performing a different selection of the time windows for the background spectrum computation (see methods subsection Background). The narrow emission feature remains clearly visible in these finer time intervals. In the four finer time-resolved spectra (8 − 21 s, 11 − 24 s, 14–27 s, 17–30 s), the ΔAIC values vary between 25.76–36.55, the \({{\rm{\ln }}}({{\rm{BF}}})\) ranges from 2.06 to 7.34, and the Δχ2 ranges from 18.53 to 34.49. These results further strongly favor adding an additional narrow emission feature (see methods subsection Model comparison). The comprehensive results of the spectral analysis for these four time intervals are presented in Table 1. The temporal evolution of Gaussian component parameters is presented in the (f), (g), and (h) panels of Fig. 1. Notably, the central energy Egauss of the narrow emission feature remains constant at around 2.1 MeV, while the width σgauss shows a possible decreasing trend over time. The flux of the narrow emission feature is approximately 10−6 erg cm−2 s−1.

Significance analysis

In order to assess the significance of narrow emission features, we created 1.00 × 107 simulated data assuming the Band model and fitted them with both the Band and Band+Gaussian models to obtain the distribution of Δχ2. Table 2 shows the chance probability values (\({p\mbox{-}{{\rm{value}}}}_{{{\rm{sim}}}}\)) calculated based on the results of 1.00 × 107 simulations. We also considered the \({p\mbox{-}{{\rm{value}}}}_{{{\rm{sim}}}}\) corrected for the number of independent search trials (p-valuesim-trial). The highest statistical significance of narrow emission features was observed in the time interval C (8–30 s), with the chance probability value \({p\mbox{-}{{\rm{value}}}}_{{{\rm{sim}}}} < \;1\times 1{0}^{-7}\) obtained from results of 1 × 107 simulations, corresponding to a Gaussian-equivalent significance  >5.32σ. Considering the correction for the number of independent search trials, the chance probability value decreases to p-valuesim-trial < 2.56 × 10−5, corresponding to a Gaussian-equivalent significance  >4.20σ (see methods subsection Significance calculation of narrow emission feature). The chance probability values for the other time intervals are shown in Table 2.

Table 2 The results of evaluating the significance
Full size table

Comparison with GRB 221009A

The spectral analysis of GRB 221023A reveals a marginally statistically significant narrow emission feature at around 2.1 MeV. This would then represent the second event following GRB 221009A with a narrow emission feature in the MeV energy range. In the case of GRB 221009A, the central energy Egauss of the narrow emission feature decreases over time (about 37 MeV to 6 MeV), while the ratio of the line width to the central energy is nearly constant (about 10%)22,23. For GRB 221023A, we observe the trend: the central energy Egauss remains steady at around 2.1 MeV throughout the observation period, while the width σgauss exhibits a possible decreasing trend as time progresses, the flux of the narrow emission feature is around 10−6 erg cm−2 s−1. The Fig. 4 displays the lightcurves of GRB 221009A and GRB 221023A within the energy range of 0.2–40 MeV. The shaded regions indicate time intervals in which narrow emission features were detected. The narrow emission feature in GRB 221023A appears during the rising and falling phases of the brightest pulse, with a duration of 22 s (time intervals: 8–30 s) and then disappears. In contrast, the narrow emission feature in GRB 221009A appears during the falling phase of the brightest pulse, with a duration of 100 s (time intervals: 246–256 s and 270–360 s)22,23. Interestingly, the Fe absorption feature previously identified during the prompt emission of GRB 990705 and GRB 011211 appears during the rising phase of the main pulse12,13. This implies a higher likelihood of detecting emission or absorption features during time intervals near the peak of the main pulse in the prompt emission phase of GRBs. It is worth noting that due to the very high photon flux of GRB 221009A, the Fermi-GBM experienced Bad Time Interval (time interval affected by saturation) between 219 and 277 seconds32. It is possible that narrow emission features also exist during the rising phase of the main emission in GRB 221009A.

Fig. 4: Comparison of light curve.
figure 4

Comparison of the light curves of GRB 221023A (red) and GRB 221009A (blue) in the energy range of 0.2–40 MeV, with the corresponding shaded colors regions showing the time intervals where narrow emission features were detected. The trigger time of GRB 221023A was shifted backward by 210 s. Source data are provided as a Source Data file.

Full size image

Discussion

In general, standard models of prompt emission in GRBs do not predict the appearance of a transient MeV narrow emission component2,33,34. To explain our potential finding, we have explored several possible scenarios. One possible explanation for the narrow emission feature is the blue-shifted annihilation line of relatively cold (kBT ≪ mec2, where kB is Boltzmann constant, T is the temperature of the medium, me is the mass of the electron, and c is the speed of light) electron-positron pairs. Within the emission region (resulting from internal shocks and/or magnetic reconnection) of GRB, electron-positron pairs are readily formed within the GRB jet (such as two-photon pair production γγ ⟶ e+e−2,33). Numerical simulations of GRB spectra indicate that the generated spectra depend on the compactness of the fireball. In scenarios with high compactness, electron-positron pairs play an essential role in shaping the GRB prompt emission spectrum. In a pair-dominated fireball, a pair annihilation line is predicted35,36. In the observer frame, a line is expected to appear at an energy of E±,line = Γmec2/(1 + z), where Γ is the bulk Lorentz factor of the emitting region and z represents the redshift. For the typical energy Eline is about 2.1 MeV of the observed lines in GRB 221023A, the bulk Lorentz factor of the emitting region is required to be Γ is about 4(1 + z). Considering a redshift of z = 2, the bulk Lorentz factor Γ is about 12. In this scenario, how such a low bulk Lorentz factor is generated and maintained for an extended period in the prompt emission of the GRB is an issue.

The second scenario involves the possibility that the narrow emission feature is an intrinsic low-energy spectral line (such as the 6.4 keV fluorescent K-α iron line). This spectral line may be emitted within the region associated with the supernova ejecta. Subsequently, the energy of the spectral line could be boosted through up-scattering by the relativistic jet. The spectral line feature identified in GRB 221023A is narrow, which implies that electrons scattering photons are cold. This form of bulk Comptonization has already been proposed to occur within blazar jets37. The boosted photon energy of the low-energy spectral line is Eline = Γ2Elow/(1 + z), where Γ is the jet bulk Lorentz factor and Elow represents the low-energy spectral line of the particular element. If the observed spectral feature in GRB 221023A arises from the 6.4 keV iron K-α line, and the typical photon energy of the observed line is around 2.1 MeV, this would require a jet bulk Lorentz factor Γ to be about \(18{(1+z)}^{\frac{1}{2}}\). Considering a redshift of z = 2, the bulk Lorentz factor Γ is about 31. This scenario faces the same issue as the first one, namely, the low bulk Lorentz factor problem.

The third possible scenario is that the narrow emission feature may originate from MeV nuclear de-excitation lines. The energetic particles interacting with ambient matter could excite heavy nuclei which can emit MeV γ-ray line emissions via de-excitation, such as the 4.44 MeV line from 12C and the 6.13 MeV line from 16O38,39,40. In fact, nuclear de-excitation line emissions from 12C and 16O have been observed in solar flares41,42,43. Moreover, the existence of such nuclear de-excitation line has long been anticipated to be found within supernova remnants44,45,46,47. The observed photon energy of the nuclear de-excitation line is Eline = Eelement/(1 + z), where z is the redshift, and Eelement corresponds to the energy of the particular element’s nuclear de-excitation line. If we assume that the narrow emission feature in GRB 221023A arises from nuclear de-excitation lines of 12C or 16O, it corresponds to redshifts of z = 1.1 or z = 1.9, respectively. However, the radiation from the nuclear de-excitation of the ambient matter is almost isotropic, while the gamma-ray burst is collimated. Therefore, the energy budget needed to generate such a MeV emission line would be larger than the energy of the prompt emission for a typical half-opening angle θj = 0.1 rad. Note that for a large θj, such as θj > 0.3 rad, this energy problem can be alleviated.

The fourth possible scenario is that the heavy nuclei, especially relativistic hydrogen-like high-atomic-number ions originating from the β decay of unstable nuclei and/or the recombination, entrained in GRB jets can produce such a narrow MeV emission line via electron de-excitation48. In this model, the reflection of the radiation from the WR star can generate enough seed photons to excite electrons. After the jet with heavy nuclei breaks through the photosphere and hydrogen-like heavy ions are generated by β decay and/or recombination, the emission line can be generated. The emission lines occur at an energy of Eline = Γϵz/(1 + z), where ϵz = meα2c2Z2/2 is the Rydberg energy, α is the fine-structure constant, and Z = 29 is the atomic number of the copper. Note that we take copper as an example here since it can satisfy the half-life requirement of the model. This model can explain the MeV emission line of GRB 221023A well with reasonable parameters. We assume the half-opening angle θj of the GRB jet is about 0.1 rad and the redshift z = 0.1, we estimate that the jet beaming-corrected gamma-ray emission energy of GRB 221023A is about 7 × 1049 erg. Considering the propagation distance d  = 4.8 × 1011 cm, which corresponds to the timescale without an observable MeV emission line at the beginning, the Lorentz factor of the jet, Γ = (1 + z)Eline/ϵz, is about 300. The total mass of heavy nuclei entrained in the GRB jet is \({M}_{{{\rm{tot,nuclei}}}}={\theta }_{j}^{2}c{m}_{i}{{{\mathcal{E}}}}_{{{\rm{line}}}}/(4{E}_{{{\rm{line}}}}{\Gamma }_{e}d{\zeta }_{i})\) approximately 1026 g. Here, \({{{\mathcal{E}}}}_{{{\rm{line}}}}\) is the observed isotropic total energy of the emission line, Γe is the total excitation rate for an electron of the high-Z ion transitioned from the ground in the lab frame, mi is the mass of the high-Z ion (i.e., the copper), and ζi = 0.1 means the mass fraction capable of producing emission lines relative to the mass of all nuclei entrained in the GRB jet48. We find the kinetic energy of heavy nuclei in the jet is about 3 × 1049 erg, which is much larger than the energy budget of the observed emission line, approximately 3 × 1048 erg.

Methods

Fermi data analysis

The GBM consists of 12 sodium iodides (NaI) detectors (8 keV–1 MeV) and two bismuth germanate (BGO) detectors (20 keV –40 MeV)25, which has three different data types: continuous time (CTIME), continuous spectroscopy (CSPEC) and time-tagged event (TTE). The CTIME data include eight energy channels and have a finer time resolution of 64 ms. The CSPEC data include 128 energy channels, with a time resolution of 1.024 s. The TTE data consits of individual detector events, each tagged with arrival time, energy (128 channels), and detector number25. We download the GBM data of GRB 221023A from the public science support center at the official Fermi Website https://heasarc.gsfc.nasa.gov/FTP/fermi/data/gbm/triggers/2022/bn221023862/.

We extracted spectrum by using the TTE data from the brightest (with the smallest angle between this detector and the source object) two NaI detectors (n0, n1) and one BGO detectors (b0). The light curves were extracted using the GBM Data Tools49. The spectral analysis of the Fermi-GBM data was performed using the Bayesian approach package, namely the Multi-Mission Maximum Likelihood Framework (3ML)50. We selected the GBM spectrum over 8–900 keV and 0.3–30 MeV for NaI detectors and BGO detector, respectively. In order to avoid the iodine K-edge at 33.17 keV25, we ignore the data for the 30–40 keV energy ranges. The background spectrum from the GBM data was extracted from the CSPEC data with two time intervals before and after the prompt emission phase and modeled with a polynomial function of order 0–4 (Selected background time intervals:  − 130–−10 s, 100–200 s). We have used the Bayesian fitting method for the spectral fitting, and the sampler is set to the dynesty-nested. And we accounted for intercalibration constant factors among NaI and BGO detectors.

Spectral fitting

Figure 1 presents the light curves for GRB 221023A at different energy band. We subdivided the light curve into five intervals labeled A (0–5 s), B (5–8 s), C (8–30 s), D (30–36 s) and E (36–60 s), respectively, which were separated by red dashed vertical lines. We fit the corresponding spectra using the empirical Band function4, formulated as follows:

$${N}_{{{\rm{Band}}}}(E)=K\left\{\begin{array}{ll}{\left(\frac{E}{100\,{{\rm{keV}}}}\right)}^{\alpha }\exp \left(-\frac{E(2+\alpha )}{{E}_{p}}\right),\hfill&\left({{\rm{if}}}\,E \, < \, (\alpha -\beta )\frac{{E}_{p}}{2+\alpha }\right).\\ {\left[\frac{(\alpha -\beta ){E}_{p}}{(2+\alpha )\,100\,{{\rm{keV}}}}\right]}^{\alpha -\beta }\exp (\beta -\alpha ){\left(\frac{E}{100\,{{\rm{keV}}}}\right)}^{\beta },&\left({{\rm{if}}}\,E\ge (\alpha -\beta )\frac{{E}_{p}}{2+\alpha }\right).\end{array}\right.$$
(1)

where K is the normalization of Band spectrum, α and β are the low and high-energy photon spectral indices, respectively. E is the observational photon energy, and Ep is the peak energy of the νFν spectrum. The maximum values of the marginalized posterior probability densities and the corresponding 1σ uncertainties for each parameter of the Band model in each time interval are presented in Table 1.

The intriguing aspect was the shape of the spectrum in the time interval 8–30 s, as shown in the a and b panels of Fig. 2, revealing a distinct narrow and bright emission feature between 1 MeV and 3 MeV, this feature did not appear in the other four spectra. We further analyzed the GRBs data detected by the Fermi satellite within ten days before and after the explosion of GRB 221023A. For each of these events, we performed time-resolved spectral analysis using different signal-to-noise ratios and Bayesian blocks, no similar narrow feature were found in these GRBs. In order to model the narrow emission feature observed at MeV energies, we incorporated a blackbody component into the Band function. However, Blackbody component is not enough narrow to properly fit the narrow emission feature. Therefore, we introduced an additional Gaussian component to fit the spectrum of the time interval 8–30 s. The Gaussian function is defined as follows:

$${N}_{{{\rm{gauss}}}}(E)=A\frac{1}{{\sigma }_{{{\rm{gauss}}}}\sqrt{2\pi }}\exp \left(\frac{-{(E-{E}_{{{\rm{gauss}}}})}^{2}}{2\,{({\sigma }_{{{\rm{gauss}}}})}^{2}}\right).$$
(2)

where A is the normalization of spectrum, Egauss and σgauss are the central energy and standard deviation of the Gaussian function. We have found that the Gaussian component is well constrained at \({E}_{{{\rm{gauss}}}}=2154.6{0}_{-65.07}^{+53.37}\,{{\rm{keV}}}\), with a width \({\sigma }_{{{\rm{gauss}}}}=229.3{6}_{-45.29}^{+93.57}\,{{\rm{keV}}}\). The fitting results of the spectrum are presented in Table 1. The c and d panels of Fig. 2 displays the counts rate and νFν spectrum, with fitting using the Band function plus a Gaussian component. From the light curve presented in the (a) and (b) panel of Fig. 1, the MeV narrow emission feature appears during the rising and falling phases of the main pulse. When compared to other time intervals (0 − 5 s, 5 − 8 s, 30 − 36 s and 36 − 60 s), the time interval 8 − 30 s exhibits the highest flux and the best signal-to-noise ratio. The evolution of the spectral parameters of the Band function in the best-fit model is shown in the (c), (d), and (e) panels of Fig. 1. The low-energy spectral index α evolved from -0.88 to -1.33, indicating an evolution from hard to soft. Additionally, the peak energy Ep varies between 397 keV and 920 keV, showing the pattern of intensity tracking51.

For the A (0 − 5 s), B (5 − 8 s), D (30 − 36 s) and E (36 − 60 s) time intervals, we fixed the line width at σgauss = 200 keV and the line central energy at Egauss = 2.1 MeV in the likelihood fit, thereby deriving the upper limits on the flux of the narrow emission feature, which are Fluxgauss < 5.1 × 10−7 erg cm−2 s−1, Fluxgauss < 3.2 × 10−7 erg cm−2 s−1, Fluxgauss < 2.4 × 10−7 erg cm−2 s−1, and Fluxgauss < 9.9 × 10−8 erg cm−2 s−1, respectively.

Model comparison

We employed three different methods to assess the necessity of adding a Gaussian component to the prompt gamma-ray spectrum of GRB 221023A.

The Akaike Information Criterion (AIC) is employed for model comparison when penalizing additional free parameters is necessary to prevent overfitting. The AIC is formulated as the logarithm of the likelihood with a penalty term52,53:

$${{\rm{AIC}}}=-2{{\rm{ln }}}({{\mathcal{L}}}(d| \theta ))+2\theta .$$
(3)

where \({{\mathcal{L}}}(d| \theta )\) is the likelihood of the model, θ is the number of free parameters of a particular model. The model with the smallest AIC is favored. ΔAIC = AICBand − AICBand+Gaussian provides a numerical assessment of the evidence that model Band+Gaussian is to be preferred over model Band. When ΔAIC > 10, it strongly favors the model Band+Gassian. As shown in Table 2, our results reveal that during the time interval 8 − 30 s, the ΔAIC value reaches its maximum at 51.87, strongly favoring the Gaussian+Band model over the simpler Band model. In the four finer time-resolved spectra (8 − 21 s (C.1), 11 − 24 s (C.2), 14 − 27 s (C.3), 17 − 30 s (C.4)), the ΔAIC values vary between 25.76 and 36.55, further strongly favoring the addition of the Gaussian component.

When evaluating the significance of emission or absorption features in spectrum analysis, the Bayesian factor is also a commonly used tool20,54,55. The Bayesian factor is utilized to compare the relative support for different models, serving as a measure to evaluate the strength of evidence in favor of one model over another. Bayesian evidence (\({{\mathcal{Z}}}\)) is calculated for model selection and can be formulated as follows:

$${{\mathcal{Z}}}=\int{{\mathcal{L}}}(d| \theta )\pi (\theta )d\theta,$$
(4)

where π(θ) represents the prior distribution for θ. The ratio of the Bayesian evidence for two different models is called the Bayes factor (BF). In this paper, the BF is formulated as follows:

$${{\rm{BF}}}=\frac{{{{\mathcal{Z}}}}_{{{\rm{Band+Gaussian}}}}}{{{{\mathcal{Z}}}}_{{{\rm{Band}}}}},$$
(5)

The corresponding logarithmic expression is as follows:

$${{\rm{ln(BF)}}}={{\rm{ln }}}({{{\mathcal{Z}}}}_{{{\rm{Band+Gaussian}}}})-{{\rm{ln }}}({{{\mathcal{Z}}}}_{{{\rm{Band}}}}).$$
(6)

If ln(BF) >8, it indicates strong evidence in favor of the Band+Gaussian model56,57. We calculated the Bayes factors for time intervals with narrow emission features (as shown in Table 2), and the results shown that the Band+Gaussian model was preferred in finer time intervals (8 − 21 s, 11 − 24 s, 14 − 27 s, 17 − 30 s) with \({{\rm{\ln }}}({{\rm{BF}}})\) between 2.06 − 7.34. Remarkably, during the entire time interval of 8 − 30 s, the \({{\rm{\ln }}}({{\rm{BF}}})=9.99\) providing strong statistical support for the addition of the Gaussian component, suggesting the presence of the narrow emission feature.

We also employed the alternative analysis software GTBURST to extract the corresponding spectra from the time intervals exhibiting a narrow emission feature. The extracted spectra were fitted using the XSPEC 12.11.158, and the fitting results similarly indicate the presence of distinct narrow and bright emission feature between 1 MeV and 3 MeV. Δχ2 represents the statistical difference in the goodness-of-fit between the models Band and Band+Gaussian, the Δχ2 values are displayed in Table 2. The highest Δχ2 value of 40.14 was observed in the time interval 8 − 30 s, while the Δχ2 values for the other time intervals ranged from 18.53 to 34.49.

Background

The selection of time intervals for background subtraction can also impact the analysis of the source spectrum. In order to assess the impact of background subtraction on extracted spectrum.

In time interval 8 − 30 s, we calculated the background spectrum by selecting several different time windows. Even with this approach, the narrow emission features are still clearly visible. We performed both Band and Band + Gaussian fittings in the spectra extracted by subtracting different backgrounds in time interval 8 − 30 s. As shown in Table 3, The central energy Egauss of the narrow emission feature are all around 2.1 MeV and the widths σgauss are all around 200 keV, and the values of the ΔAIC are around 50. The result of the narrow Gaussian feature is substantially unaffected.

Table 3 Energy spectrum resulted from different time windows selected for the background region in time interval C (8 − 30 s)
Full size table

In four subintervals (8 – 21 s (C.1), (11 – 24 s (C.2), 14 – 27 s (C.3), 17 – 30 s (C.4)), we extracted the spectra by performing a different selection of the time windows for the background spectrum computation, the background time intervals selected for each time intervals, for 8 – 21 s:  −200 – −40 s, 120 – 250 s; for 11 – 24 s:  −90 – −10 s, 100 – 150 s; for 14 – 27 s:  −90 – −20 s, 180 – 250 s; for 17 – 30 s:  −200 – −50 s, 120 –250 s.

Significance calculation of narrow emission feature

We calculated the chance probability value (p-value) of the narrow emission feature through spectral simulation. The spectral simulation across the entire energy range (10 keV − 30 MeV) is performed using the fakeit command in XSPEC. These simulations are based on the parameters obtained from fitting the actual data using the Band model. The tclout simpars (based on the covariance matrix at the best fit) command in XSPEC is used to generated randomized model parameters before each simulation. The total number of spectral simulations N is 1.00 × 107. For each simulated spectra, we perform both Band and Band+Gaussian fittings (search for Gaussian components across the entire energy range of 10 keV to 30 MeV) and record the maximum Δχ2 value20,55. Finally, we assess the significance of the narrow emission feature by analyzing the Δχ2 values recorded in Tables 2. The p-value represents the fraction of simulated \(\Delta {\chi }_{i}^{2}\) values that exceeds the observed Δχ2 value:

$$p\mbox{-}{{{\rm{value}}}}_{{{\rm{sim}}}}=n[\Delta {\chi }_{i}^{2}\ge \Delta {\chi }^{2}]/N.$$
(7)

If after N simulations we still do not obtain a Δχ2 value exceeding the observed value, we report p-\({{{\rm{value}}}}_{{{\rm{sim}}}}\, < 1/N\). The probability distribution function (PDF) of Δχ2 values obtained from 1 × 107 simulations for different time intervals are shown in Fig. 5.

Fig. 5: Simulated Δχ2 distribution.
figure 5

Δχ2 is the statistical difference in the goodness-of-fit between the models Band and Band+Gaussian. The Panels a, b, c, d, and e show the probability distribution function (PDF) of Δχ2 values obtained from 1 × 107 simulations for time intervals 8 − 30 s, 8 − 21 s, 11 − 24 s, 14 − 27 s, 17 − 30 s, respectively. The red dashed line represents the observed Δχ2 value. The p-\({{{\rm{value}}}}_{{{\rm{sim}}}}\) is the chance probability value obtained from 1.00 × 107 simulations. If after N simulations no Δχ2 value exceeds the actual fitting result, we report p-\({{{\rm{value}}}}_{{{\rm{sim}}}}\, < 1/N\). Source data are provided as a Source Data file.

Full size image

In the process of calibrating the Δχ2 test distribution through simulation, the intensity, ___location and width of the line, are not fixed to predetermined values but are allowed to vary freely during the fit. This is a standard setup when performing the simulation. The number of independent search trials conducted by dividing multiple time intervals in the time series of different GRBs must be considered (the look-elsewhere effect59). The chance probability value p-valuesim-trial after considering the correction for the number of independent search trials on the basis of the p-\({{{\rm{value}}}}_{{{\rm{sim}}}}\) is59,60,61:

$$p{\mbox{-}}{{{\rm{value}}}}_{{{{\rm{sim}}}{\mbox{-}}} {{\rm{trial}}}}=1-{(1-p {\mbox{-}}{{{\rm{value}}}}_{{{\rm{sim}}}})}^{t}.$$
(8)

where t is the number of independent search trials. We searched for GRBs spectral lines from the Fermi-GBM catalog (https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermigbrst.html) in descending order of fluence29. We excluded GRB 221009A, which already has identified narrow emission features22,23. The extreme brightness of GRB 130427A and GRB 230307A caused detector pile-up effects, so we excluded the saturated time intervals of 4.5 − 11.5 s for GRB 130427A and 3 − 7 s for GRB 230307A62,63. We searched a total of 9 GRBs, for each burst, time intervals were divided based on BGO light curve signal-to-noise ratio greater than 40. This resulted in a total of 256 searches. Therefore, the number of independent search trials t = 256.

We found the highest statistical significance of narrow emission feature in the time interval 8 − 30 s, with the chance probability value \(p\mbox{-}{{{\rm{value}}}}_{{{\rm{sim}}}} < \;1\times 1{0}^{-7}\) obtained from results of 1.00 × 107 simulations, corresponding to the Gaussian-equivalent significance  >5.32σ. Considering the correction for the number of independent search trials, the chance probability value decreases to p-valuesim-trial < 2.56 × 10−5, corresponding to the Gaussian-equivalent significance  >4.20σ. The chance probability values for the other time intervals are shown in Table 2.