We present radio observations of the longduration gammaray 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 postburst 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 wellstudied GRB radio afterglows to date. Our observations are best explained by three synchrotronemitting 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 latetime spectral or temporal changes that could result from a jet break but note that any lateral structure could significantly affect a jet break signature.
['Astrophysics - High Energy Astrophysical Phenomena']
Rocking the BOAT the ups and downs of the longterm radio light curve for GRB 221009A
{'Rocking the BOAT: the ups and downs of the long-term radio light curve for GRB 221009A': "- L. Rhodes, 1 ★ A. J. van der Horst, 2 , J. S. Bright 1 , J. K. Leung, 3 , 4 , 5 , G. E. Anderson 6 , R. Fender 1 , 7 ,\n- J. F. Agüí Fernandez 8 , M. Bremer, 9 P. Chandra 10 , D. Dobie 11 , 12 , W. Farah 13 , 14 , S. Giarratana, 15 , K. Gourdji 16 ,\n- D. A. Green 17 , E. Lenc 18 , M. J. Michałowski 19 ,T. Murphy 11 , 12 , A. J. Nayana 20 , A. W. Pollak 13 , 14 ,\n- A. Rowlinson 21 , 22 , F. Schussler 23 , A. Siemion 1 , 13 , 14 , 24 , 25 R. L. C. Starling 26 , P. Scott 17 , C. C. Thöne 27 ,\n- D. Titterington 17 A. de Ugarte Postigo 9 , 28\n- 1 Astrophysics, Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford, OX1 3RH, UK\n- 2 Department of Physics, George Washington University, 725 21st St NW, Washington, DC, 20052, USA\n- 3 David A. Dunlap Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada\n- 4 Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada\n- 5 Racah Institute of Physics, The Hebrew University of Jerusalem, Jerusalem 91904, Israel\n- 6 International Centre for Radio Astronomy Research, Curtin University, GPO Box U1987, Perth, WA 6845, Australia\n- 7 Department of Astronomy, University of Cape Town, Private Bag X3, Rondebosch 7701, South Africa\n- 8 Centro Astronómico Hispano en Andalucía, Observatorio de Calar Alto, Sierra de los Filabres, Gérgal, Almería, 04550, Spain\n- 9 Observatoire de la Côte d'Azur, Université Côte d'Azur, Boulevard de l'Observatoire, 06304 Nice, France\n- 9 Institut de Radioastronomie Millimétrique, 300 Rue de la Piscine, 38400 Saint-Martin-d'Hères, France\n- 10 National Radio Astronomy Observatory, 520 Edgemont Rd, Charlottesville VA 22903, USA\n- 11 Sydney Institute for Astronomy, School of Physics, The University of Sydney, New South Wales 2006, Australia\n- 12 ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Hawthorn, Victoria, Australia\n- 13 SETI Institute, 339 Bernardo Ave, Suite 200 Mountain View, CA 94043, USA\n- 14 Berkeley SETI Research Centre, University of California, Berkeley, CA 94720, USA\n- 15 Istituto Nazionale di Astrofisica, Osservatorio Astronomico di Brera, via E. Bianchi 46, 23807 Merate, (LC), Italy\n- 16 Centre for Astrophysics and Supercomputing, Swinburne University of Technology, P.O. Box 218, Hawthorn, VIC 3122, Australia\n- 17 Astrophysics Group, Cavendish Laboratory, 19 J.J. Thomson Avenue, Cambridge, CB3 0HE, UK\n- 18 Australian Telescope National Facility, CSIRO Astronomy and Space Science, PO Box 76, Epping, NSW 1710, Australia  \n19  \nAstronomical Observatory Institute, Faculty of Physics, Adam Mickiewicz University, ul. Słoneczna 36, 60-286 Poznań, Poland  \n- 20 Department of Astronomy, University of California, Berkeley, USA\n- 21 Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH, Amsterdam, the Netherlands\n- 22 ASTRON, the Netherlands Institute for Radio Astronomy, Oude Hoogeveensedĳk 4, 7991 PD Dwingeloo, the Netherlands\n- 23 IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France\n- 24 Department of Physics and Astronomy, University of Manchester, UK\n- 25 University of Malta, Institute of Space Sciences and Astronomy, Msida, MSD2080, Malta\n- 26 University of Leicester, School of Physics and Astronomy, University Road, Leicester LE1 7RH\n- 27 Astronomical Institute of the Czech Academy of Sciences (ASU-CAS), Fricova 298, Ondřejov, 251 65, Czech Republic\n- 28 Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France  \nAccepted XXX. Received YYY; in original form ZZZ", 'ABSTRACT': 'Wepresent radio observations of the long-duration gamma-ray burst (GRB) 221009A which has become known to the community as the Brightest Of All Time or the BOAT. Our observations span the first 475 days 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 which 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.  \nKey words: 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 seconds, 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-HXMT and GECAM-C (An et al. 2023), Konus Wind and SRG/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 × 10 55 erg, 1 . 5 × 10 55 erg and 1 . 2 × 10 55 erg by Lesage et al. (2023); Anetal. (2023); Frederiks et al. (2023, between 1-10,000 keV, 10 keV - 6 MeV and 20 keV - 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 year event (Burns et al. 2023). In fact, GRB 221009A was so bright that the prompt emission caused disturbances in the ionosphere (Hayes & Gallagher 2022).  \nThe afterglow to GRB 221009A has been detected consistently between 0.4 GHz and 20 TeV (Laskar et al. 2023; LHAASO Collaboration et al. 2023). In terms of spectral coverage, it exceeds the all other TeV afterglows with radio detections like GRB 190114C or GRB 190829A (MAGIC Collaboration et al. 2019a; H. E. S. S. Collaboration et al. 2021). In terms of data quantity and quality, it exceeds the GHz-to-GeV afterglow of GRB 130427A (e.g. van der Horst et al. 2014; Levan et al. 2014; Ackermann 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 & Rhodes 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 timescales of weeks to months (Frail et al. 2000; van der Horst et al. 2005; Resmi et al. 2005). The origin of the early time radio peaks are 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 et al. 1998; Wĳers & Galama 1999).  \nThe 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 et al. 2019a; H. E. S. S. Collaboration et al. 2021; Abe et al. 2024), the radio observations provide the necessary context for understanding the physics of the jetted GRB outflow, together with multi-wavelengths 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 do not seem to follow expectations of the standard model that is typically used to describe radio afterglows (e.g. Wĳers & Galama 1999; Granot & Sari 2002). The extremely dense sampling of the light curves at various radio frequencies as presented in this paper is unprecedented  \nand allows for detailed modelling that will lead to better descriptions of GRB jets and the relevant emission processes.  \nThe 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 which 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 et al. 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 et al. 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 tophat behaviour (van der Horst et al. 2014; Bright et al. 2019; Rhodes et al. 2022).  \nBesides these macrophysical considerations, high-quality multiwavelength data as presented here reveals 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 multi-wavelength 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).  \nGiven 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 𝐹 𝜈 ∝ 𝑡 𝛼 𝜈 𝛽 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 summarised in Table 1. In addition to the datasets we present here, our work also incorporates the previously published radio data from Laskar et al. (2023); Giarratana et al.  \n(2023) and Bright & Rhodes et al. (2023), and the X-ray data from Williams et al. (2023).', '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 & Rhodes 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 days post-burst when the first non-detection occurred. Between 210 and 320 days post-burst, we concatenated separate non-detections to obtain deeper limits.  \nAMI-LA data is reduced using a custom software package: /r.pc/e.pc/d.pc/u.pc/c.pc/e.pc\\_/d.pc/c.pc (Perrott et al. 2013). The software performs bandpass and flux scaling using 3C286 and complex gain calibration using J1925+2106. Flagging and imaging is done in /c.pc/a.pc/s.pc/a.pc 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 & Rhodes 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 221009AfieldwiththeAustralianSquare 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 figure 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 1 ).  \nObservations 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, fullpolarisation 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 ∼ 50 𝜇 Jy beam -1 . We evaluated and corrected for the systematic fluxscale 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 ∼ 200 miles north of San Fransisco, the Allen Telescope Array is a 42-element radio interferometer hosted at the Hat Creek Radio Observatory. Mounted on the focus of each element is a dualpolarization, log-periodic feed that is cryogenically cooled and sensitive to radiation in the range of 1 to 12 GHz. Analogue signals from the array are transmitted through fibre to a centralised signal processing room and are split into 4 independent chains that get multiplexed by 4 tunable local oscillators in a super-heterodyne system. The current correlator backend supports the digitisation of 2 out of  \nthe 4 available tunings for 20 of the 42 antennas, where each tuning can be placed anywhere in the available RF range of the log periodic feed, with ∼ 700 MHz of usable bandwidth for each.  \nThe radio counterpart of GRB 221009A was observed extensively with the ATA beginning just a few hours after the burst as reported in Bright & Rhodes et al. (2023). Here we build on that work and utilised the flexible frequency tunability of the ATA to monitor the 1-10 GHz spectral evolution over its entire outburst. Either 3C147, 3C48, or 3C286 was observed as flux calibrator at the beginning of each observing block, and a 10 minute observation of the phase calibrator J1925 + 2106 was interleaved for every 30 minutes 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': "Wecarried 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), C3542 (PI: G. E. Anderson). These observations were carried out using the 5.5/9, 16.7/21.2, 33/35, and 43/45 GHzreceiver configurations, with a bandwidth of 2048 MHz for each intermediate frequency.  \nFor each observation, we reduced the visibility data using standard procedures in M/i.pc/r.pc/i.pc/a.pc/d.pc (Sault et al. 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, we split the 2048 MHzbandwidth into further sub-bands 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 M/i.pc/r.pc/i.pc/a.pc/d.pc 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  \nspread across the UK. With a maximum baseline of 217 km, whilst observing at 5 GHz, it can resolve angular scales of 0.05'. 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-band. Our L- and C-band observations were centered 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 & Rhodes et al. (2023).  \nAll observations were reduced using the e -MERLIN pipeline within /c.pc/a.pc/s.pc/a.pc (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 OQ208 and J1905+1943, respectively, along with flux scaling using 3C286. Further flagging of the target field is conducted. We performed interactive cleaning and deconvolution using the casa task tclean .", '2.6 LOFAR': "Eight hours of Director's Discretionary Time with the Low Frequency Array (LOFAR; DDT20\\_003) were awarded to observe GRB 221009A. The allocated time was split into two observing runs of 4-hours, which took place on 18 and 20 July 2023 at matching local sidereal times. Each observing run was preceded by a 10-minute calibrator scan of 3C295. 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 second. 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 preprocessing. Both target observations were calibrated for direction independent effects using LINC 2 with default settings, a pipeline developed by the observatory to correct for various instrumental and ionospheric effects present in interferometric LOFAR data (de Gasperin et al. 2019; van Weeren et al. 2016; Williams et al. 2016). Due to its relative proximity, Cygnus A was subtracted from the visibilities using the 'demixing' step in LINC. The data were further averaged to 4 channels of 48.82 kHz per subband and 4 seconds during calibration. The resulting calibrated data were concatenated into groups of 20 subbands and averaged in time to 8 seconds. These data products from both observations were subsequently jointly put through /d.pc/d.pc/f.pc/p.pc/i.pc/p.pc/e.pc/l.pc/i.pc/n.pc/e.pc 3 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 arcseconds and 1.5 arcseconds pixel resolution.", '2.7 NOEMA': 'The NOrthern Extended Millimetre Array (NOEMA, situated in the southern French Alps) monitored GRB 221009A between October 10 th 2022 and April 25 th 2023 in the 3, 2 and 1mm bands. Interferometer configurations were medium-extended C and extended A  \n2  \nhttps://linc.readthedocs.io/en/latest/index.html  \nconfigurations with up to 12 antennas, primary flux calibrators were MWC349 and LKHA101. The data were reduced with the /c.pc/l.pc/i.pc/c.pc and /m.pc/a.pc/p.pc/p.pc/i.pc/n.pc/g.pc software packages that are part of the /g.pc/i.pc/l.pc/d.pc/a.pc/s.pc 4 package. Fluxes and their errors were derived from point-source 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 January 2023 and then in March 2023. We recorded the data in 2048 frequency channels covering a bandwidth of 400 MHz with an integration time of ∼ 10 s. We used 3C286 and 3C48 as flux density and bandpass calibrators. J1924+3329 was used as a phase calibrator.  \nThe data were analysed using the /c.pc/a.pc/s.pc/a.pc 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 & MODEL': "There have been several GRB 221009A afterglow modelling efforts whichhaveusedasubsetoftheradiodata published to date (including but not limited to Laskar et al. 2023; O'Connor et al. 2023; Levan et al. 2023; Gill & Granot 2023). Here, we present the results of our observing campaigns and described out modelling of the radio and X-ray afterglow.", '3.1 Light Curves and SEDs': "Theradio data presented in this paper spans three orders of magnitude in frequency space, from 0.15 MHz to 230 GHz, and lasts out to 475 days post-burst. Figure 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.  \nAbove 19 GHz, the afterglow is decaying at all times, with observations obtained between 1 and 200 days post-burst (the top two rows of Figure 1). The light curves between 90 and 105 GHz in Figure 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 & Rhodes 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 days. 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.  \nFigure 2 shows the broadband radio spectral energy distributions (SEDs) throughout our campaign. For the first 30 days, a lowfrequency turnover is visible and the below-turnover spectral index is consistent with 𝛽 ∼ 5 / 2 below the turnover. Above the turnover, we  \nTable 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 𝜎 upper limit in the uncertainty column. The full list of radio observations are presented in supplementary material online.  \nfindaflatspectrumextendingtothehighestfrequencies ( ∼ 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 days post-burst does the spectrum steepen with typical optically thin spectral indices ( 𝛽 ∼ -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 𝜈 c , sits in the XRT band. The broken power law is favoured but the fits are only performed on data up to one day post-burst whereas the X-ray light curve itself extends out to 200 days post-burst.", '3.2 Modelling': 'Here, we build on previous modelling efforts by combining our new observations from AMI-LA, ATA, ATCA,ASKAP, e -MERLIN,LOFAR,NOEMAanduGMRTwithradiodataavailableintheliterature. 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 VHE component whose origin and emission mechanism is still debated (Aharonian et al. 2023; LHAASO Collaboration et al. 2023; Savchenko et al. 2024).  \nWe 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 (Figure 1) and SEDs (Figure 2). Each synchrotron spectrum is constructed of four power-law slopes divided by three frequency breaks: the synchrotron self-absorption break ( 𝜈 sa ), the characteristic or minimum electron energy break ( 𝜈 m ), and the cooling break ( 𝜈 c , above which radiative cooling is important). The peak of the spectrum, 𝐹 𝜈, max , is at whichever frequency break of 𝜈 sa or 𝜈 m is higher. The spectral index of each branch depends on the order of the frequency breaks. In the regime where 𝜈 sa < 𝜈 m < 𝜈 c , the spectral indices are 𝐹 𝜈<𝜈 sa ∝ 𝜈 2 , 𝐹 𝜈 sa <𝜈<𝜈 m ∝ 𝜈 1 / 3 , 𝐹 𝜈 m <𝜈<𝜈 c ∝ 𝜈 ( 1 -𝑝 )/ 2 and 𝐹 𝜈 c <𝜈 ∝ 𝜈 -𝑝 / 2 , where 𝑝 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 et al. 2013). In the regime where 𝜈 m < 𝜈 sa < 𝜈 c , the spectral indices are 𝐹 𝜈<𝜈 m ∝ 𝜈 2 , 𝐹 𝜈 m <𝜈<𝜈 sa ∝ 𝜈 5 / 2 , 𝐹 𝜈 sa <𝜈<𝜈 c ∝ 𝜈 ( 1 -𝑝 )/ 2 and 𝐹 𝜈 c <𝜈 ∝ 𝜈 -𝑝 / 2 . As the jet expands and evolves, the spectral  \nbreaks 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, 𝜌 ∝ 𝑟 -𝑘 , where 𝑘 = 0 for a homogeneous medium and 𝑘 = 2 represents a stellar wind (Granot & Sari 2002; Granot & van der Horst 2014).  \nWe use /e.pc/m.pc/c.pc/e.pc/e.pc to fit our respective models to the data (ForemanMackey et al. 2013). Each model uses 40 walkers and runs for at least 70000 steps or until convergence. All priors are uniform, and the only priors with fixed bounds were 𝑝 ∈ [ 1 . 5 , 3 . 5 ] to help rule out unphysical solutions. The best fit value for each parameter is the 50 th percentile post burn in of the posterior distribution, and the 84 th and 16 th 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 & Rhodes et al. (2023), we find that the peak of the synchrotron spectrum is produced by 𝜈 sa and fit for the normalisation and evolution of the spectrum as well as 𝑝 . The second component is a forward shock that appears to dominate the optical and X-rays (e.g. Williams et al. 2023; Fulton et al. 2023; Shrestha et al. 2023), and also the late-time radio emission. Here we allow both 𝜈 sa and 𝜈 m to vary freely. We fit for the normalisation and evolution of 𝐹 𝜈, max , 𝜈 sa and 𝜈 m as well as 𝑝 . The resulting model parameters are provided in Table 2.  \nWe find that the two-component model cannot reproduce the flat spectrum observed shown in Figure 2, the posterior distribution of 𝑝 for the reverse shock always ends up at the lower bound of the prior with values for 𝑝 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 which was needed in the two-component model to explain the flat spectrum that is present during the first ∼ 150 days (Figure 2) and the additional bumps in the 5 and 15.5 GHz light curves around 5-10 days post-burst (see Figure 1).  \nTo best explore the parameter space of the third component, first, we test both 𝜈 sa and 𝜈 m as the peak frequency of the third component and find that 𝜈 sa provides a better fit. Then we consider two different iterations of this extra shock with differing degrees of freedom, which are summarised 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  \nFigure 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 downwards-facing triangles for 3 𝜎 upper limits.  \n<!-- image -->  \nFigure 2. Broadband radio SEDs for GRB 221009A as a function of time. As in Figure 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 Figure 1 are also shown here.  \n<!-- image -->  \nTable 2. The parameter values (50 th percentile) and their associated uncertainties (18 th and 64 th percentiles) derived for our best-fit two-component model. Any 𝛼 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 𝜈, max and 𝜈 sa are normalised to 1 day and 6.5 days, respectively. For each shock, p is the value of the electron energy spectral index.  \ncomponents follows a smoothly broken power law (Rhodes et al. 2020):  \n𝐹 𝜈 = 𝐹 𝜈, max GLYPH<18> 0 . 5 GLYPH<18> 𝑡 𝑡 𝑏 GLYPH<19> -𝛼 1 𝑠 + 0 . 5 GLYPH<18> 𝑡 𝑡 𝑏 GLYPH<19> -𝛼 2 𝑠 GLYPH<19> -1 𝑠 (1)  \nwhere 𝐹 𝜈, max is the flux density at the break time 𝑡 𝑏 , 𝛼 1 and 𝛼 2 are the power-law indices, and 𝑠 is the smoothing parameter which we set to be 0.5. In Model 1, the synchrotron self-absorption break follows a single power law: 𝜈 sa = 𝜈 sa , 0 𝑡 𝛼 sa where 𝜈 sa , 0 is the location of the self-absorption break at 1 day post-burst. We set 𝛼 𝐹 𝜈, max , 1 = 3 and both 𝛼 𝐹 𝜈, max , 2 (defined in Table 3 as 𝛼 ) and 𝛼 sa can vary freely. We invoke a 𝛼 𝐹 𝜈, max , 1 = 3 as done in Peng et al. (2005) which is used in the regime where a blastwave that is initially off-axis has undergone significant deceleration and so the radiation begins to enter the observers' line of sight. In the 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.  \nIn model 2, both the peak flux density and 𝜈 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.  \nFigure 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 𝑝 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. Figures 4 and 5 show our best-  \nTable 3. Summary of the different iterations of the three-component model which explore the possible evolution of the third shock component. Each 𝛼 corresponds to a temporal index of the subscripted value, e.g., 𝛼 𝐹 𝜈, max , 1 corresponds to the first temporal index used to describe the behaviour of 𝐹 𝜈, 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 Figures 4 and 5. The best fit parameters are shown in Table 4.  \nTable 4. The parameter values (50 th percentile) and their associated uncertainties (18 th and 64 th percentiles) derived for our best-fit three-component model (model 1). Any 𝛼 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 𝜈, max and 𝜈 sa are normalised to 1 day, 6.5 days and t dec , respectively, where t dec is a parameter we fitted for. For each shock, p is the value of the electron energy spectral index.  \nting model overlaid on the light curves and SEDs. Figure 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 over predicts the late time 0.8 GHz flux. Figure 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 𝜈 sa much lower than the observed position.  \nFigure 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 (50 th percentile value) of the posterior distribution for clarity. The lefthand vertical axis of each plot corresponds to the evolution of the frequency breaks (dotted and dashed lines for 𝜈 sa and 𝜈 m , respectively). The righthand vertical axis shows the evolution of the peak flux (solid lines) of each shock component.  \n<!-- image -->", '4 DISCUSSION': 'Here we discuss the implications of our best-fitting three-component model and place them in context of other detailed radio studies of GRBs.', '4.1 Reverse shock': 'The dashed lines in Figures 4 and 5 show the contribution of the reverse shock from our model. Bright & Rhodes et al. (2023) used radio observations in the first five days post-burst to measure the evolution of F max and 𝜈 sa with time. They found that 𝐹 𝜈, max ∝ 𝑡 -0 . 70 ± 0 . 02 and 𝜈 sa ∝ 𝑡 -1 . 08 ± 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 max ∝ 𝑡 -0 . 59 ± 0 . 05 and 𝜈 sa ∝ 𝑡 -0 . 86 ± 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.  \nTo contextualise these findings, we compare our results to both thin and thick reverse shock models summarised 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).  \nCompared to the number of detailed forward shock studies, there are very few GRBs where the reverse shock is observed in sufficient  \ndetail 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 GRB221009A; Ackermann et al. 2014; MAGIC Collaboration et al. 2019b; H. E. S. S. Collaboration et al. 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 Figures 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 (Figure 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 𝜈 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 onto the optical data from Fulton et al. (2023) in Figure 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. Figure 7 reinforces that there is significant extinction affecting the optical emission from GRB 221009A (Fulton et al. 2023; Levan et al. 2023; Kann et al. 2023; Tiengo et al. 2023). We find that nearly 2 magnitudes of extinction are needed in the r-band, decreasing to ∼ 0 . 1 -0 . 2 magnitudes in the y-band.  \nTraditional forward shock spectral models take the three frequency breaks and the peak flux density and calculate four afterglow param-  \n2  \nFigure 4. The multi-frequency radio light curves for GRB 221009A overlaid with our best-fit three-component model (model 1).  \n<!-- image -->  \nFigure 5. Broadband radio SEDs for GRB 221009A as a function of time with our best-fit three-component model (model 1) overlaid.  \n<!-- image -->  \nFigure 6. The X-ray light curve for GRB 221009A at 10 keV, overlaid with our best-fit afterglow model.  \n<!-- image -->  \nFigure 7. Optical light curves of GRB 221009A from Fulton et al. (2023), overlaid with our best-fit 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 magnitudes in the r-band, is needed to get the correct normalisation of our model with respect to the data.  \n<!-- image -->  \neters: the total kinetic energy, the circumburst density, 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 et al. 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 𝐹 𝜈, max , 𝜈 m and 𝜈 sa  \nfor the forward shock, with the data we use in this work we are unable to localise 𝜈 c since it appears to be above the X-ray band (Williams et al. 2023), and 𝜈 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.  \nOur model finds that 𝜈 m ∝ 𝑡 -1 . 06 ± 0 . 06 , whereas theoretically it is expected that 𝜈 m ∝ 𝑡 -1 . 5 independent of circumburst environment density profile, strongly in disagreement with our findings. We also find that 𝐹 𝜈, peak and 𝜈 sa do not evolve in agreement with expectations from the standard afterglow model, instead we find that 𝐹 𝜈, peak ∝ 𝑡 -0 . 97 ± 0 . 03 and 𝜈 sa ∝ 𝑡 -1 . 4 + 0 . 2 -0 . 1 (we note that the temporal index for 𝜈 sa is pushing up on the bounds set for the priors in the /e.pc/m.pc/c.pc/e.pc/e.pc fit). Comparatively, for a stellar wind ( 𝑘 = 2) and homogeneous ( 𝑘 = 0) environment, 𝐹 𝜈, peak is expected to evolve as 𝑡 -0 . 5 and 𝑡 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 𝑡 0 and 𝑡 -0 . 6 for 𝑘 = 0 and 𝑘 = 2, respectively, again the temporal indices are too slow to match our model.  \nUsing the relations from table 5 in van der Horst et al. (2014), we can derive individual circumburst density profiles from the evolution of both 𝜈 sa and F 𝜈, max . We find that 𝜈 sa ∝ 𝑡 -1 . 4 + 0 . 2 -0 . 1 and F 𝜈, max ∝ 𝑡 -0 . 97 ± 0 . 03 corresponds to 𝑘 = 2 . 8 + 0 . 10 -0 . 06 and 𝑘 = 2 . 64 ± 0 . 03, respectively. Both the evolution of 𝐹 𝜈, max and 𝜈 sa strongly favour a steeper circumburst density profile over a 𝑘 = 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 𝜈 m is independent of the circumburst environment, therefore we cannot assume that the slow evolution of 𝜈 m is due to environmental effects. In other GRBs, (e.g. Bright et al. 2019) the unexpected evolution of 𝜈 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 𝐹 𝜈, max and 𝜈 sa , which could potentially provide an alternative explanation, other than a steep 𝑘 value, for the observed behaviour.", '4.2.1 Late-time evolution': "Our latest observations were made with ATCA at 475 days 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 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 ( 𝑡 -3 𝑝 / 4 or 𝑡 -𝑝 , 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.  \nJet breaks have been observed at many different times post-burst,  \nfrom 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 years post-burst (De Pasquale et al. 2016). Comparatively, we rule out the presence of a break in the light curve out to 1.3 years 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 which is predicted for top hat jets (Sari et al. 1999; Gao et al. 2013).  \nWhilst 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 a it has a wide jet opening angle. As already suggested by O'Connor et al. (2023), GRB 221009A may belong to a class of hyper-energetic 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.  \nThere 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. 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 years post-burst. In that case, it was estimated that the jet had travelled between 50 and 105 parsecs, 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 years, or that the stellar wind bubble produced by the stellar progenitor exists in a very low preexisting 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 hyper-energetic 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.  \nStudies of GRB progenitor systems predict termination shock radii to be less than 20 parsec (Fryer et al. 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 · , the jet should have propagated ∼ 10, 4 and 2 parsecs, 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 Figure 3 for the results). The dash-dotted lines in Figures 4 and 5 denote the contribution of this component. The most important aspect of the third spectral component is the delayed deceleration timescale 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 ± 0.02 days, the break time in our 𝐹 𝜈, max broken power law evolution. The delayed deceleration timescale 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.  \nTo ensure that the data needs the 𝐹 𝜈, max ∝ 𝑡 3 rise, we also ran a separate model iteration which 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 𝐹 𝜈, max ∝ 𝑡 1 . 6 to the edge of the prior which is 𝐹 𝜈, max ∝ 𝑡 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).  \nAfter the peak, for a decelerating shock, afterglow models predict 𝐹 𝜈, max to decay between 𝑡 -1 . 7 and 𝑡 -1 . 8 , for 𝑝 = 3 . 1for a stellar wind and homogeneous medium, respectively (Granot & Sari 2002). Our observations find 𝐹 𝜈, max ∝ 𝑡 -0 . 71 ± 0 . 03 , significantly slower than the models predict. The break frequency 𝜈 sa is expected to decay as 𝑡 -1 . 1 and 𝑡 -1 . 3 , for 𝑘 = 0 and 2, respectively, whereas we find 𝑡 -0 . 46 + 0 . 03 -0 . 04 . Therefore, we find that the evolution of 𝜈 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.  \nWecanalsousetheobservedevolution to extract the density profile of the circumburst environment and 𝑝 , independently of the spectral fit (van der Horst et al. 2014). In this case, we take the observed 𝐹 𝜈, max and 𝜈 sa behaviour as a function of time and solve for 𝑝 and 𝑘 . However, solving for 𝑝 and 𝑘 does not provide physical solutions for either value, i.e., a negative value of 𝑝 .  \nGiven 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 a slower outflow component such a circumstellar interaction from the supernova. The peak luminosity of the extra shock is around 10 30 erg s -1 Hz -1 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.  \nTherefore, 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.  \nCocoons have been invoked in previous GRB systems (e.g. Mooley et al. 2018a; Izzo et al. 2019) where sufficiently high-quality data has 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 which 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 behavior that is inconsistent with theoretical afterglow models. Given the high signal-to-noise 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 the thank the anonymous referee for their helpful comments. L.R. thanks Eric Burns, Courey Elliott, Geoffrey Ryan and Ashley Villar for useful conversations in the writing of this paper. K.G. acknowledges support through the Australian Research Council Discovery Project DP200102243. M.J.M. 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.  \nParts 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 ERC2012StG-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 US Naval Observatory and the US National Science Foundation. The Australia Telescope Compact Array is part of the Australia Telescope National Facility 5 which is funded by the Australian Government for operation as a National Facility managed by CSIRO. We acknowledge the Gomeroi  \npeople as the Traditional Owners of the Observatory site. This scientific work uses 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 5 . 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 MurchisonRadio-astronomyObservatoryandthePawseySupercomputing 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, that 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: CNRS-INSU, Observatoire de Paris and Université d'Orléans, France; BMBF, MIWF-NRW, 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 cooperative agreement by Associated Universities, Inc.", 'DATA AVAILABILITY': 'All the new data presented in this paper is given the supplementary material.', 'REFERENCES': '```\nAbe H., et al., 2024, MNRAS, 527, 5856 Achterberg A., Gallant Y. A., Kirk J. G., Guthmann A. W., 2001, MNRAS, 328, 393 Ackermann M., et al., 2014, Science, 343, 42 Aharonian F., et al., 2023, ApJ, 946, L27 An Z.-H., et al., 2023, arXiv e-prints, p. arXiv:2303.01203 Anderson G. E., et al., 2014, MNRAS, 440, 2059 Blanchard P. K., et al., 2024, Nature Astronomy, Bright J. S. & Rhodes L. et al., 2023, Nature Astronomy, 7, 986 Bright J. S., et al., 2019, MNRAS, 486, 2721 Burns E., et al., 2023, ApJ, 946, L31 Cenko S. B., et al., 2011, ApJ, 732, 29 Chandra P., et al., 2008, ApJ, 683, 924 Chevalier R. A., Li Z.-Y., 1999, ApJ, 520, L29 Clark B. G., 1980, A&A, 89, 377\n```'}