Rocking the BOAT: the ups and downs of the long-term radio light curve for GRB 221009A

ABSTRACT

We present radio observations of the long-duration gamma-ray burst (GRB) 221009A that has become known to the community as the Brightest Of All Time or the BOAT. Our observations span the first 475 d post-burst and three orders of magnitude in observing frequency, from 0.15 to 230 GHz. By combining our new observations with those available in the literature, we have the most detailed radio data set in terms of cadence and spectral coverage of any GRB to date, which we use to explore the spectral and temporal evolution of the afterglow. By testing a series of phenomenological models, we find that three separate synchrotron components best explain the afterglow. The high temporal and spectral resolution allows us to conclude that standard analytical afterglow models are unable to explain the observed evolution of GRB 221009A. We explore where the discrepancies between the observations and the models are most significant and place our findings in the context of the most well-studied GRB radio afterglows to date. Our observations are best explained by three synchrotron-emitting regions that we interpret as a forward shock, a reverse shock, and an additional shock potentially from a cocoon or wider outflow. Finally, we find that our observations do not show any evidence of any late-time spectral or temporal changes that could result from a jet break but note that any lateral structure could significantly affect a jet break signature.

gamma-ray burst: individual: GRB 221009A, ISM: jets and outflows, radio continuum: transients

1 INTRODUCTION

Long-duration gamma-ray bursts (GRBs) are produced in highly relativistic jets, launched during the collapse of massive stars, and they are the most powerful explosions in the Universe. GRB 221009A has been dubbed the Brightest Of All Time or the BOAT (Burns et al. 2023). Lasting about 600 s, the variable, high-energy, prompt emission was detected by the Neil Gehrels Swift Observatory – Burst Alert Telescope and X-Ray Telescope (BAT and XRT, respectively; Williams et al. 2023), Insight-Hard X-ray Modulation Telescope (HXMT) and Gravitational wave high-energy Electromagnetic Counterpart All-sky Monitor C (GECAM-C; An et al. 2023), Konus-Wind and Spectrum-Roentgen-Gamma (SRG)/Astronomical Roentgen Telescope – X-ray Concentrator (ART-X; Frederiks et al. 2023), and the Fermi Gamma-ray Space Telescope (Lesage et al. 2023). Placed at a redshift of 0.151 (de Ugarte Postigo et al. 2022; Malesani et al. 2023), the isotropic gamma-ray energy output has been measured as |$1\times 10^{55}$|⁠, |$1.5\times 10^{55}$|⁠, and |$1.2\times 10^{55}$| erg by An et al. (2023), Frederiks et al. (2023), and Lesage et al. (2023), between 1 and 10 000 keV, 10 keV and 6 MeV, and 20 keV and 10 MeV, respectively, nearly twice the value of the next most energetic, GRB 080916C (Greiner et al. 2009). Given its prompt emission properties, it has been established as a once in 10 000 yr event (Burns et al. 2023). In fact, GRB 221009A was so bright that the prompt emission caused disturbances in the ionosphere (Hayes & Gallagher 2022).

The afterglow to GRB 221009A has been detected consistently between 0.4 GHz and 20 TeV (Laskar et al. 2023; LHAASO Collaboration 2023). In terms of spectral coverage, it exceeds all other TeV afterglows with radio detections like GRB 190114C or GRB 190829A (MAGIC Collaboration 2019a; H. E. S. S. Collaboration 2021). In terms of data quantity and quality, it exceeds the GHz-to-GeV afterglow of GRB 130427A (e.g. Ackermann et al. 2014; Levan et al. 2014; van der Horst et al. 2014; de Pasquale et al. 2016), although the latter had a much better sampling of optical light curves since it did not suffer from extinction in the way that GRB 221009A did (Fulton et al. 2023; Levan et al. 2023). Similar to GRB 130427A (Anderson et al. 2014), Bright et al. (2023) showed that GRB 221009A had a bright light-curve peak at 15 GHz within the first day, followed by an overall decline at radio frequencies. This behaviour is quite different from the ‘classical’ well-sampled radio afterglows of, for instance, GRB 970508 and GRB 030329, which have peaks at time-scales of weeks to months (Frail, Waxman & Kulkarni 2000; Resmi et al. 2005; van der Horst et al. 2005). The origin of the early time radio peaks is thought to be from reverse shock, produced from a shock front propagating back through the jet. Details of the light-curve behaviour, in particular over a wide frequency range, give important insights into the underlying physics at various scales, from the jetted explosion outflow to the accelerated particles generating the observed emission (e.g. Sari, Piran & Narayan 1998; Wijers & Galama 1999).

The focus of this paper is the radio emission from GRB 221009A, covering three orders of magnitude in both observing frequency and days post-burst. While the TeV emission leads to various questions regarding possible emission processes at these high energies and the potential for detecting GRBs at TeV energies more frequently (MAGIC Collaboration 2019a; H. E. S. S. Collaboration 2021; Abe et al. 2024), the radio observations provide the necessary context for understanding the physics of the jetted GRB outflow, together with multiwavelengths observations in the optical and X-rays (e.g. Gill & Granot 2023; O’Connor et al. 2023). The light-curve behaviour of GRB 221009A in the first days to weeks does not seem to follow expectations of the standard model that is typically used to describe radio afterglows (e.g. Wijers & Galama 1999; Granot & Sari 2002). The extremely dense sampling of the light curves at various radio frequencies as presented in this paper is unprecedented and allows for detailed modelling that will lead to better descriptions of GRB jets and the relevant emission processes.

The dominant emission mechanism in GRB afterglows at radio frequencies is synchrotron radiation from extremely relativistic electrons accelerated by shocks at the front of a relativistic collimated outflow (Meszaros & Rees 1993; Sari et al. 1998). This is also the emission mechanism assumed to be at play in the GRB 221009A afterglow. While we are only considering one emission mechanism, i.e. synchrotron, there can be multiple emission sites. For instance, the jet sweeping up particles in the ambient medium leads to a forward shock, but will also lead to the formation of the aforementioned reverse shock that can be dominant at early times given the right conditions (Kobayashi & Sari 2000). Besides this shock structure in the radial direction, there can also be structure in the lateral direction. This structure could be smooth, for instance, a structured energy profile as a function of distance to the jet axis instead of a homogeneous energy profile (Rossi, Lazzati & Rees 2002; Lamb et al. 2021; Salafia & Ghirlanda 2022); but there could also be multiple jet components (Starling et al. 2005), and potentially a cocoon around the jet (Ramirez-Ruiz, Celotti & Rees 2002; Nakar & Piran 2017; Izzo et al. 2019). This could lead to multiple synchrotron emission components or emission components that evolve differently from the canonical top-hat behaviour (van der Horst et al. 2014; Bright et al. 2019; Rhodes et al. 2022).

Besides these macrophysical considerations, high-quality multiwavelength data as presented here reveal nuances in the microphysics of GRB afterglows. Afterglow modelling can lead to insights into the magnetic field strength and energetics, but also the total energetics, acceleration efficiency, and energy distribution of the accelerated electrons (Granot & Sari 2002; Eichler & Waxman 2005). To complicate this further, detailed simulations of particle acceleration and magnetic field amplification by relativistic shocks indicate that there is potentially a time dependence of the energies in magnetic fields and electrons (Rossi & Rees 2003), and this has also been adopted in multiwavelength modelling of some GRB afterglows with peculiar behaviour (van der Horst et al. 2014; Bright et al. 2019; Misra et al. 2021; Salafia et al. 2022).

