GRB 221009A: Discovery of an Exceptionally Rare Nearby and Energetic Gamma-Ray Burst

Figure 1 shows the BAT raw light curves (i.e., not background subtracted), summed over all detectors, from T₀ − 500 s to T₀ + 5000 s. The location of GRB 221009A was occulted by the Earth until T₀ + 1870 s. At ∼T₀ + 1100 s, the overall count rate began to rise due to increased particle background as Swift approached the South Atlantic Anomaly (SAA). From T₀ + 1317–2183 s, BAT data collection was disabled as Swift transited the SAA. The count rate remained elevated while exiting the SAA until the spacecraft slew beginning at T₀ + 2496 s, when we attribute the linearly decaying enhanced count rate to emission from GRB 221009A. However, the source location was not in the coded portion of the BAT field of view until a slew beginning at T₀ + 3095 s.

Zoom In Zoom Out Reset image size

Finally after this slew completed, BAT triggered on GRB 221009A, at 14:10:18 (T₀ + 3199 s; trigger ID 1126853) and 14:17:06 (T₀ + 3607 s; trigger ID 1126854) UTC. The event data from these two triggers cover a time range from T₀ + 2960 s to T₀ + 4570 s. The mask-weighted light curve shows steadily declining emission present when the burst location came into the BAT field of view at T₀ + 3173 s, and extending beyond the available event data range. The time-averaged spectrum from T₀ + 3302 s to T₀ + 4538 s is best fit by a simple power-law model, with Γ = 2.08 ± 0.03. The fluence in the 15–150 keV band is (7.4 ± 0.1) × 10⁻⁵ erg cm⁻².

The smooth temporal evolution observed by BAT, together with the large temporal offset between the BAT and GBM triggers, indicates that BAT triggered on the afterglow of GRB 221009A. This marks the first such occurrence of a BAT afterglow trigger in the 18 yr of Swift operations.

Given the exceptionally bright afterglow, we searched for even later emission in the BAT survey mode data using the BatAnalysis Python package (T. Parsotan et al. 2023, in preparation). In individual pointings of survey mode data, the afterglow was detected until 2022 October 9 21:55:38 UT (T₀ + 31 ks). We attempted to fit the spectra of each survey data set with a power-law (cflux*po) model in XSPEC (Arnaud 1996) to obtain fluxes and photon indices. Nondetections were then analyzed to obtain 5σ upper limits following the procedure outlined in Laha et al. (2022).

At later times, the survey data were binned daily and then mosaiced together. The spectrum from each mosaiced image was fitted, and the flux and photon indices were derived similar to the procedure above. The results of this analysis are plotted in the top panel of Figure 2, and listed in Appendix A.

Zoom In Zoom Out Reset image size

The XRT began observing GRB 221009A at 14:13:09 UT (T₀ + 3370 s, 170 s after the first BAT trigger) and located a bright afterglow in the initial 0.1 s image-mode exposure, after which observations began in windowed timing (WT) mode. The initial WT count rate was 910 ± 40 ct s⁻¹ (all XRT count rates are corrected for the effects of pile-up and hot columns; see, e.g., Evans et al. 2007); making it more than an order of magnitude brighter at this time—in observed flux—than any other GRB observed by XRT. Due to the high count rate, the XRT remained initially in WT mode, in which only 1D spatial information is collected. Significant structures are present in the 1D spatial profiles, with a clear excess compared to the expected point-spread function (PSF), which were evolving with time (see Appendix B). This resembles the behavior expected when dust clouds in our Galaxy scatter X-rays from the GRB prompt emission, which were not initially traveling toward Earth, back into our line of sight.

At T₀ + 89 ks, the GRB had faded sufficiently that the XRT automatically switched to photon counting (PC) mode. The 2D image from this observation, shown in Figure 3, confirmed the presence of a complex series of expanding, bright rings associated with a dust-scattering echo (see also Tiengo et al. 2022; Vasilopoulos et al. 2023), the properties of which are discussed in Section 3.1.

Zoom In Zoom Out Reset image size

The presence of scattered emission complicates the data analysis, as it can contribute events to the regions over which source and background counts are accumulated. Both the intensity and spectrum of the rings are spatially variable; thus the selected background region may in fact not be representative of the background within the source region. As a result, the automated XRT analysis (Evans et al. 2007, 2009) needs some modification. For PC mode, this is relatively simple, as we have full 2D imaging, and the innermost rings were reasonably well separated from the GRB itself. We restricted the source region to a radius of 20 pixels (47''), and extracted a combined background from four regions ³⁷  across the detector, which are as free from dust contamination as possible (identified from a stacked image of all PC-mode data). In order to quantify the impact of dust contamination in the background, we calculated the background rate in these regions, and found an average value of 1.6 × 10⁻⁶ ct s⁻¹, which is consistent with historical values (Evans et al. 2023). During the first ∼5 days, dust ring contamination caused a 4–10 times elevation of the background rate. However, during this period, the afterglow of GRB 221009A is 30–50 times brighter than the background level, meaning the effect is negligible.

For WT mode, this is more problematic since the CCD is read out in columns, rather than pixels. After investigating the variation in the echo contribution to WT-mode data, we modified the default WT extraction regions to minimize its impact (see Appendix B). We found that the WT flux is accurate to ∼6%. We further verified this by extracting a light curve using only data above 4 keV, where dust scattering becomes less efficient, given the roughly ν⁻² dependence of the scattering cross section; only minimal changes in the light-curve shape are observed. On the other hand, the dust has more significant impact on the WT spectra (Appendix B), resulting in an increase in the best-fit photon index up to ∼7% and a significant increase to the best-fit absorption column (up to ∼27%).