Given the extremely high quality of the radio data presented in this paper and the dynamics of the synchrotron spectrum that is likely quite different from the standard behaviour, we take a fairly cautious approach in the modelling presented here. While a standard GRB synchrotron spectrum is still assumed, the temporal evolution of the spectrum is kept free of constraints where possible, to provide input on detailed modelling and theoretical efforts, and get a better handle on the interpretation of the wealth of these data from this unique source. We highlight here the use of the convention |$F_{\nu } \propto t^{\alpha }\nu ^{\beta }$| throughout this work to describe the temporal and spectral evolution. This paper is laid out in the following manner: in Section 2, we present the new radio observations and the data reduction methods used; in Section 3, we lay out the results of our observing campaigns and describe the model used to explain the data; in Section 4, we put our results in a broader context and interpret the data using various models; and we conclude in Section 5.

2 OBSERVATIONS

Here, we present the data reduction processes for the observations used in this work. The flux density measurements and upper limits for our new observations are summarized in Table 1. In addition to the data sets we present here, our work also incorporates the previously published radio data from Giarratana et al. (2023), Laskar et al. (2023), and Bright et al. (2023), and the X-ray data from Williams et al. (2023).

Table 1.

A table of the new radio observations presented in this work. All non-detections are indicated by a ‘–’ in the flux density column followed by the |$3\sigma$| upper limit in the uncertainty column. The full list of radio observations is presented in supplementary material online.

Obs date (MJD)	Observing frequency	Flux density (mJy)	Uncertainty (mJy)	Telescope	T − T0 (d)
59866.65	15.50	7.18	0.36	AMI-LA	5.097
59866.84	15.50	7.05	0.36	AMI-LA	5.284
59867.66	15.50	6.68	0.34	AMI-LA	6.106
59867.84	15.50	6.84	0.35	AMI-LA	6.283
59868.62	15.50	5.99	0.31	AMI-LA	7.065
59869.82	15.50	5.59	0.28	AMI-LA	8.266
...	...	...	...	...	...

2.1 AMI-LA

The Arcminute Microkelvin Imager – Large Array (AMI-LA) is an eight-dish interferometer based in Cambridge, UK (Zwart et al. 2008). It observes at a central frequency of 15.5 GHz with a bandwidth of 5 GHz, achieving an angular resolution of about 30 arcsec (Hickish et al. 2018). Bright et al. (2023) presented the first five days of observations from AMI-LA, and here we present the rest of the observing campaign. We continued to observe the position of GRB 221009A almost daily until 210 d post-burst when the first non-detection occurred. Between 210 and 320 d post-burst, we concatenated separate non-detections to obtain deeper limits.

AMI-LA data are reduced using a custom software package: reduce_dc (Perrott et al. 2013). The software performs bandpass and flux scaling using 3C 286 and complex gain calibration using J1925+2106. Flagging and imaging are done in casa using the tasks rflag, tfcrop, and clean (McMullin et al. 2007). The details of observing times and measured flux densities are provided in Table 1. We note that unlike in Bright et al. (2023), we do not split each observation up, because the duration of a given epoch is a negligible fraction of the total time since the burst was first detected, so no significant evolution is expected within an observation.

2.2 ASKAP

We obtained target-of-opportunity observations of the GRB 221009A field with the Australian Square Kilometre Array Pathfinder (ASKAP; Johnston et al. 2007). Our observations were centred on 888 MHz, with a bandwidth of 288 MHz, taken using the square_6x6 beam footprint (see fig. 20 of Hotan et al. 2021). The data products for these observations can be found under the project code AS113 with SBIDs: 44780, 44857, 44918, 45060, 45086, 45416, 46350, 46419, 46492, 46554, and 48611 in the CSIRO ASKAP Science Data Archive (CASDA).¹

Observations of PKS B1934−638 were used to calibrate the antenna gains, bandpass, and the absolute flux density scale. Flagging of radio frequency interference, calibration of raw visibilities, full-polarization imaging, and source finding on total intensity images were all performed through the standard askapsoft pipeline (Guzman et al. 2019). The resulting image reached a typical rms of |${\sim}50\, \mu \mathrm{ Jy}$| beam⁻¹. We evaluated and corrected for the systematic flux scale offset by comparing the flux density of field sources in each observation against those in the Rapid ASKAP Continuum Survey (RACS) catalogue (Hale et al. 2021).

2.3 ATA

Located |${\sim}200$| miles north of San Francisco, the Allen Telescope Array (ATA) is a 42-element radio interferometer hosted at the Hat Creek Radio Observatory. Mounted on the focus of each element is a dual-polarization, log-periodic feed that is cryogenically cooled and sensitive to radiation in the range of 1–12 GHz. Analogue signals from the array are transmitted through a fibre to a centralized signal processing room and are split into four independent chains that get multiplexed by four tunable local oscillators in a superheterodyne system. The current correlator backend supports the digitization of two out of the four available tunings for 20 of the 42 antennas, where each tuning can be placed anywhere in the available radio frequency range of the log-periodic feed, with |${\sim}700$| MHz of usable bandwidth for each.

The radio counterpart of GRB 221009A was observed extensively with the ATA beginning just a few hours after the burst as reported in Bright et al. (2023). Here, we build on that work and utilize the flexible frequency tunability of the ATA to monitor the |$1\!-\!10\, \rm {GHz}$| spectral evolution over its entire outburst. Either 3C 147, 3C 48, or 3C 286 was observed as flux calibrator at the beginning of each observing block, and a 10 min observation of the phase calibrator J1925+2106 was interleaved for every 30 min of science target recording (regardless of observing frequency). We evolved our total integration time on source over the course of the follow-up campaign to account for the fading of GRB 221009A. Visibilities from each observation block were reduced using a custom pipeline using aoflagger (Offringa 2010) and casa (McMullin et al. 2007). Images for the flux, phase, and science targets were formed using standard casa tasks and by deconvolving with the clean algorithm (Högbom 1974; Clark 1980; Sault & Wieringa 1994). We used two Taylor terms to account for the high fractional bandwidth (especially at low frequencies) and a Briggs robust value of 0.5 when imaging. Finally, flux densities for GRB 221009A were derived by fitting a point source (i.e. with a source size fixed to the dimensions of the main lobe of the dirty beam) to the science target.

2.4 ATCA

We carried out multiple observations of the radio counterpart to GRB 221009A using the Australia Telescope Compact Array (ATCA) under the project codes: CX515 (Director’s Discretionary Time), C3374 (PI: G. E. Anderson), and C3542 (PI: G. E. Anderson). These observations were carried out using the 5.5/9, 16.7/21.2, 33/35, and 43/45 GHz receiver configurations, with a bandwidth of 2048 MHz for each intermediate frequency.

For each observation, we reduced the visibility data using standard procedures in miriad (Sault, Teuben & Wright 1995). We used a combination of manual and automatic radio frequency interference flagging before calibration. For bandpass calibration, we used PKS B1934−638 at 5.5/9 GHz, while at higher frequencies (16.7/21.2, 33/35, and 43/45 GHz) we used either B1921−293 or B1253−055; the spectral shape of B1921−293 and B1253−055 was accounted for by fitting to first order the measured flux densities of these calibrators at each intermediate frequency for each of the higher frequency observing bands. The flux density scale was set using B1934−638 for all observing frequency bands. For all observations, we used B1923+210 to calibrate for the time-variable complex gains. After calibration, where there was sufficient signal-to-noise ratio, we split the 2048 MHz bandwidth into further subbands to obtain higher spectral resolution. We then inverted the visibilities and applied the multifrequency synthesis clean algorithm (Högbom 1974; Clark 1980; Sault & Wieringa 1994) to the target source field using standard tasks in miriad to obtain our final images. The flux densities of the radio afterglow candidate were extracted by fitting a point source to the radio source, in the case of a detection, while, in the case of a non-detection, the limits were obtained using the rms sensitivity in the residual image.

2.5 e-MERLIN

The enhanced Multi-Element Remotely Linked Interferometer Network (e-MERLIN) is a radio interferometer made up of seven dishes spread across the UK. With a maximum baseline of 217 km, whilst observing at 5 GHz, it can resolve angular scales of 0.05 arcmin. We observed the position of GRB 221009A with e-MERLIN through a combination of rapid response time requests (PI: L. Rhodes; RR14001) and open time proposals (PI: L. Rhodes; CY13003, CY14001, and CY15206) at both L and C bands. Our L- and C-band observations were centred at 1.51 and 5.08 GHz, respectively, both with a bandwidth of 512 MHz. We note that the first two epochs obtained at L-band have previously been published in Bright et al. (2023).

All observations were reduced using the e-MERLIN pipeline within casa (McMullin et al. 2007; Moldon 2021). The pipeline performs preliminary flagging for radio frequency interference and known observatory issues. It then performs two rounds of bandpass calibration and complex gain calibration, using OQ 208 and J1905+1943, respectively, along with flux scaling using 3C 286. Further flagging of the target field is conducted. We performed interactive cleaning and deconvolution using the casa task tclean.

2.6 LOFAR

8 h of Director’s Discretionary Time with the LOw Frequency ARray (LOFAR; DDT|$20\_003$|⁠) were awarded to observe GRB 221009A. The allocated time was split into two observing runs of 4 h, which took place on 2023 July 18 and 20 at matching local sidereal times. Each observing run was preceded by a 10 min calibrator scan of 3C 295. All observations were conducted in the HBA_dual_inner configuration where, in addition to the 22 core stations available, the inner tiles of 14 remote stations were also used. The single-beam observations were centred at 152.05 MHz with 380 subbands, and data were recorded with an integration time of 1 s. Each subband consisted of 64 frequency channels of width 3.051 kHz. The data were subsequently averaged to 16 channels of 12.21 kHz per subband by the observatory during data pre-processing. Both target observations were calibrated for direction-independent effects using linc² with default settings, a pipeline developed by the observatory to correct for various instrumental and ionospheric effects present in interferometric LOFAR data (van Weeren et al. 2016; Williams et al. 2016; de Gasperin et al. 2019). Due to its relative proximity, Cygnus A was subtracted from the visibilities using the ‘demixing’ step in linc. The data were further averaged to four channels of 48.82 kHz per subband and 4 s during calibration. The resulting calibrated data were concatenated into groups of 20 subbands and averaged in time to 8 s. These data products from both observations were subsequently jointly put through ddf-pipeline³ for direction-dependent calibration and imaging (Shimwell et al. 2019; Tasse et al. 2021). This resulted in a final image generated using a circular restoring beam of radius 3 and 1.5 arcsec pixel resolution.

2.7 NOEMA

The NOrthern Extended Millimetre Array (NOEMA, situated in the southern French Alps) monitored GRB 221009A between 2022 October 10 and 2023 April 25 in the 3, 2, and 1 mm bands. Interferometer configurations were medium-extended C and extended A configurations with up to 12 antennas, primary flux calibrators were MWC 349 and LKHA 101. The data were reduced with the clic and mapping software packages that are part of the gildas⁴ package. Fluxes and their errors were derived from point-source ultraviolet (UV) plane fits to the calibrated interferometric visibilities.

2.8 uGMRT

We observed GRB 221009A with the upgraded Giant Metrewave Radio Telescope (uGMRT) in bands 5 (1000–1450 MHz) and 4 (550–900 MHz) under a DDT proposal (ddtC251, PI: P. Chandra). The observations were made at two epochs in both bands, once in 2023 January and then in 2023 March. We recorded the data in 2048 frequency channels covering a bandwidth of 400 MHz with an integration time of |${\sim}10 \, \mathrm{ s}$|⁠. We used 3C 286 and 3C 48 as flux density and bandpass calibrators. J1924+3329 was used as a phase calibrator.

The data were analysed using the casa package (McMullin et al. 2007) following the procedure in Nayana et al. (2022). We also performed a few rounds of phase only and one round of amplitude and phase self-calibration to improve the image quality. The final flux densities were obtained by fitting a Gaussian at the GRB position.

3 RESULTS AND MODEL

There have been several GRB 221009A afterglow modelling efforts that have used a subset of the radio data published to date (including but not limited to Gill & Granot 2023; Laskar et al. 2023; Levan et al. 2023; O’Connor et al. 2023). Here, we present the results of our observing campaigns and describe out modelling of the radio and X-ray afterglow.

3.1 Light curves and SEDs

The radio data presented in this paper span three orders of magnitude in frequency space, from 0.15 MHz to 230 GHz, and lasts out to 475 d post-burst. Fig. 1 shows the radio afterglow light curves split by observing frequency. Symbols with lower opacity denote all previously published data, whereas the solid symbols mark data presented in this paper. We include all previous and newly published data to extract the clearest scenario of the afterglow.

Figure 1.

Radio afterglow light curves of GRB 221009A split by observing frequency (or frequency range). Any low-opacity data points are from previously published observations. All observations presented in this paper are shown with solid circles for detections and downward-facing triangles for |$3\sigma$| upper limits.

Above 19 GHz, the afterglow is decaying at all times, with observations obtained between 1 and 200 d post-burst (the top two rows of Fig. 1). The light curves between 90 and 105 GHz in Fig. 1 show that the decay rate slowly steepens with time like a very smooth broken power law. Below 16 GHz, we observe the light-curve peak in almost each observing band, except at 9–10 and 0.4 GHz since we were not observing early enough at those frequencies. Bright et al. (2023) interpreted this peak as emanating from the reverse shock, which we are tracking from 17.7 to below 1 GHz. The data between 1.3 and 3 GHz also show a second, distinct bump at around 50 d. In addition to the early peaks caught at 5 and 15.5 GHz, we also see evidence of further bumps during the decay phase. It is possible that the additional bumps originate from different spectral components.

Fig. 2 shows the broad-band radio spectral energy distributions (SEDs) throughout our campaign. For the first 30 d, a low-frequency turnover is visible and the below-turnover spectral index is consistent with |$\beta \sim 5/2$| below the turnover. Above the turnover, we find a flat spectrum extending to the highest frequencies (⁠|${\sim}200$| GHz). A flat spectrum is inconsistent with optically thin synchrotron emission from a single component and so provides further evidence of multiple spectral components, similar to GRB 130427A (Perley et al. 2014). Only after 150 d post-burst does the spectrum steepen with typical optically thin spectral indices (⁠|$\beta \sim -0.5$| to |$-1$|⁠), more consistent with that from the late-time X-ray data (Williams et al. 2023). Williams et al. (2023) performed a joint fit to the UV, X-ray, and gamma-ray data, which shows that the high-energy spectra can be described by either a single power law or a broken power law where the break, interpreted as the synchrotron cooling break |$\nu _{\textrm {c}}$|⁠, sits in the XRT band. The broken power law is favoured but the fits are only performed on data up to 1 d post-burst, whereas the X-ray light curve itself extends out to 200 d post-burst.