We modified the settings for the automated light curve ³⁸  to limit the source extraction region to 20 pixels radius, and to use the PC-mode background region (described in the Appendix B); due to the unusual brightness of the source, we also increased the number of counts per bin by a factor of 10 compared to the default (see Evans et al. 2007 for details). We fit the time-averaged PC-mode spectrum using XSPEC with a power-law model and two absorption components. The first of these was a tbabs fixed to the Galactic value of 5.38 × 10²¹ cm⁻² (Willingale et al. 2013); the second was a ztbabs with N_H free and redshift fixed at 0.151. ³⁹  The fit yields N_H = (1.4 ± 0.4) × 10²² cm⁻², with a photon index Γ = 1.8 ± 0.2. We used these results to convert the light curve into observed 0.3–10 keV flux. The resulting flux measurements are shown in Figure 2.

To investigate possible spectral evolution, we extracted a series of time-resolved spectra using the "Add time-sliced spectrum" option on the UKSSDC website ⁴⁰   (Evans et al. 2009), which ensures that the modifications to the default processing parameters, described above, are employed. From the WT-mode data, we created one spectrum from the first spacecraft orbit (T₀ + 3.4–4.5 ks), and then spectra in 10 ks chunks from T₀ + 10–50 ks. For PC mode, we extracted spectra over the intervals T₀ + 67–200, 200–400, 400–1000, and >1000 ks. For WT mode, we used only grade 0 events, as is recommended for very absorbed objects. ⁴¹  We fitted these spectra independently with the model defined above; the results are shown in Table 1. The photon index ranges between ∼1.6 and ∼1.9 with the harder profiles observed during the first block of observations (from 3.4 to 4.5 ks) and around 1–2 days (from 68 to 175 ks) after T₀.

Table 1. Spectral Parameters Determined from Fitting a Power-law Model to the Time-resolved X-Ray Data from Swift, MAXI, and NICER

T − T0	Exposure	Intrinsic NH	Photon Index	Fit Stat./d.o.f.	Observatory/Detector
(ks)	(ks)	(1022 cm−2)
2.5	0.047	⋯	1.75 ± 0.09	307/322	MAXI/GSC
3.4–4.5	1.2	1.28 ± 0.04	1.61 ± 0.02	1249/959	XRT/WT
8.0	0.047	⋯	${2.07}_{-0.26}^{+0.28}$	145/166	MAXI/GSC
10.0–20.0	0.5	1.39 ± 0.09	1.89 ± 0.05	859/840	XRT/WT
14.0	0.029	⋯	1.89 ± 0.16	26.9/35	NICER/XTI
19.5	0.062	⋯	1.84 ± 0.15	23.6/36	NICER/XTI
20.0–30.0	2.9	1.36 ± 0.05	${1.84}_{-0.04}^{+0.02}$	995/922	XRT/WT
30.0–40.0	1.7	1.25 ± 0.07	1.85 ± 0.04	800/840	XRT/WT
40.0–50.0	2.8	1.20 ± 0.06	1.84 ± 0.03	875/873	XRT/WT
52.6–98.0	2.1	⋯	1.85 ± 0.07	131/153	NICER/XTI
68–175	17	${1.27}_{-0.15}^{+0.16}$	1.65 ± 0.08	655/727	XRT/PC
103.0–198.0	6.5	⋯	1.79 ± + 0.06	144/168	NICER/XTI
203.3–360.5	9.8	⋯	1.61${}_{-0.09}^{+0.10}$	94/94	NICER/XTI
216–394	12	${1.26}_{-0.17}^{+0.18}$	1.80 ± 0.10	632/665	XRT/PC
404–954	26	1.20 ± 0.18	1.92 ± 0.10	501/622	XRT/PC
505.0–961.6	2.7	⋯	1.82 ± 0.19	58.1/67	NICER/XTI
1005–6285	0.24	1.00 ± 0.15	1.91 ± 0.09	723/668	XRT/PC

Note. NICER/XTI and MAXI/GSC data are fit in the energy ranges 4–10 and 4–20 keV, respectively, to avoid contamination from unresolved dust echo contamination, with intrinsic N_H fixed to 1.29 × 10²² cm⁻². XRT spectra are fit in the range 0.3–10 keV, allowing for measurement of intrinsic N_H. Note that fit statistics are C-stat (Cash 1979) for MAXI and Swift data; NICER are χ². Galactic N_H is fixed to a value of 5.38 × 10²¹ cm⁻² (Willingale et al. 2013).

Download table as:  ASCIITypeset image

The UVOT began settled observations of the field of GRB 221009A at 2022 October 9 14:13:17 UT, T₀ + 3.4 ks (179 s after the first BAT trigger; Dichiara et al. 2022; Kuin et al. 2022). A nearby star, 5'' away from the GRB, contaminates the photometry when using the standard circular aperture of 5'' radius. In order to minimize the star's contribution to the measurement of the afterglow, we use a 25 aperture to extract source counts. To be consistent with the UVOT calibration, these count rates were then corrected to 5'' using the curve of growth contained in the calibration files. Background counts were extracted using two circular regions of radius 15'' located in source-free regions identified from stacked deep images in each filter. The count rates were obtained from the image lists using the Swift tool uvotsource.