Figure 2.

Broad-band radio SEDs for GRB 221009A as a function of time. As in Fig. 1, low-opacity data points denote previously published data, while solid points are observations presented in this paper. Because epochs have been chosen to demonstrate the spectral evolution, we note that not all data presented in Fig. 1 are also shown here.

3.2 Modelling

Here, we build on previous modelling efforts by combining our new observations from AMI-LA, ATA, ATCA, ASKAP, e-MERLIN, LOFAR, NOEMA, and uGMRT with radio data available in the literature. We also include the full Swift-XRT light curve (in flux densities at 10 keV; Williams et al. 2023). We do not include any optical or other high-energy data in our modelling work as there are too many contaminating components in these bands. At optical frequencies, there is significant extinction (Tiengo et al. 2023; Vasilopoulos et al. 2023) both from the Milky Way and the host galaxy, plus a contribution from the associated supernova. Above keV energies, there is an increasing contribution from the additional very high energy component whose origin and emission mechanism is still debated (Aharonian et al. 2023; LHAASO Collaboration 2023; Savchenko et al. 2024).

We consider models that use either two or three synchrotron spectral components that can evolve independently in time to explain the behaviour shown in the light curves (Fig. 1) and SEDs (Fig. 2). Each synchrotron spectrum is constructed of four power-law slopes divided by three frequency breaks: the synchrotron self-absorption break (⁠|$\nu _{\textrm {sa}}$|⁠), the characteristic or minimum electron energy break (⁠|$\nu _{\textrm {m}}$|⁠), and the cooling break (⁠|$\nu _{\textrm {c}}$|⁠, above which radiative cooling is important). The peak of the spectrum, |$F_{\nu , \textrm {max}}$|⁠, is at whichever frequency break of |$\nu _{\textrm {sa}}$| or |$\nu _{\textrm {m}}$| is higher. The spectral index of each branch depends on the order of the frequency breaks. In the regime where |$\nu _{\textrm {sa}} \lt \nu _{\textrm {m}} \lt \nu _{\textrm {c}}$|⁠, the spectral indices are |$F_{\nu \lt \nu _{\textrm {sa}}} \propto \nu ^{2}$|⁠, |$F_{\nu _{\textrm {sa}} \lt \nu \lt \nu _{\textrm {m}}} \propto \nu ^{1/3}$|⁠, |$F_{\nu _{\textrm {m}} \lt \nu \lt \nu _{\textrm {c}}} \propto \nu ^{(1-p)/2}$|⁠, and |$F_{\nu _{\textrm {c}} \lt \nu } \propto \nu ^{-p/2}$|⁠, where p is the electron energy distribution index and is typically expected to be between 2 and 3 (although values slightly below 2 and above 3 have been reported; Kirk et al. 2000; Achterberg et al. 2001; Sironi, Spitkovsky & Arons 2013). In the regime where |$\nu _{\textrm {m}} \lt \nu _{\textrm {sa}} \lt \nu _{\textrm {c}}$|⁠, the spectral indices are |$F_{\nu \lt \nu _{\textrm {m}}} \propto \nu ^{2}$|⁠, |$F_{\nu _{\textrm {m}} \lt \nu \lt \nu _{\textrm {sa}}} \propto \nu ^{5/2}$|⁠, |$F_{\nu _{\textrm {sa}} \lt \nu \lt \nu _{\textrm {c}}} \propto \nu ^{(1-p)/2}$|⁠, and |$F_{\nu _{\textrm {c}} \lt \nu } \propto \nu ^{-p/2}$|⁠. As the jet expands and evolves, the spectral breaks are expected to change as a power-law function of time, which depends on the jet dynamics and the density profile through which the jet is propagating, |$\rho \propto r^{-k}$|⁠, where |$k = 0$| for a homogeneous medium and |$k = 2$| represents a stellar wind (Granot & Sari 2002; Granot & van der Horst 2014).

We use emcee to fit our respective models to the data (Foreman-Mackey et al. 2013). Each model uses 40 walkers and runs for at least 70 000 steps or until convergence. All priors are uniform, and the only priors with fixed bounds were |$p \in [1.5,3.5]$| to help rule out unphysical solutions. The best-fitting value for each parameter is the 50th percentile post-burn-in of the posterior distribution, and the 84th and 16th percentiles are quoted as the upper and lower uncertainties, respectively.

3.2.1 Two-component model

First, we fit the data with two separate synchrotron spectra. The first is the reverse shock identified in Bright et al. (2023), we find that the peak of the synchrotron spectrum is produced by |$\nu _{\textrm {sa}}$| and fit for the normalization and evolution of the spectrum as well as p. The second component is a forward shock that appears to dominate the optical and X-rays (e.g. Fulton et al. 2023; Shrestha et al. 2023; Williams et al. 2023), and also the late-time radio emission. Here, we allow both |$\nu _{\textrm {sa}}$| and |$\nu _{\textrm {m}}$| to vary freely. We fit for the normalization and evolution of |$F_{\nu , \textrm {max}}$|⁠, |$\nu _{\textrm {sa}}$|⁠, and |$\nu _{\textrm {m}}$| as well as p. The resulting model parameters are provided in Table 2.

Table 2.

The parameter values (50th percentile) and their associated uncertainties (18th and 64th percentiles) derived for our best-fitting two-component model. Any |$\alpha$| parameter refers to the temporal power-law index of the parameter written in the subscript, as described in Section 3. For the reverse, forward, and extra shock component, |$F_{\nu , \textrm {max}}$| and |$\nu _{\textrm {sa}}$| are normalized to 1 and 6.5 d, respectively. For each shock, p is the value of the electron energy spectral index.

Parameter	Value
Reverse shock
|$F_{\nu , \textrm {max}}$| (mJy) [1 d]	24.0 |$\pm$| 0.8
|$\nu _{\textrm {sa}}$| (GHz) [1 d]	6.3 |$\pm$| 0.1
|$\alpha _{{F}_{\nu , \textrm {max}}}$|	−0.84 |$\pm$| 0.02
|$\alpha _{\textrm {sa}}$|	−0.957 |$\pm$| 0.008
p	|${\lt} 1.5$|
Forward shock
|$F_{\nu , \textrm {max}}$| (mJy) [6.5 d]	3.10 |$\pm$| 0.06
|$\log (\nu _{\textrm {sa}})$| (GHz) [6.5 d]	−0.5 |$\pm$| 0.1
|$\log (\nu _{\textrm {m}})$| (GHz) [6.5 d]	2.20 |$\pm$| 0.04
|$\alpha _{{F}_{\nu , \textrm {max}}}$|	−0.63 |$\pm$| 0.02
|$\alpha _{\textrm {m}}$|	−1.67 |$\pm$| 0.03
|$\alpha _{\textrm {sa}}$|	−0.11 |$\pm$| 0.07
p	2.32 |$\pm$| 0.03

We find that the two-component model cannot reproduce the flat spectrum observed shown in Fig. 2, the posterior distribution of p for the reverse shock always ends up at the lower bound of the prior with values for p below 1.5 or even below 1, and such a low value is unphysical and so we no longer consider this scenario.

3.2.2 Three-component model

Given the issues with a two-component model, we include a third component to alleviate the shallow value of p that was needed in the two-component model to explain the flat spectrum that is present during the first |${\sim}150$| d (Fig. 2) and the additional bumps in the 5 and 15.5 GHz light curves around 5–10 d post-burst (see Fig. 1).

To best explore the parameter space of the third component, first, we test both |$\nu _{\textrm {sa}}$| and |$\nu _{\textrm {m}}$| as the peak frequency of the third component and find that |$\nu _{\textrm {sa}}$| provides a better fit. Then, we consider two different iterations of this extra shock with differing degrees of freedom, which are summarized in Table 3, in addition to the two shock components described in the previous section. In both iterations of our three-component model, the peak flux density of each of the three components follows a smoothly broken power law (Rhodes et al. 2020):

$$\begin{eqnarray} F_{\nu } = F_{\nu , \textrm {max}} \left(0.5\left(\frac{t}{t_{\mathrm{ b}}} \right)^{-\alpha _{1}s} +0.5\left(\frac{t}{t_{\mathrm{ b}}} \right)^{-\alpha _{2}s} \right)^{-\frac{1}{s}}, \end{eqnarray}$$

(1)

Table 3.

Summary of the different iterations of the three-component model that explore the possible evolution of the third shock component. Each |$\alpha$| corresponds to a temporal index of the subscripted value, e.g. |$\alpha _{F_{\nu , \textrm {max}, 1}}$| corresponds to the first temporal index used to describe the behaviour of |$F_{\nu , \textrm {max}}$|⁠. We find that model 1, combined with a forward and reverse shock, describes the data best. The model is shown compared to the data in Figs 4 and 5. The best-fitting parameters are shown in Table 4.

Model |$\#$|	|$\alpha _{F_{\nu , \textrm {max}, 1}}$|	|$\alpha _{F_{\nu , \textrm {max}, 2}}$|	|$\alpha _{\nu _{\textrm {sa}}, 1}$|	|$\alpha _{\nu _{\textrm {sa}}, 2}$|
1	3	|$\alpha$|	|$\alpha _{\textrm {sa}}$|	–
2	|$\alpha _{1}$|	|$\alpha _{2}$|	|$\alpha _{\textrm {sa}, 1}$|	|$\alpha _{\textrm {sa}, 2}$|

where |$F_{\nu , \textrm {max}}$| is the flux density at the break time |$t_{\mathrm{ b}}$|⁠, |$\alpha _{1}$| and |$\alpha _{2}$| are the power-law indices, and s is the smoothing parameter that we set to be 0.5. In model 1, the synchrotron self-absorption break follows a single power law: |$\nu _{\textrm {sa}} = \nu _{\textrm {sa}, 0} t^{\alpha _{\textrm {sa}}}$|⁠, where |$\nu _{\textrm {sa}, 0}$| is the location of the self-absorption break at 1 d post-burst. We set |$\alpha _{F_{\nu , \textrm {max}, 1}} = 3$| and both |$\alpha _{F_{\nu , \textrm {max}, 2}}$| (defined in Table 3 as |$\alpha$|⁠) and |$\alpha _{\textrm {sa}}$| can vary freely. We invoke a |$\alpha _{F_{\nu , \textrm {max}, 1}} = 3$| as done in Peng, Königl & Granot (2005), which is used in the regime where a blast wave that is initially off-axis has undergone significant deceleration and so the radiation begins to enter the observers’ line of sight. In their paper, they do not consider the self-absorption break, but we find it fits well within the constraints of our work. Ryan et al. (2020) also consider off-axis afterglows from a numerical perspective and find steeper rise rates for ‘far off-axis events’. We choose to be more conservative and use Peng et al. (2005) value.

In model 2, both the peak flux density and |$\nu _{\textrm {sa}}$| are both described with broken power laws where all the indices are fit for but the break time is the same. A full summary of our models to explain the extra forward shock is shown in Table 3.

Fig. 3 shows the results of the different iterations of our models. Unfortunately, not one of our models provides a perfect fit to the data, this may be due to a combination of unknown systematic uncertainties, and perhaps more importantly, this exquisite data set is showing evidence of more complicated physics and emission mechanisms that cannot be accounted for by the basic synchrotron models. As a result, we find quoting Bayesian evidence values inappropriate. However, we do find that model 1 provides the best fit. This is because our posterior distributions for all values of p sit between 2 and 3 and do not require such uncomfortably large temporal index values. We present the parameters of this fit in Table 4. Figs 4 and 5 show our best-fitting model overlaid on the light curves and SEDs. Fig. 4 shows that our model describes the long-term evolution at all frequencies well. However, it cannot replicate the bumps and wiggles observed at 15.5, 5, and 0.4 GHz, despite that being one of the motivations for the three-component model. Furthermore, it marginally overpredicts the late time 0.8 GHz flux. Fig. 5 demonstrates that the superposition of multiple components recreates the high-frequency emission accurately and describes well the flat spectrum and broad turnover at earlier times post-burst. On the other hand, we find that it tends to place the |$\nu _{\textrm {sa}}$| much lower than the observed position.

Figure 3.

Evolution of the break frequencies and peak flux for the three-component model. Each panel corresponds to a different iteration of our model as described in Section 3 and Table 3. For each iteration, we show only the average value (50th percentile value) of the posterior distribution for clarity. The left-hand vertical axis of each plot corresponds to the evolution of the frequency breaks (dotted and dashed lines for |$\nu _{\textrm {sa}}$| and |$\nu _{\textrm {m}}$|⁠, respectively). The right-hand vertical axis shows the evolution of the peak flux (solid lines) of each shock component.

Figure 4.

The multifrequency radio light curves for GRB 221009A overlaid with our best-fitting three-component model (model 1).

Figure 5.

Broad-band radio SEDs for GRB 221009A as a function of time with our best-fitting three-component model (model 1) overlaid.

Table 4.

The parameter values (50th percentile) and their associated uncertainties (18th and 64th percentiles) derived for our best-fitting three-component model (model 1). Any |$\alpha$| parameter refers to the temporal power-law index of the parameter written in the subscript, as described in Section 3 and Table 3. For the reverse, forward, and extra shock component, |$F_{\nu , \textrm {max}}$| and |$\nu _{\textrm {sa}}$| are normalized to 1 d, 6.5 d, and |$t_{\textrm {dec}}$|⁠, respectively, where |$t_{\textrm {dec}}$| is a parameter we fitted for. For each shock, p is the value of the electron energy spectral index.

Parameter	Value
Reverse shock
|$F_{\nu , \textrm {max}}$| (mJy) [1 d]	|$9.6^{+1.6}_{-1.5}$|
|$\nu _{\textrm {sa}}$| (GHz) [1 d]	|$4.4^{+0.2}_{-0.3}$|
|$\alpha _{{F}_{\nu , \textrm {max}}}$|	−0.59 |$\pm$| 0.05
|$\alpha _{\textrm {sa}}$|	−0.86 |$\pm$| 0.03
p	|$2.2^{+0.4}_{-0.3}$|
Forward shock
|$F_{\nu , \textrm {max}}$| (mJy) [6.5 d]	4.2 |$\pm$| 0.2
|$\log (\nu _{\textrm {sa}})$| (GHz) [6.5 d]	0.3 |$\pm$| 0.2
|$\log (\nu _{\textrm {m}})$| (GHz) [6.5 d]	2.71 |$\pm$| 0.08
|$\alpha _{{F}_{\nu , \textrm {max}}}$|	−0.97 |$\pm$| 0.03
|$\alpha _{\textrm {m}}$|	−1.06 |$\pm$| 0.06
|$\alpha _{\textrm {sa}}$|	|$-1.4^{+0.2}_{-0.1}$|
p	2.32 |$\pm$| 0.03
Extra shock
|$F_{\nu , \textrm {max}}$| (mJy) [|$t_{\textrm {dec}}$| d]	17 |$\pm$| 2
|$\nu _{\textrm {sa}}$| (GHz) [|$t_{\textrm {dec}}$| d]	|$1.03^{+0.05}_{-0.04}$|
|$\alpha$|	−0.71 |$\pm$| 0.02
|$\alpha _{\textrm {sa}}$|	|$-0.46^{+0.03}_{-0.04}$|
|$t_{\textrm {dec}}$| (d)	0.27 |$\pm$| 0.02
p	3.1 |$\pm$| 0.3