As GRB 221009A fades, even using a small aperture, the contribution from the nearby star increases in dominance. In order to estimate the level of contamination, for each filter, we combined the late-time exposures between T₀ + (3.4 × 10⁶) s and T₀ + (4.4 × 10⁶) s. We extracted the count rate in the late combined exposures using the same 25 aperture and applied an aperture correction to 5''. These were subtracted from the source count rates to obtain the afterglow count rates. The afterglow count rates were converted to magnitudes using the UVOT photometric zero-points (Poole et al. 2008; Breeveld et al. 2011), as well as to fluxes. To improve the signal-to-noise ratio (S/N), the count rates in each filter were binned using Δt/t = 0.2. The final photometry is plotted in the bottom panel of Figure 2 and listed in Appendix C. We caution that the photometry in the uvw1 and uvw2 is likely impacted by red leak; see Siegel et al. (2014), Brown et al. (2010) for additional details.

MAXI, which is on board the International Space Station (ISS), performs scanning observations of about 80% of the sky every 92 minutes. Prior to the GBM detection of the prompt emission from GRB 221009A, the two gas slit camera (GSC; Mihara et al. 2011) units, GSC_4 and GSC_5, scanned the source region at 12:25 UT on 2022 October 9 (T₀−3.1 ks), and no enhancement was seen at that time. The 2–10 keV 3σ upper limit is 0.074 ct cm⁻² s⁻¹, corresponding to a flux of approximately 7.9 × 10⁻¹⁰ erg cm⁻² s⁻¹. At 13:58 UT (T₀+2.5 ks) the GSCs detected a very bright transient with a 2–10 keV flux of 5.97 ± 0.19 ct cm⁻² s⁻¹ (∼6.5 × 10⁻⁸ erg cm⁻² s⁻¹; Negoro et al. 2022). Unfortunately, at the time the data were collected, there was a loss of the signal to the ground, and the data were stored in the High Rate Communication Outage Recorder in the ISS, leading to a delay in reporting the transient. Before downloading these stored data at around 15:55 UT, the MAXI/GSC Nova-Alert System (Negoro et al. 2016) triggered on the source at 15:31:03 UT (T₀+8.0 ks), about 81 minutes after the BAT trigger. In this second observation by MAXI, the 2–10 keV flux had dropped to 1.00 ± 0.08 ct cm⁻² s⁻¹ (about 1 × 10⁻⁸ erg cm⁻² s⁻¹).

Thanks to the exceptional brightness, the GSCs significantly detected the afterglow during four scans that followed the second detection, to 21:41 UT, six scans in total (e.g., Kobayashi et al. 2022). Typically, GRBs are not detected beyond the first one or two scans, taken at ∼92 minute intervals. After the six detections, the sky region containing the source was not observable by MAXI until 15:31 UT on 2022 October 11.

We performed a spectral analysis using data obtained from GSC_4 and GSC_5 (for details see Appendix D). We first attempted a standard analysis using the same model applied to the XRT spectra, i.e., a power law with two absorption components, with a fixed Galactic absorption column density of N_H = 5.38 × 10²¹ cm⁻². However, we obtained steeper photon indices and/or lower intrinsic absorption column densities than those of the XRT (see Table 7). If we fixed N_H,zabs at the weighted mean value obtained with the XRT at T₀ + (3.4 − 50) ks, 1.29 × 10²² cm⁻², even steeper photon indices were obtained.

These discrepancies can be explained by the GSC spectra containing dust scattered soft X-rays, which cannot be spatially resolved by the GSC. Mitigation of the dust scattering process is discussed in Appendix D; to summarize, an additional spectral component to account for the scattered emission was included with a power-law index derived from XRT PC-mode data and the associated flux allowed to vary as a free parameter. The afterglow photon index was fit as a free parameter in the data from the first two scans, but fixed at the weighted mean XRT value at T₀ + (10 − 50) ks due to poor statistics. The resulting spectral parameters are displayed in Table 1, and the flux measurements (extrapolated to the 0.3–10.0 keV band) are plotted in Figure 2.

NICER received notification of this GRB from MAXI on the ISS via On-orbit Hookup of MAXI and NICER at 14:10:57 UT, but could not begin regular observations due to poor visibility until T₀ + 52.870 ks. Only two observations were possible during this early phase. The first and second observations were at T₀ + 14.003 ks and T₀ + 19.560 ks with an exposure time of 29 and 62 s, respectively. After T₀ + 52.870 ks, NICER intermittently made observations over a period of about 42 days.

NICER's X-ray Timing Instrument (XTI) consists of 56 modules, each comprising an X-ray concentrator optic and silicon drift detector (Gendreau et al. 2012). During these observations, 52 XTI modules were operational. We reduced the data with the nicerl2 command and generated level (2) cleaned events. To avoid potential contamination of instrumental noise, we further excluded XTI module numbers (14) and (34). We defined a time interval in which an exposure time of more than 200 s can be continuously secured as a one-time interval and created an energy spectrum at each time interval (except for the first two observations) with good statistics.

To select a time period with less background influence due to charged particles, we also checked the distribution of the rate of 12–15 keV photons, where the optics have less effective area. Most data are distributed around 0.048 ct s⁻¹ and can be represented by a normal distribution with a standard deviation of 0.018 ct s⁻¹. Therefore, we removed data with a rate greater than 2σ from the mean in the 12–15 keV range after T₀ + 100 ks seconds where the background contribution becomes large (>0.084 ct s⁻¹). We also did not use data after T₀ + 1000 ks because the background rate increased due to a geomagnetic storm. The background spectral model is estimated by the 3C 50 background model (Remillard et al. 2022) using the nibackgen3c50 command.