4 DISCUSSION

Here, we discuss the implications of our best-fitting three-component model and place them in the context of other detailed radio studies of GRBs.

4.1 Reverse shock

The dashed lines in Figs 4 and 5 show the contribution of the reverse shock from our model. Bright et al. (2023) used radio observations in the first five days post-burst to measure the evolution of F_max and |$\nu _{\textrm {sa}}$| with time. They found that |$F_{\nu , \textrm {max}} \propto t^{-0.70\pm 0.02}$| and |$\nu _{\textrm {sa}} \propto t^{-1.08\pm 0.04}$|⁠, and concluded that the evolution of the spectral peak was too slow to match theoretical predictions and most likely a superposition of multiple emitting regions. When considering the full radio data set, we find a different, even slower reverse shock evolution: |$F_{\rm max} \propto t^{-0.59\pm 0.05}$| and |$\nu _{\textrm {sa}} \propto t^{-0.86\pm 0.03}$|⁠, and that multiple shocks are contributing to the early 15.5 GHz observation. We find that the slow evolution of the reverse shock means that it contributes significantly to the low-frequency emission at all times.

To contextualize these findings, we compare our results to both thin and thick reverse shock models summarized in van der Horst et al. (2014). The distinction between thin and thick shell models refers to the depth and velocity spread of the shell that the shock is moving through. The reverse shock emission is produced as it propagates back through the shell at the front of the jet. In a thick shell scenario, the velocity spread of the ejected material is large enough such that the shock can accelerate to become relativistic, and the resulting light curves depend on the circumburst environment profile, as does the forward shock. In the thin shell scenario, the reverse shock remains Newtonian, and reverse shock light curves are dependent on the deceleration profile of the jet (Sari & Piran 1995; Mészáros & Rees 1999). With the results of our model, we cannot recreate our observations with physically realistic parameter values for either a thick or thin shell reverse shock model. We find that the reverse shock evolution that we measure is too slow compared to analytical models such as those in van der Horst et al. (2014).

Compared to the number of detailed forward shock studies, there are very few GRBs where the reverse shock is observed in sufficient detail to confidently examine certain reverse shock models. GRBs 130427A, 190114C, and 190829A are the three most well-studied GRBs with bright reverse shock components (they also happen to all have – at least tentative – very high energy components like GRB 221009A; Ackermann et al. 2014; MAGIC Collaboration 2019b; H. E. S. S. Collaboration 2021). The reverse shock component from GRB 190114C appears to match with theoretical models for reasonable physical parameters (Laskar et al. 2019). However, GRBs 130427A and 190829A could not be explained by analytical reverse shock models (van der Horst et al. 2014; Salafia et al. 2022). In the case of GRB 190829A, the best fit came from assuming a rapid decay in the magnetic field strength post-shock crossing (Salafia et al. 2022). It is possible that GRB 221009A requires a similarly complex model to explain the observed behaviour but that is beyond the scope of this work.

4.2 Forward shock

The dotted lines in Figs 4 and 5 denote the contribution from the forward shock. The forward shock component of our model dominates all of the high-frequency light curves (above 33 GHz) at all times. Moving to lower observing frequencies the forward shock contributes less, and below 10 GHz the forward shock component is always subdominant. At X-ray energies (Fig. 6), the emission is always dominated by the forward shock component (the dotted line). Given how well our model fits the X-ray data, the cooling break |$\nu _{\textrm {c}}$| seems to be situated above the X-ray regime throughout the observations. Although we do not fit our model to the optical data, we have overlaid our model on to the optical data from Fulton et al. (2023) in Fig. 7. The decay rate of our model matches that of the data except for the late time y-band data, which Fulton et al. (2023) suggested was due to a supernova component. Fig. 7 reinforces that there is significant extinction affecting the optical emission from GRB 221009A (Fulton et al. 2023; Kann et al. 2023; Levan et al. 2023; Tiengo et al. 2023). We find that nearly 2 mag of extinction are needed in the r band, decreasing to |${\sim}0.1\!-\!0.2$| mag in the y band.

Figure 6.

The X-ray light curve for GRB 221009A at 10 keV, overlaid with our best-fitting afterglow model.

Figure 7.

Optical light curves of GRB 221009A from Fulton et al. (2023), overlaid with our best-fitting afterglow model. While we do not fit our model to the optical data due to the large and mostly unconstrained extinction contribution as well as the supernova (Kann et al. 2023; Blanchard et al. 2024), our model reproduces the decay rate of the optical data well. It is clear that significant extinction, 2 mag in the r band, is needed to get the correct normalization of our model with respect to the data.

Traditional forward shock spectral models take the three frequency breaks and the peak flux density and calculate four afterglow parameters: the total kinetic energy, the circumburst density, and the fraction of kinetic energy that goes into the electrons and magnetic fields (Sari et al. 1998; Chevalier & Li 1999). From there, if a jet break is detected (an achromatic break in the light curves), the opening angle of the jet can be calculated (Sari, Piran & Halpern 1999). For GRB 221009A, we cannot calculate these parameters for two main reasons. The first is that whilst we are able to track the evolution of |$F_{\nu , \textrm {max}}$|⁠, |$\nu _{\textrm {m}}$|⁠, and |$\nu _{\textrm {sa}}$| for the forward shock, with the data we use in this work we are unable to localize |$\nu _{\textrm {c}}$| since it appears to be above the X-ray band (Williams et al. 2023), and |$\nu _{\textrm {c}}$| is needed to break the degeneracy between the different afterglow parameters. Secondly, to calculate the afterglow parameters, the observed evolution must match the model’s prediction. Otherwise, the afterglow parameters derived at each time-step will have different values.

Our model finds that |$\nu _{\textrm {m}} \propto t^{-1.06\pm 0.06}$|⁠, whereas theoretically it is expected that |$\nu _{\textrm {m}} \propto t^{-1.5}$| independent of circumburst environment density profile, strongly in disagreement with our findings. We also find that |$F_{\nu , \textrm {peak}}$| and |$\nu _{\textrm {sa}}$| do not evolve in agreement with expectations from the standard afterglow model, instead we find that |$F_{\nu , \textrm {peak}} \propto t^{-0.97\pm 0.03}$| and |$\nu _{\textrm {sa}} \propto t^{-1.4^{+0.2}_{-0.1}}$| (we note that the temporal index for |$\nu _{\textrm {sa}}$| is pushing up on the bounds set for the priors in the emcee fit). Comparatively, for a stellar wind (⁠|$k = 2$|⁠) and homogeneous (⁠|$k = 0$|⁠) environment, |$F_{\nu , \textrm {peak}}$| is expected to evolve as |$t^{-0.5}$| and |$t^{0}$| (Granot & Sari 2002), respectively, which is far slower than what we observe. The expected evolution of the synchrotron self-absorption break is also dependent on the circumburst environment’s density profile: with |$t^{0}$| and |$t^{-0.6}$| for |$k = 0$| and |$k = 2$|⁠, respectively, again the temporal indices are too slow to match our model.

Using the relations from table 5 in van der Horst et al. (2014), we can derive individual circumburst density profiles from the evolution of both |$\nu _{\textrm {sa}}$| and |$F_{\nu , \textrm {max}}$|⁠. We find that |$\nu _{\textrm {sa}} \propto t^{-1.4^{+0.2}_{-0.1}}$| and |$F_{\nu , \textrm {max}} \propto t^{-0.97\pm 0.03}$| correspond to |$k = 2.8^{+0.10}_{-0.06}$| and |$k = 2.64\pm 0.03$|⁠, respectively. Both the evolution of |$F_{\nu , \textrm {max}}$| and |$\nu _{\textrm {sa}}$| strongly favour a steeper circumburst density profile over a |$k = 2$| stellar wind profile. Such a density profile could arise from a changing mass-loss rate of the progenitor star as it reaches the end stages of its life. Standard afterglow models predict that the evolution of |$\nu _{\textrm {m}}$| is independent of the circumburst environment, therefore, we cannot assume that the slow evolution of |$\nu _{\textrm {m}}$| is due to environmental effects. In other GRBs (e.g. Bright et al. 2019), the unexpected evolution of |$\nu _{\textrm {m}}$| is considered as a result of time-varying microphysical parameters or scintillation. In the case of GRB 221009A, we find no evidence for significant scintillation effects, and time-varying microphysical parameters would cause further changes in the evolution of |$F_{\nu , \textrm {max}}$| and |$\nu _{\textrm {sa}}$|⁠, which could potentially provide an alternative explanation, other than a steep k value, for the observed behaviour.

4.2.1 Late-time evolution

Our latest observations were made with ATCA at 475 d post-burst at 5.5 and 9 GHz. Our model finds that at such late times, the forward shock is the brightest emission component during the decay phase of the light curve at these radio frequencies. Many late-time radio and X-ray light curves extending out to hundreds of days show achromatic behaviour referred to as a jet break (e.g. Tanvir et al. 2010; Kangas & Fruchter 2021). As the jet decelerates, the beaming angle, dictating the fraction of the jet that the observer can see, increases. Before the jet break, the light curve at a given frequency will decay at a shallower rate than the intrinsic evolution because a greater fraction of the jet is visible at every new time-step. At the point where the opening angle is equal to the inverse of the bulk Lorentz factor, the jet break, the whole jet is within the beaming angle, so the light curve at all wavelengths will begin to decay at a steeper rate (⁠|$t^{-3p/4}$| or |$t^{-p}$|⁠, depending on whether lateral spreading is assumed or not; Sari et al. 1999; Gao et al. 2013), which matches the intrinsic evolution of the shock. By observing the jet break, it is possible to measure the opening angle of the jet.

Jet breaks have been observed at many different times post-burst, from a fraction of a day to tens of days or even later. For most GRBs, the afterglow quickly fades below detection limits before a jet break can be observed. In some long-lasting afterglows, no jet break is observed at all for a very long time, the best example being GRB 130427A where no jet break was observed out to at least 1.9 yr post-burst (de Pasquale et al. 2016). Comparatively, we rule out the presence of a break in the light curve out to 1.3 yr based on our latest ATCA observations. We note that the presence of lateral structure, as indicated by the need for a third shock component which is discussed in Section 4.3, could disguise the jet break signature that is predicted for top-hat jets (Sari et al. 1999; Gao et al. 2013).

Whilst the presence of the jet break is used to measure the jet opening angle, the measurement is also dependent on the jet’s kinetic energy and the density of the circumburst environment. The fact that there has been no change in light-curve behaviour out to over a year post-burst due to a jet break indicates that the kinetic energy of the jet could be higher than what is deemed ‘normal’ for a regular GRB jet, the circumburst density is very low, or it has a wide jet opening angle. As already suggested by O’Connor et al. (2023), GRB 221009A may belong to a class of hyperenergetic GRBs (Chandra et al. 2008; Cenko et al. 2011; Martin-Carrillo et al. 2014), events whose kinetic energies are greater than |$10^{51}$| erg. Given the large isotropic equivalent kinetic energies inferred from modelling so far, a large jet opening angle is unlikely as it would require the beaming-corrected kinetic energy to be physically challenging, approaching that of the isotopic equivalent kinetic energy. It has been suggested (e.g. Levan et al. 2023; O’Connor et al. 2023) that a jet break occurred within the first-day post-burst. Our observations and modelling provide no evidence that such a jet break occurred.

There is also expected to be a change in the observed light-curve behaviour as the jet leaves the stellar wind bubble produced by the progenitor star and enters the surrounding homogeneous interstellar medium (ISM). The stellar wind bubble is expected to be several tens of parsecs in size (Dwarkadas 2005; Eldridge et al. 2006). For GRB 130427A, a stellar wind to homogeneous transition is ruled out to 1.9 yr post-burst. In that case, it was estimated that the jet had travelled between 50 and 105 pc, putting strong constraints on the presence/size of a termination shock, other nearby stars, etc. Our model for GRB 221009A disfavours any change in the structure of the circumburst environment out to 1.3 yr, or that the stellar wind bubble produced by the stellar progenitor exists in a very low pre-existing ISM density for the stellar wind to expand into. However, if the circumburst density profile is very steep, as our forward shock model suggests, it may be very difficult to observe such a transition. Cenko et al. (2011) suggested that the hyperenergetic events can occur in lower metallicity environments where the progenitor star maintains a higher angular momentum for longer and therefore evacuates a larger cavity with its stellar wind, therefore, delaying any change in temporal behaviour.

Studies of GRB progenitor systems predict termination shock radii to be less than 20 pc (Fryer, Rockefeller & Young 2006; Schulze et al. 2011). Using the radio source size growth rate from Giarratana et al. (2023), we estimate the distance travelled by the jet for three different assumed opening angles. For opening angles of 2, 5, and 10|$^{\circ }$|⁠, the jet should have propagated |${\sim}10$|⁠, 4, and 2 pc, respectively. At the current epoch, our observations are still consistent with the sizes of termination shocks found in the literature (Fryer et al. 2006). Therefore, we can treat these values as lower limits on the termination shock size. Continued low-frequency radio observations will be vital in tracking the jet as it continues to expand into the surrounding medium.

4.3 Extra shock

As explained in Section 3, we ran two different iterations of the third shock component in our model to test different theoretical predictions (see Table 3 for a summary, and Fig. 3 for the results). The dash–dotted lines in Figs 4 and 5 denote the contribution of this component. The most important aspect of the third spectral component is the delayed deceleration time-scale over which the component comes into the observer’s line of sight (Peng et al. 2005). We find that the deceleration time for the third component is 0.27 |$\pm$| 0.02 d, the break time in our |$F_{\nu , \textrm {max}}$| broken power-law evolution. The delayed deceleration time-scale is used to show that there is a possibility that the third component is either off-axis and therefore takes time to enter our line of sight, or that it is less relativistic than the main jet component and so needs longer to shock sufficient mass such that it undergoes significant deceleration.

To ensure that the data need the |$F_{\nu , \textrm {max}} \propto t^{3}$| rise, we also ran a separate model iteration that allows the rise index to vary (model 2 in Table 3). In this iteration, we find a broad posterior distribution, i.e. not a Gaussian posterior, extending from |$F_{\nu , \textrm {max}} \propto t^{1.6}$| to the edge of the prior which is |$F_{\nu , \textrm {max}} \propto t^{3}$|⁠. Such a broad posterior could be indicative of some lateral structure in the outflow such that the whole shock front does not enter our line of sight at once (Mooley et al. 2018b; Ryan et al. 2020).