To explore the spectral evolution as in Section 2.2, we divided the NICER data into six time intervals and fitted the time-resolved spectra with a power-law model and the same two absorption components. Since NICER is a nonimaging detector, it cannot remove photons originating from the dust echo within its 5' diameter field of view. To evaluate the impact of the dust echo, we fitted the data in two energy ranges, 1–8 and 4–8 keV bands. Since the data fitting in the 4–8 keV range does not allow us to constrain the intrinsic N_H value, we fixed it to 1.29 × 10²² cm⁻², which is the weighted mean value obtained by XRT at T₀ + (10.0–961.6) ks (Table 1). Comparing the fitting results in the two energy bands, the results for the 4–8 keV band show that the photon indices are systematically harder than those for the 1–8 keV band results and consistent with the XRT results. Therefore, we conclude that 4–8 keV light curve is much less affected by the dust echo, and report these results in Table 1.

While complicating the point-source analysis of the afterglow emission, the dust scattering echo contains a rich amount of information about the intervening Galactic dust along the line of sight to GRB 221009A. Dust echoes from GRBs offer the unique opportunity to measure distances to dust clouds geometrically, with a precision limited only by the angular resolution of the telescope. This is because the distance D_GRB to the X-ray source is effectively infinite, relative to the distance to the dust. This allows us to eliminate the source distance from the scattering angle equation (e.g., Heinz et al. 2015) and solve it directly for the distance D_dust to the dust (using the small angle approximation):

$$\begin{eqnarray}&&{D}_{\mathrm{dust}}=\displaystyle \frac{{D}_{\mathrm{GRB}}}{1+\tfrac{{D}_{\mathrm{GRB}}{\theta }^{2}}{2c{\rm{\Delta }}t}}\approx \displaystyle \frac{2c{\rm{\Delta }}t}{{\theta }^{2}}\end{eqnarray} \tag{ 1 }$$

where Δt is the time delay between the GRB and the exposure, and θ is the off-axis angle of the ring relative to the position of the GRB.

Because the vast majority of the fluence of the GRB was concentrated in the prompt GRB emission, to the lowest order, we can approximate the input light curve that caused the echo as a delta function in time. This implies that there is a single, unique Δt for each photon detected by the XRT, and thus, each photon can be referenced to a specific dust distance in Equation (1). This makes the analysis of this echo substantially simpler than in comparable echoes from Galactic X-ray transients.

Details of our analysis of the dust echoes from GRB 221009A are provided in Appendix E. Figure 4 shows the best-fit column density histogram for a standard Mathis–Rumpl–Nordsieck (MRN) model (Mathis et al. 1977), assuming a soft X-ray fluence of ${{ \mathcal F }}_{0.8-5\,\mathrm{keV}}=2.1\times {10}^{-3}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}$ (S. Lesage et al. 2023, in preparation). The figure shows clear evidence for several dust components on the near side of the Galaxy, with the largest column densities found within distances of about 1 kpc, as expected given the Galactic latitude of GRB 221009A.

Zoom In Zoom Out Reset image size

However, there is also clear evidence for dust at larger distances, located well above the Galactic disk (see also Negro et al. 2023). While the column densities of these dust concentrations are small, the large cross section at the correspondingly small scattering angles of this more distant dust makes the scattering from these clouds detectable in the radial intensity profiles, as can be seen in Figure 13, for example. As such, the dust echo from GRB 221009A proves that echo tomography is exquisitely sensitive to small dust concentrations at large distances that are not typically measurable with other techniques.

The largest column density cloud is located approximately between 400 and 600 pc distance, with a very nearby dust cloud at a distance of about 200 pc and a third cloud at distance approximately 700 pc. Naturally, these structure are very local compared to the large scale Galactic structure. The very rapid detection of the echo by the XRT allows us to map very nearby dust typically hard to study using dust echoes.

We find a total Galactic column density of N_H ∼ 5.9 × 10²¹ cm⁻², broadly consistent with the Galactic column density of 5.38 × 10²¹ cm⁻² found by Willingale et al. (2013). Given that the largest uncertainty in this value derives from the poorly constrained estimate of the soft X-ray fluence of GRB 221009A, we can infer that our assumed fluence of ${{ \mathcal F }}_{0.8-5\,\mathrm{keV}}=2.1\times {10}^{-3}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}$ is likely correct within roughly 10%.

We find that a simple MRN distribution provides an excellent fit to the temporal evolution of the radial intensity profile, without the need for more complex dust chemistry (such as ice mantles). In this sense, the dust toward GRB 221009A appears to be typical for interstellar dust seen along sightlines to other recent X-ray light echoes, such as the echo from V404 Cyg (Beardmore et al. 2016; Heinz et al. 2016; Vasilopoulos & Petropoulou 2016).

The combined XRT+MAXI+NICER data set allows us to explore the soft X-ray afterglow evolution in great detail. In the earliest soft X-ray spectra (T₀ + 2.5–10 ks), the afterglow exhibits strong spectral evolution, with the XRT WT-mode power-law photon index increasing (i.e., getting softer) by ∼17% (Table 1). This change is much larger than the inferred WT-mode systematic uncertainty (Appendix B), and is corroborated by the 2–20 keV MAXI spectra derived around this time (Appendix D).

From T₀ + 10–50 ks, the spectra show no definitive signs of evolution; there are hints at a reduction in the host galaxy absorption throughout the XRT WT-mode data, although this trend is reversed when XRT switched to PC mode. Given the differences between the WT- and PC-mode results are much less than the dust-induced systematics on the WT results (Appendix B), and no similar spectral evolution is observed in the NICER spectra, we do not believe there is supportable evidence for spectral evolution at these times.

On the longer timescales probed by the XRT PC-mode data, clear spectral evolution is seen and takes the form of an increasing (i.e., softer) photon index, which is confirmed by the 4–8 keV NICER spectra (Table 1).

To constrain the broadband spectral behavior, we built and fit two spectral energy distributions (SEDs) at T₀ + 4.2 ks and T₀ + 43 ks, following the procedure outlined in Schady et al. (2007, 2010), using data from the BAT survey mode, XRT, and UVOT. For the 4.2 ks SED, we use data in the range 2.5–4.7 ks. For the 43 ks SED, we used data in the range 23–103 ks. Details of the SED construction are provided in Appendix F.

The SEDs were fit using XSPEC. We tested two different models for the continuum: a single power law, and a broken power law with the change in spectral slope fixed to be Δβ = 0.5 (corresponding to the expected change in spectral slope caused by the synchrotron cooling frequency; Sari et al. 1998). In each of these models, we also included two dust and gas multiplicative components to account for the Milky Way and host galaxy dust extinction and photoelectric absorption (tbabs, ztbabs, zdust). ⁴²  We also include a component for attenuation by the intergalactic medium (zigm). The Galactic components were frozen to the reddening and column density values from Schlafly & Finkbeiner (2011), Willingale et al. (2013), respectively. In the analysis, we found consistent results between Milky Way, Small Magellanic Cloud (SMC), and Large Magellanic Cloud dust extinction laws; therefore we only report the fits using the SMC dust extinction law.

Due to the uncertain contribution to the XRT WT spectra by the dust rings, we extracted WT spectra using the 5–20 pixel annular extraction region identified in Appendix B for the spectrum for the early SED, and a 10 pixel radius circular extraction region for the later time SED, chosen to minimize the fraction of dust scattered emission in the extraction region, with a background spectrum taken from a 75–125 pixel annular region. ⁴³  Since the intrinsic N_H in the WT fits is sensitive to the dust contamination (see Appendix B and Table 5), we also provide fits using BAT and UVOT alone. The results of the analysis are given in Table 2 and plotted in Figure 5.

Zoom In Zoom Out Reset image size

Table 2. Fits to the Two Spectral Energy Distributions, at T₀ + 4.2 ks and T₀ + 43 ks, Built Using UV–Optical (U), X-Ray (X), and Gamma-Ray (B) Data

SED	Model	NH,X,host	E(B − V)host	Γ	Ebk	χ2 (d.o.f)	Null Hypothesis Probability
		(1022 cm−2)	(mag)		(keV)
B+U 4.2 ks	POW	⋯	0.26 ± 0.05	1.74 ± 0.02	⋯	16 (9)	7.3 × 10−2
B+U 4.2 ks	BKPOW	⋯	0.23 ± 0.05	${1.70}_{-0.12}^{+0.03}$	${33.6}_{-28.6}^{+16.1}$	4 (8)	8.4 × 10−1
B+X+U 4.2 ks	POW	1.75 ± 0.03	1.04 ± 0.02	1.93 ± 0.01	⋯	1372 (681)	6.9 × 10−49
B+X+U 4.2 ks	BKPOW	1.35 ± 0.03	0.51 ± 0.03	1.69 ± 0.01	6.8 ± 0.3	861 (680)	2.7 × 10−6
B+U 43 ks	POW	⋯	0.38 ± 0.06	1.89 ± 0.04	⋯	14 (9)	1.2 × 10−1
B+U 43 ks	BKPOW	⋯	0.38 ± 0.06	1.39 ± 0.04	${7.3}_{-7.3}^{+17}\times {10}^{-4}$	14 (8)	8.3 × 10−2
B+X+U 43 ks	POW	1.27 ± 0.03	0.49 ± 0.03	1.83 ± 0.01	⋯	778 (678)	4.2 × 10−3
B+X+U 43 ks	BKPOW	1.14 ± 0.08	0.28 ± 0.03	1.72 ± 0.02	4.5 ± 0.2	748 (677)	3.0 × 10−2

Note. Both epochs were fit with a power-law (POW) continuum and a broken power-law (BKPOW) continuum, accounting for Galactic and host galaxy gas and dust. The columns are time of the SED; model; host galaxy equivalent column density N_H,X,host; host galaxy reddening, E(B − V)_host; power-law photon index or the first photon index of the broken power-law model Γ (Γ₂ is fixed to be Γ + 0.5); break energy for the broken power-law spectral models E_bk; χ², degrees of freedom (d.o.f); and the null hypothesis probability.

Download table as:  ASCIITypeset image

Examining the SED at T₀ + 4.2 ks including BAT, XRT, and UVOT data, we find that a broken power law is preferred, with the F-test finding a statistical improvement of ≫5σ. The break energy required is at ∼7 keV. At the same epoch, excluding the XRT data, we find again a broken power-law model is preferred, with the F-test suggesting the improvement is >3σ. The photon index for the broken power-law model is consistent with the fits including the XRT spectral files. However, the intrinsic dust extinction is lower: 0.23 ± 0.05 mag compared to 0.51 ± 0.03 mag including all data. Not surprisingly, excluding the XRT data greatly increases the uncertainty in the break frequency.