After the peak, for a decelerating shock, afterglow models predict |$F_{\nu , \textrm {max}}$| to decay between |$t^{-1.7}$| and |$t^{-1.8}$|⁠, for |$p = 3.1$| for a stellar wind and homogeneous medium, respectively (Granot & Sari 2002). Our observations find |$F_{\nu , \textrm {max}} \propto t^{-0.71\pm 0.03}$|⁠, significantly slower than the models predict. The break frequency |$\nu _{\textrm {sa}}$| is expected to decay as |$t^{-1.1}$| and |$t^{-1.3}$|⁠, for |$k = 0$| and |$k=2$|⁠, respectively, whereas we find |$t^{-0.46^{+0.03}_{-0.04}}$|⁠. Therefore, we find that the evolution of |$\nu _{\textrm {sa}}$| for this extra component is far slower than predicted by analytical blast wave models, contrary to the evolution in the forward shock case, which is too fast.

We can also use the observed evolution to extract the density profile of the circumburst environment and p, independently of the spectral fit (van der Horst et al. 2014). In this case, we take the observed |$F_{\nu , \textrm {max}}$| and |$\nu _{\textrm {sa}}$| behaviour as a function of time and solve for p and k. However, solving for p and k does not provide physical solutions for either value, i.e. a negative value of p.

Given the clear disagreement between our modelling results using three components and expectations from analytical shock models, it is possible that this third additional component is not produced by a relativistic shock but by a slower outflow component such as a circumstellar interaction from the supernova. The peak luminosity of the extra shock is around |$10^{30}$| erg s⁻¹ Hz⁻¹ which is still an order of magnitude higher than the most luminous radio-detected supernovae (e.g. Palliyaguru et al. 2019; Dong et al. 2021), and reaches such high luminosities within a day as opposed to 100–1000 s of days later.

Therefore, we find that the origin of the additional spectral component is most likely a wider outflow or cocoon-like component, as opposed to circumstellar interaction from a supernova. Being slightly less relativistic than the jet, the cocoon will take less time to sweep up mass whose rest-mass energy is equal to that of the outflow and therefore will experience delayed deceleration. It is also likely to be slightly off-axis compared to the forward and reverse shock-emitting jet.

Cocoons have been invoked in previous GRB systems (e.g. Mooley et al. 2018a; Izzo et al. 2019) where sufficiently high-quality data have been used to infer their presence. It is possible that cocoons are a more universal component of GRBs but our observations have been too sparse to find them.

5 CONCLUSIONS

In this work, we have collated and presented the most detailed radio study of any long GRB to date. When combined with the published X-ray data, we find that the radio observations are best described with three synchrotron spectra, each evolving individually. A reverse shock component dominates the early-time low-frequency data below 20 GHz. The higher frequency radio emission and X-ray data can be ascribed to a forward shock. Due to the high temporal and spectral coverage, we are also able to constrain the evolution and properties of a third component that we attribute to a potential cocoon-like outflow. Whilst it is possible to match the different spectra with different shock components, we find that in all cases the evolution of the self-absorbed regions of the afterglow does not match up with the models currently in the literature. Also the peak frequency and peak flux show temporal behaviour that is inconsistent with theoretical afterglow models. Given the high signal-to-noise ratio of our latest observations, we aim to continue observing the afterglow of GRB 221009A for years to come to detect a potential jet break, track the jet into the non-relativistic regime, and constrain the size of the wind bubble in which this GRB resides.

ACKNOWLEDGEMENTS

The authors would like to thank the anonymous referee for their helpful comments. LR thanks Eric Burns, Courey Elliott, Geoffrey Ryan, and Ashley Villar for useful conversations in the writing of this paper. KG acknowledges support through the Australian Research Council Discovery Project DP200102243. MJM acknowledges the support of the National Science Centre, Poland through the SONATA BIS grant 2018/30/E/ST9/00208 and the Polish National Agency for Academic Exchange Bekker grant BPN/BEK/2022/1/00110. AR acknowledges funding from the NOW Aspasia grant (number: 015.016.033). RLCS acknowledges support from a Leverhulme Trust Research Project Grant.

Parts of this research were conducted by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), project numbers CE170100004 and CE230100016. We thank the Mullard Radio Astronomy Observatory staff for scheduling and carrying out the AMI-LA observations. The AMI telescope is supported by the European Research Council under grant ERC2012-StG-307215 LODESTONE, the UK Science and Technology Facilities Council, and the University of Cambridge. The Allen Telescope Array refurbishment program and its ongoing operations are being substantially funded through the Franklin Antonio Bequest. Additional contributions from Frank Levinson, Greg Papadopoulos, the Breakthrough Listen Initiative, and other private donors have been instrumental in the renewal of the ATA. Breakthrough Listen is managed by the Breakthrough Initiatives, sponsored by the Breakthrough Prize Foundation. The Paul G. Allen Family Foundation provided major support for the design and construction of the ATA, alongside contributions from Nathan Myhrvold, Xilinx Corporation, Sun Microsystems, and other private donors. The ATA has also been supported by contributions from the U.S. Naval Observatory and the U.S. National Science Foundation. The Australia Telescope Compact Array is part of the Australia Telescope National Facility⁵ that is funded by the Australian Government for operation as a National Facility managed by CSIRO. We acknowledge the Gomeroi People as the Traditional Owners of the Observatory site. This scientific work used data obtained from Inyarrimanha Ilgari Bundara/the Murchison Radio-astronomy Observatory. We acknowledge the Wajarri Yamaji People as the Traditional Owners and native title holders of the Observatory site. CSIRO’s ASKAP radio telescope is part of the Australia Telescope National Facility. Operation of ASKAP is funded by the Australian Government with support from the National Collaborative Research Infrastructure Strategy. ASKAP uses the resources of the Pawsey Supercomputing Research Centre. Establishment of ASKAP, Inyarrimanha Ilgari Bundara, the CSIRO Murchison Radio-astronomy Observatory, and the Pawsey Supercomputing Research Centre are initiatives of the Australian Government, with support from the Government of Western Australia and the Science and Industry Endowment Fund. e-MERLIN is a National Facility operated by the University of Manchester at Jodrell Bank Observatory on behalf of STFC. This paper is based (in part) on data obtained with the International LOFAR Telescope (ILT) under project code DDT20_003. LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and that are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefited from the following recent major funding sources: INSU,CNRS, Observatoire de Paris, and Université d'Orléans, France; BMBF, MIWF-NRW, and MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, the Netherlands; the Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland. Based on observations carried out under project numbers S22BC, S22BF, W22BU, and S22BE with the IRAM NOEMA Interferometer. IRAM is supported by INSU,CNRS (France), MPG (Germany), and IGN (Spain). We thank the staff of the GMRT that made these observations possible. GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under a cooperative agreement by Associated Universities, Inc.

DATA AVAILABILITY

All the new data presented in this paper are given as the supplementary material.

Footnotes

1

https://research.csiro.au/casda/

2

https://linc.readthedocs.io/en/latest/index.html

3

Second data release version: https://github.com/mhardcastle/ ddf-pipeline. The tier1-july2018.cfg pipeline configuration was used.

4

https://www.iram.fr/IRAMFR/GILDAS

5

https://ror.org/05qajvd42