For the SED at T₀ + 43 ks, a broken power law is again preferred, with the F-test suggesting a break is statistically required at >3σ. The break frequency remains in the X-ray bandpass, with evidence for a decline in energy over this period. The fits still prefer a significant host galaxy reddening, but the exact value of E(B − V)_host again varies depending on the inclusion or exclusion of the XRT data.

To summarize, both at T₀ + 4.2 ks and T₀ + 43 ks, we find strong evidence favoring a broken power-law model, with the lower energy photon index Γ ≈ 1.7. Both epochs favor a break energy in the X-ray bandpass, and a modest amount of intrinsic absorption and/or reddening in the host galaxy: E(B − V)_host ≈ 0.3–0.5 mag; N_H ≈ (1.1–1.4) × 10²² cm⁻². These results are broadly consistent with those reported in Levan et al. (2023), Kann et al. (2023). More precise constraints, however, are precluded due to the complications resulting from the dust scattering echoes.

Finally, we note that, while we have assumed here that any spectral break must be due to the synchrotron cooling frequency (and thus Δβ = 0.5), alternative explanations for complex spectral behavior aside from standard forward shock afterglow evolution have been suggested previously in the literature. For example, the possibility of dust destruction by the intense early-time UV flash could lead to spectral evolution due to a time-variable column density (Morgan et al. 2014). Alternatively, thermal emission has been claimed in the X-ray spectra of a number of previous GRB afterglows (e.g., Campana et al. 2006). Given the complexities induced by the dust echoes, however, a detailed investigation of the X-ray spectral evolution is beyond the scope of this work.

The combined X-ray and UV–optical light curves of GRB 221009A are shown in Figure 2 (top and bottom panels, respectively). The XRT, MAXI, and NICER measurements are placed in a common bandpass from 0.3–10.0 keV. The BAT survey mode data, however, are not extrapolated to this band and instead plotted from 14 to 195 keV—due to the inference of a spectral break below this band (Section 3.2), these points are instead meant to be illustrative.

While soft X-ray observations of the afterglow did not begin until T₀ + 2.5 ks, the light curve clearly lacks signatures of the canonical X-ray afterglow behavior at early times, including the prompt decay phase, the plateau phase, and any prominent flaring (Nousek et al. 2006). This relatively simple monotonic decay has been noted previously in other GRBs detected at high (≥GeV) energies (Yamazaki et al. 2020).

We fit the joint 0.3–10.0 keV (XRT+MAXI+NICER) X-ray light curve using the Markov Chain Monte Carlo software package emcee (Foreman-Mackey et al. 2019). An additional fixed fractional error was included as a nuisance parameter to account for cross-calibration uncertainty. We find a broken power-law model is strongly preferred over a single power law, although still formally a more complex behavior is indicated. The best-fit parameters for the broken power-law model are α_X,1 = −1.498 ± 0.004, α_X,2 = −1.672 ± 0.008, and ${t}_{\mathrm{break},{\rm{X}}}=({7.9}_{-1.0}^{+1.1})\times {10}^{4}$ s (68% confidence intervals). The initial MAXI detection at T₀ + 2.5 ks clearly indicates a shallower decay at early times (T₀ + ≲ 3.3 ks), while a late-time excess (T₀ + 2 Ms) may indicate some energy injection (e.g., Zhang et al. 2006), or alternatively could be due to contamination from, e.g., an unrelated foreground object.

For the UV–optical data, we constructed a combined single-band light curve following the procedure outlined in Oates et al. (2009). We then fit this light curve using the identical procedure described above for the X-ray data (including the cross-calibration nuisance parameter). We find that a single power-law model and broken power-law model provide comparable quality fits to the data—the statistically preferred model depends sensitively on how late-time (>10⁵ ks) measurements are included in the fit. The best-fit index derived for the single power-law model is α_O = −1.13 ± 0.01, while for the broken power-law model the parameters are ${\alpha }_{{\rm{O}},1}=-{0.98}_{-0.05}^{+0.11}$, ${\alpha }_{{\rm{O}},2}=-{1.31}_{-0.07}^{+0.05}$, and ${t}_{\mathrm{break},{\rm{O}}}=({2.2}_{-1.1}^{+1.7})\times {10}^{4}$ s. Regardless of the model selected, the UV–optical light curve clearly declines more slowly than the X-ray emission. Given the large foreground extinction, these wavelengths are unaffected by contamination from an emerging supernova or an underlying host galaxy, and thus this accurately reflects the afterglow behavior.

The top panel of Figure 6 plots the observed 0.3–10.0 keV flux light curve of GRB 221009A, overplotted on the entire sample of XRT light curves since the Swift launch (all uncorrected for absorption). At the time of the first XRT observations, GRB 221009A is approximately an order of magnitude brighter than any previous X-ray afterglow. Put differently, forming the distribution of X-ray afterglow flux at T₀ + 4.5 ks (a time selected to minimize the number of events requiring interpolation due to orbital gaps in light-curve coverage), GRB 221009A falls 5.5σ above the mean.

Zoom In Zoom Out Reset image size

However, GRB 221009A is also one of the nearest GRBs detected since the launch of Swift—12 events (out of ≈400) have lower redshifts reported in the online Swift GRB Table. ⁴⁴  To compare to the intrinsic properties of the X-ray afterglow sample, we converted the observed count-rate light curve into intrinsic luminosity in the 0.3–10 keV energy band in the GRB comoving frame. The luminosity in a given bin is given by the following:

$$\begin{eqnarray}&&L=4\pi {D}_{L}^{2}{{RC}}_{u}k\end{eqnarray} \tag{ 2 }$$

where D_L is the luminosity distance (749.3 Mpc), R is the measured 0.3–10 keV count rate (in the observer frame), C_u is the conversion from count rate to unabsorbed 0.3–10 keV flux (also in the observer frame) obtained from the spectral fits above (Section 2.2), and k is the k-correction of Bloom et al. (2001), which corrects from observed 0.3–10 keV flux to the 0.3–10 keV flux in the GRB comoving frame. To calculate k, we assume the spectrum is an unbroken power law with the photon index Γ = 1.8. The time axis must also be corrected to the observer frame, by dividing the time by 1 + z = 1.151. The resultant light curve is shown in the bottom panel of Figure 6.

While the lack of early-time observations somewhat complicates the comparison, it is clear that GRB 221009A is at the high end of the X-ray afterglow luminosity distribution (1.8σ above the mean). But unlike when considering this comparison in flux space, GRB 221009A is by no means an outlier. Rather we conclude that the exceptional observed brightness of the X-ray afterglow of GRB 221009A results from the combination of being very nearby and very intrinsically luminous: while neither the distance nor the luminosity are unprecedented, together they make GRB 221009A unique among the Swift afterglow sample.

Comparing GRB 221009A to the broader population of UV–optical afterglows is further complicated by the large (and uncertain) line-of-sight extinction, both foreground and in the host galaxy. With ${m}_{U}={17.59}_{-0.12}^{+0.13}$ mag at T₀ + 3574 s, GRB 221009A is significantly fainter than the brightest UVOT-detected afterglows, such as GRB 080319B (m_U ≈ 15.4 mag; Bloom et al. 2009; Racusin et al. 2009) and GRB 130427A (m_U ≈ 13.9 mag; Maselli et al. 2014) at comparable observer-frame times post-burst. But even applying only a correction for foreground (i.e., Milky Way) extinction of A_U ≈ 6.8 mag (Schlafly & Finkbeiner 2011), the UV–optical emission from GRB 221009A would likely have been the brightest observed to date had it occurred at high Galactic latitude (see also Kann et al. 2023).

If we adopt E(B − V)_host = 0.4 mag and an SMC-like extinction law (Section 3.2), we infer a U-band host extinction of 1.8 mag. Neglecting k-corrections (which are likely to be much smaller than the uncertainty in the extinction correction), we find an absolute magnitude of M_U ≈ − 29.3 mag (AB). Comparing to Figure 7 in Perley et al. (2014), we find that the UV–optical light curve of GRB 221009A occupies a similar location in luminosity space as the X-ray afterglow: toward the most luminous end, but entirely consistent with the existing distribution.

Although Swift was behind the Earth at T₀, we use the GBM light curve and spectrum (S. Lesage et al. 2023, in preparation) to estimate what GRB 221009A would have looked like to the BAT, as well as to determine the occurrence rate of such events. Adopting the time-averaged spectrum of the main pulse from T₀ + 218.5–277.9 s (10–1000 keV fluence of 8.293 × 10⁻² erg cm⁻²), we derive a BAT (15–150 keV band) fluence of 1.86 × 10⁻² erg cm⁻² s⁻¹ and an average flux of 3.13 × 10⁻⁴ erg cm⁻² s⁻¹. This fluence from the main pulse alone is 50 × higher than the largest fluence detected to date by BAT (GRB 130427A; Maselli et al. 2014).

At z = 0.151, the 15–150 keV (rest-frame) isotropic prompt energy release for the main pulse is E_γ,iso = 9 × 10⁵³ erg. We employ the 15–150 keV energy range to facilitate a direct comparison with the entire population of BAT GRBs with measured redshifts.

To estimate the relative occurrence rate of such an energetic event, we compare GRB 221009A with the intrinsic long GRB luminosity function derived in Lien et al. (2014). We consider only the main pulse, as this dominates the burst energetics, and it is difficult to ascertain when the prompt emission ends, given the afterglow was bright enough to trigger the BAT at T₀ + 3.4 ks.

We randomly generated 10⁴ GRBs using the intrinsic GRB rate and luminosity distribution in Lien et al. (2014), and each simulated burst was randomly assigned a pulse structure drawn from the real BAT GRB sample. We calculated the E_γ,iso of these simulated bursts, and found that only one burst had an energy release (slightly) higher than the main pulse of GRB 221009A. Therefore, we conclude only ∼1/10⁴ long GRBs are as energetic as GRB 221009A.

Using this value, we crudely estimate the rate at which such energetic GRBs are detectable by BAT:

$$\begin{eqnarray}&&{N}_{\mathrm{BAT}}\approx {R}_{\mathrm{GRB}}\ {f}_{\mathrm{Eiso}}\ {f}_{\mathrm{FOV}}\ {f}_{\mathrm{survey}}=0.06\ {\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 3 }$$

where R_GRB is the intrinsic all-sky long-GRB rate from Lien et al. (2014), f_Eiso = 1.0 × 10⁻⁴ is the fraction of GRBs in this intrinsic sample that have E_γ,iso larger than GRB 221009A, f_FOV = 1/6 is the fraction of sky covered by the BAT field of view, and f_survey = 0.8 is the fraction of time BAT is able to trigger on GRBs (neglecting spacecraft slews, SAA passage, etc.). In other words, we need to wait ∼1/0.06 = 17 yr for an event like the main pulse of GRB 221009A to occur in the BAT field of view.

We emphasize, however, that the above estimate simply indicates that BAT should have detected ∼1 GRB more energetic than GRB 221009A over its lifetime, independent of distance. But in addition to being highly energetic, GRB 221009A is also one of the most nearby GRBs detected by Swift. To estimate the intrinsic (volumetric) rate of comparable events, we utilize the BAT trigger simulator (Lien et al. 2014) to calculate the detectability of the prompt emission at different redshifts under different representative geometries. The full details of the simulation are provided in Appendix G.

Using this framework, we derive an upper limit on the local volumetric rate of GRB 221009A–like events of R_GRB,comov(z = 0) ≤ 6.1 × 10⁻⁴ Gpc⁻³ yr⁻¹. Integrating this flat comoving rate from z = 0, to z = 10, we obtain an upper limit on the all-sky rate of such events of ≤0.5 yr⁻¹. Comparing to the all-sky intrinsic long GRB rate of ∼4571 yr⁻¹ in Lien et al. (2014), the fraction of GRB 221009A–like events is roughly 0.5/4571 ≤ 1.0 × 10⁻⁴, which is similar to the relative rate derived above.

Even more remarkable, we can use these results to derive the rate of GRB 221009A–like GRBs within the volume out to z = 0.151 (1.1 Gpc⁻³). This implies we would need to wait over ≈10³ yr to detect another GRB 221009A–like event within this volume. The combination of the large energy release and the small distance make GRB 221009A truly a once in a lifetime phenomenon.

Interpreting GRB 221009A in the context of the standard afterglow synchrotron model (e.g., Sari et al. 1998) presents a number of challenges, in particular if the X-ray and optical emission is assumed to arise from a common origin (see Ghisellini et al. 2007; De Pasquale et al. 2009). To begin with, the broadband spectral fitting performed in Section 3.2 strongly favors the presence of a spectral break around the soft X-ray band, with the change in spectral slope being consistent with the cooling frequency ν_c . However, if we adopt relatively standard afterglow parameters (_B = 0.01, p = 2.5, η_γ = 0.15), then we find that a cooling frequency ≈5 keV (10¹⁸ Hz) at ≈0.1 day post-trigger implies a jet expanding into an extremely low density circumburst environment: n₀ ≈ 10⁻³ cm⁻³ for a constant-density (interstellar-medium-like) environment. Such low densities have been inferred previously for highly energetic events (e.g., Cenko et al. 2011), but this remains difficult to reconcile with the massive star progenitors of long GRBs.

Furthermore, standard afterglow closure relations (e.g., Racusin et al. 2009) are unable to reproduce the observed spectral and temporal power-law indices. For example, the broadband SED fits presented in Section 3.2 find an optical spectral index β_O ≈ 0.7. This is consistent with the temporal decay observed for a jet expansion into a constant-density medium (α₀ ≈ 1.1; p ≈ 2.4). However, this would predict a significantly shallower temporal decay in the X-rays (α_X ≈ 1.3) than what is observed. While expansion into a wind-like medium would predict a steeper X-ray temporal decay, the optical emission is inconsistent with this picture.

Given the exceptionally large E_γ,iso, GRB 221009A requires a narrow opening angle (and hence an early jet break time) in order to be compatible with the geometry-corrected energy release inferred from the pre-Swift sample (e.g., Frail et al. 2001). Assuming an efficiency of ≈15% in converting bulk kinetic energy into prompt γ-ray emission, a beaming-corrected energy release of ≤10⁵² erg requires an opening angle of θ_j ≲ 4°. For a top-hat jet expanding into a constant-density medium, this implies a jet break time of the following:

$$\begin{eqnarray}&&{t}_{j}\leqslant 3.5{\left(\displaystyle \frac{{n}_{0}}{0.1\,{\mathrm{cm}}^{-3}}\right)}^{-1/3}\,\mathrm{days}\end{eqnarray} \tag{ 4 }$$

where n₀ is the circumburst density in units of 0.1 cm⁻³ (Sari et al. 1999). Note that, due to the −1/3 exponent, this result is relatively insensitive to the exact value of the circumburst density.

With X-ray coverage extending out to T₀ + 70 days, the Swift+MAXI+NICER X-ray light curve is more than adequate to search for such a geometric signature. But out to late times, the temporal decay does not steepen beyond α_X = 1.7, significantly shallower than the α > 2 behavior anticipated for a simple on-axis top-hat jet.

The X-ray afterglow steepening observed at T₀ + 8 × 10⁴ s could conceivably be attributed to a jet break (see also D'Avanzo et al. 2022). In this case, the corresponding jet would indeed be extremely narrow, with an implied opening angle of ≈2°. However, as described above, the post-break decay index is significantly shallower than expected for an on-axis top-hat jet. Interestingly, similar behavior was observed in GRB 130427A and attributed to time-variable microphysical parameters (Maselli et al. 2014). One alternative explanation is a complex angular structure for the jet (i.e., compared to the simple top-hat model assumed above, with a constant energy as a function of angle). The presence of energetic material outside the jet "edges" would naturally account for the shallower decline initially (e.g., O'Connor et al. 2023). However, eventually when viewed sufficiently far off-axis, the decay should steepen (assuming the material is still relativistic), something not observed in the X-ray light curve out to T₀ + 70 days. A detailed exploration of the jet structure would require additional multiwavelength observations, in particular radio coverage, and is thus beyond the scope of this work.