L Rhodes, A J van der Horst, R Fender, D R Aguilera-Dena, J S Bright, S Vergani, D R A Williams, Jet-cocoon geometry in the optically dark, very high energy gamma-ray burst 201216C, Monthly Notices of the Royal Astronomical Society, Volume 513, Issue 2, June 2022, Pages 1895–1909, https://doi.org/10.1093/mnras/stac1057

We present the results of a radio observing campaign on GRB 201216C, combined with publicly available optical and X-ray data. The detection of very high energy (VHE, >100 GeV) emission by MAGIC makes this the fifth VHE GRB at the time of publication. Comparison between the optical and X-ray light curves show that GRB 201216C is a dark GRB, i.e. the optical emission is significantly absorbed and is fainter than expected from the X-ray detections. Our e-MERLIN data also shows evidence of diffractive interstellar scintillation. We can study the column density along the line of sight to the GRB in both the host galaxy, from the damped optical light curve, and the Milky Way, via scintillation studies. We find that the afterglow is best modelled using a jet-cocoon geometry within a stellar wind environment. Fitting the data with a multicomponent model, we estimate that the optical, X-ray, and higher frequency radio data before ∼25 d originates from an ultrarelativistic jet with an isotropic equivalent kinetic energy of (0.6–10) × 1052 erg and an opening angle of ∼1–9°. The lower frequency radio emission detected by MeerKAT, from day 28 onwards, is produced by the cocoon with a kinetic energy that is between two and seven orders of magnitude lower (0.02–50) × 1048 erg. The energies of the two components are comparable to those derived in simulations of such scenarios.

Long Gamma-ray Bursts (lGRBs) are flashes of gamma radiation usually lasting upwards of two seconds (Kouveliotou et al. 1993). At cosmological distances, they are among the most energetic transients known. lGRBs are thought to be produced via internal processes in jets launched during the core collapse of a subpopulation of fast rotating massive stars (see Levan et al. 2016, for a recent review). Evidence for this connection has come from the presence of Type Ic broad-line supernovae signatures in the optical afterglow spectra of lGRBs (e.g. Hjorth et al. 2003).

Following the lGRB prompt emission, there is a broad-band afterglow, often visible from radio frequencies to X-ray energies. The afterglow is produced when the jet interacts with the circumburst medium, creating two shocks: a forward and a reverse shock. The forward shock propagates into the surrounding medium whereas the reverse shock travels back into the jet and towards the newly formed compact object. The forward shock accelerates the electrons in the circumburst medium into a power-law distribution in energy: N(E)dE ∝ E−pdE, which subsequently cools emitting synchrotron and inverse-Compton radiation. The synchrotron emission is interpreted in terms of the fireball model, the standard model used for GRB afterglows (Rees & Meszaros 1992).

The synchrotron spectrum is described using four parameters: three break frequencies (Sari, Piran & Narayan 1998) i.e. the self-absorption frequency (νSA), the frequency corresponding to the electrons with the minimum energy (νm), and the cooling frequency (νc), and |${\rm F}_{\nu ,\, {\rm max}}$|⁠ : the flux density at the peak of the spectrum, the higher of νm or νSA. The three break frequencies are connected by power laws of flux as a function of frequency. The aforementioned four parameters evolve with time and they are dependent on the circumburst environment, its density and radial profile (n or A*), the isotropic equivalent kinetic energy of the jet (EK,ISO), and the jet microphysical parameters, i.e. the fraction of the shock energy given to the electrons and magnetic fields (ϵe and ϵB, respectively; Chevalier & Li 1999; Wijers & Galama 1999; Granot & Sari 2002).

The reverse shock usually dominates between radio to optical wavelengths at early times, often decaying quickly within hours and days in the optical and radio wavebands, respectively (e.g. Kulkarni et al. 1999; Chandra & Frail 2012; Huang et al. 2016; Laskar et al. 2016; Alexander et al. 2017). The forward shock emission is visible at X-ray energies from minutes to a few hours post burst (e.g. Racusin et al. 2009; Oates et al. 2011). After the reverse shock fades, the forward shock becomes visible at lower frequencies and can be visible for up to hundreds of days (Van der Horst et al. 2008).

As the peak of the synchrotron spectrum moves to lower frequencies with time, the afterglow light curves evolve chromatically as the frequency breaks pass through different observing bands at different times. At later times, some light curves also show achromatic breaks, caused by changes due to the jet’s geometry and relativistic beaming of the emission. This is called a jet break and occurs in the regime where Γ < 1/θj, where Γ is bulk Lorentz factor of the forward shock and θj is the opening angle of the jet. The deceleration of the jet is such that the beaming cone widens allowing the observer to see the edges of the jet. As a result the observed flux of the jet begins to decay rapidly. Large sample studies of lGRB jet opening angles show jets to be highly collimated with θj = 7|$^{\textrm {+11}}_{\textrm {-4}}{}^{\circ }$| on average (Laskar et al. 2014).

Observations of some GRB afterglows have shown evidence of a second forward shock component originating from a wider outflow (Resmi et al. 2005; Starling et al. 2005; Racusin et al. 2008; Kamble et al. 2009; Filgas et al. 2011; Van der Horst et al. 2014; Lan, Wu & Dai 2018; Chen et al. 2020). This second component is sometimes interpreted as a cocoon (e.g. Chen et al. 2020). Cocoons are often suggested to explain why GRB jets are so highly collimated. Magnetohydrodynamical numerical simulations show that as the relativistic jet propagates through the collapsing star, it deposits a large amount of energy into the surrounding material forming a cocoon. In turn, the cocoon reduces the lateral expansion of the jet resulting in a high degree of collimation when the jet breaks free (Ramirez-Ruiz, Celotti & Rees 2002; Zhang, Woosley & Heger 2004; De Colle, Kumar & Aguilera-Dena 2018). The kinetic energy of the cocoon is expected to be several orders of magnitude less than the jet, more similar to relativistic supernovae, with a bulk Lorentz factor lower than the core of the jet (De Colle et al. 2018). The signature of the cocoon can appear similar to that of the jet due to the cocoon’s interaction with the circumburst medium producing synchrotron emission (Nakar & Piran 2017). However, the emission of this component is more likely observable in systems where the jet is viewed off-axis, or with a favourable combination of energetics and geometry, when the core jet cannot dominate the observed emission at all times.

The afterglow is usually visible between the radio and X-ray wavebands, even reaching GeV energies in the most luminous events (Ackermann et al. 2013). However in the past 3 yr, detections of the afterglow have been made at very high energies (VHE, >100 GeV; MAGIC Collaboration 2019; H. E. S. S. Collaboration 2021). There are now a handful of VHE detections made seconds to hours after the lGRB prompt emission. Their light curves are very similar to that observed in the X-ray band, implying a connection to the afterglow rather than the prompt emission (H. E. S. S. Collaboration 2019; Blanch et al. 2020a). Different production mechanisms have been invoked to explain the VHE emission: GRB 190114C has been best fit with a synchrotron self-Compton component (MAGIC Collaboration 2019), whereas the GRB 190829A data set has been best described using a single synchrotron emission component spanning from radio to VHE (H. E. S. S. Collaboration 2021).

VHE detections associated with lGRBs are limited in redshift due to pair production between the VHE photons and the extragalactic background light: optical and infrared photons produced by star formation processes. As a result, VHE photons from sources above redshift 1.5 are highly attenuated, and are not expected to be detectable. This is supported by the VHE GRBs detected thus far: all at redshifts of about one or below (Vreeswijk et al. 2018; Castro-Tirado et al. 2019; Valeev et al. 2019; Izzo et al. 2020a). We note that all of the VHE GRBs to date have been associated with strong radio detections (MAGIC Collaboration 2019; Marcote et al. 2020; Rhodes et al. 2020a, b).

Once the photons have left the GRB jet, they propagate through the interstellar medium (ISM) of the host galaxy, the intergalactic medium, and the ISM of the Milky Way. These media along the line of sight can dramatically affect the observed afterglow emission. Some lGRBs have appeared to be optically faint, the so called dark GRBs. Chandra & Frail ( 2012) showed that about 25 per cent of Swift GRBs do not have detected optical counterparts. There are three possible explanations for dark GRBs (Jakobsson et al. 2004; Resmi et al. 2005): (1) they occur at high redshift resulting in the Lyman break falling in the optical band, (2) dust along the line of sight absorbs the optical photons, and (3) additional emission components at X-ray energies, increasing the X-ray flux with respect to the optical emission.

For many dark GRBs, it is dust along the line of sight that causes the optical darkness. Significant dust is expected in the regions of lGRBs as they occur in areas of high star formation, near to the birth-sites of their progenitor, since their lifetime is short. In some cases the extinction is greater than 10 magnitudes (V-band; e.g. Zauderer et al. 2013). When compared to the measured neutral hydrogen column densities inferred from the X-ray spectra, such high optical extinction deviates strongly from the linear AV–NH relationship measured within the Milky Way NH ≈ 2 × 1021AV (Predehl & Schmitt 1995; Güver & Özel 2009). It is likely that this is a result of a combination of the AV–NH relation varying from galaxy-to-galaxy, in combination with a non-uniform distribution of gas within each galaxy. This is supported by observations of dark GRBs’ host galaxies which appear to have ‘normal’ colours, implying a lack of increased dust across the galaxy as a whole but rather localized to regions of increased star formation (Perley et al. 2009).

Material along the line of sight can also affect radio emission in the form of interstellar scintillation (ISS; Goodman 1997; Walker 1998). Turbulence in the Milky Way’s ISM causes flux density fluctuations at radio frequencies up to the order of unity. ISS can be divided into weak and strong scintillation. Weak scintillation occurs above some characteristic transition frequency. Strong scintillation occurs at frequencies below the transition frequency and can be further divided into diffractive (DISS) and refractive scintillation (RISS). DISS is a narrow-band effect resulting from multipath propagation whereas RISS, a broad-band effect, occurs due to the focusing and defocusing of rays as they propagate through the ISM. If observed, DISS will dominate at early times when the GRB jet is more compact, but as the size of the jet on the sky grows, the effects of DISS fade away leaving RISS (Frail et al. 1997). The effects of RISS will also quench at some time when the jet has expanded beyond a certain angular size on the sky. Depending on the observing frequency and angular size of the source, the variability can be up to 100 per cent. The angular size dependence of DISS and RISS can be used to place on constraints on the size of the jet at different epochs (Frail et al. 1997, 2000; Chandra et al. 2008; Alexander et al. 2019).

The prompt emission from GRB 201216C was detected on 2020 December 16 at 23:07:31 ut by the Neil Gehrels Swift Observatory (here after Swift) Burst Alert Telescope (BAT; Beardmore et al. 2020). Three optical observatories also reported detections of a counterpart from early-time observations. A team searching for the afterglow with the Very Large Telescope (VLT) detected a source within the BAT error region at 21.81 ± 0.05 magnitudes (r’-band) 2.19 h after the burst (Izzo, Malesani & Kann 2020b). The VLT also measured a very steep optical spectral index (ν−4.1 ± 0.2) and placed GRB 201216C at redshift z = 1.1 (Vielfaure et al. 2020). Jelinek et al. ( 2020) and Shrestha et al. ( 2020) confirmed the optical source as the afterglow. A number of other observatories reported deep upper limits (Belkin et al. 2020; Oates, Beardmore & Swift/UVOT Team 2020; Gokuldass et al. 2021). There was no report of a detection of a supernova component. The reported optical detections and upper limits are shown in Fig.  1 as the squares and downwards facing grey triangles, respectively.

X-ray, optical, and radio observations from GRB 201216C. The flux densities for the radio data points are given in Table  1. The optical flux densities and upper limits are from the Gamma-ray burst Coordinates Network Circulars (Belkin et al. 2020; Izzo et al. 2020b; Shrestha et al. 2020; Gokuldass et al. 2021). The Swift-XRT light curve for GRB 201216C has been rebinned into 5 min bins. The inset shows a clearer view of the radio data set.

The Swift-X-ray Telescope (XRT) started observing ∼50 min post burst. The XRT unabsorbed fluxes were very high with respect to the optical counterpart. When combined with the steep optical spectral index, GRB 201216C was classified as a dark GRB (Vielfaure et al. 2020). It is unlikely that such a steep optical spectral index is a result of galactic extinction as the reddening in the direction of the burst is E(B–V) = 0.05 (Oates et al. 2020). We discuss this further in Section  3.4.

The Major Atmospheric Gamma Imaging Cherenkov (MAGIC) telescope, observing between 50 GeV and 50 TeV, reported the detection of a significant VHE counterpart less than a minute after the initial burst detection (Blanch et al. 2020b), making GRB 201216C is the highest redshift VHE GRB to date. Upon the notification of a VHE detection from the MAGIC Collaboration, we began a multifrequency radio campaign with a series of successful Director’s discretionary time observations (DDTs) with e-MERLIN, the Karl G. Jansky Very Large Array (VLA) and MeerKAT. We also applied for late-time Target of Opportunity (ToO) observations with Swift-XRT. All of these observations are further discussed in Section  2. In Section  3, we present the temporal and spectral evolution of the afterglow of GRB 201216C. In Section  4, we discuss the possible interpretations of the data using different jet models and discuss the ISM in both the Milky Way and the GRB’s host galaxy.

Here, we present the multifrequency observations obtained and the data reduction process. A list of observing dates, peak flux density measurements, and uncertainties are given in Table  1. Spectral indices are only calculated in epochs where we have high enough signal-to-noise ratios.

A table of the radio observations made with e-MERLIN, the VLA, and MeerKAT. The columns are the following: T–T0, the time between the burst detection and the centre of the observation, in days; ΔT, the duration of the observation, in hours; Δν, the observing frequency range; Sν, the peak flux density (or 3σ upper limit); β, the in-band spectral index. The uncertainties on each flux density measurement are a combination of the fitting error and a calibration error (5 per cent for e-MERLIN and VLA, 10 per cent for MeerKAT) added in quadrature. We only give values for β for epochs when the source is bright enough to be detected in at least one half of the band.

A table of the radio observations made with e-MERLIN, the VLA, and MeerKAT. The columns are the following: T–T0, the time between the burst detection and the centre of the observation, in days; ΔT, the duration of the observation, in hours; Δν, the observing frequency range; Sν, the peak flux density (or 3σ upper limit); β, the in-band spectral index. The uncertainties on each flux density measurement are a combination of the fitting error and a calibration error (5 per cent for e-MERLIN and VLA, 10 per cent for MeerKAT) added in quadrature. We only give values for β for epochs when the source is bright enough to be detected in at least one half of the band.

GRB 201216C was observed by e-MERLIN three times through successful DDT proposals at 5, 12, and 29 d post burst (PI: Rhodes, project codes: DD10010 and DD11001). We observed at 5 GHz with a bandwidth of 512 MHz. Each epoch consisted of 60 min on the flux calibrator (3C 286) and 90 min on the bandpass calibrator (OQ208) followed by 8 h of interleaved target and phase calibrators cycles: 6 min on the target and 2 min on the phase calibrator (J0056+1625).

We used the e-MERLIN pipeline to reduce the observations 1 (Moldon 2021). The pipeline performs flagging, delay, and bandpass calibration, and calculates phase- and frequency-dependent amplitude gain corrections which are all applied to the target field, along with flux density scaling from the flux calibrator. The calibrated measurement set was imaged in casa (Version 5.3.0) using the tclean task (McMullin et al. 2007). The uncertainties associated with the flux density measurements combine the statistical uncertainty and a 5 per cent calibration error.

Six VLA observations were obtained through a DDT proposal (PI: Rhodes, project ID: 20B-456). We spread the observations out between 12 and 53 d post burst. The observations were made at 10 GHz with a bandwidth of 4 GHz. For each epoch, we observed the target field for 10 min, book-ended with the phase calibrator (J0121+1149) and the primary calibrator (3C 147). We reduced the observations using the VLA pipeline in casa (Version 5.3.0; Kent et al. 2018). The pipeline performs flagging, creates a model of the flux calibration, and performs initial calibration including antenna position corrections. Delay, bandpass, and gain corrections are derived and applied to the data after which further flagging is performed. Imaging was also performed in casa . The uncertainties on the flux densities were calculated by combining the statistical error and 5 per cent calibration uncertainty added in quadrature.

We obtained four DDT observations with MeerKAT (PI: Rhodes, DDT-20210107-LR-01) at 22, 29, 40, and 54 d post burst. Each observation lasted 140 min, made up of a 5 min scan of a primary calibrator (J0408−6545) preceded by a series of 20 min scans of the target interleaved with 2 min scans of the secondary calibrator (J1808+0134). The observations were made at a central frequency of 1.28 GHz with a bandwidth of 856 MHz, split into 4096 channels.

The MeerKAT data were reduced using oxkat , a set of python scripts used for semi-automatic processing (Heywood 2020). First, the calibrator fields were flagged for RFI as well as the first and last 100 spectral channels. A spectral model from the primary calibrator was applied to the secondary. Delay, bandpass, and complex gain calibration was performed on the primary and secondary calibrators and applied to the target field. Finally the target field was flagged using tricolour . 2 The data were imaged with wsclean using a Briggs weighting with robust parameter of −0.7 (Offringa et al. 2014). We derived a model from the image and used it to reimage after a round of phase-only self-calibration. The flux uncertainties include statistical uncertainties and a 10 per cent calibration error.

The Swift X-ray telescope (XRT) observed the field of GRB 201216C from 3000 s until 22 d after the initial burst (Evans & Swift-XRT Team 2021). This included two ToO observations that we obtained between days 20 and 27 post burst. Each epoch was automatically fitted with a power-law spectrum. The light curve and spectra are made public on the Swift Burst Analyser (Evans et al. 2007, 2009, 2010). The X-ray flux densities used in our analysis are calculated at 5 keV to avoid systematic under or overestimations in calculating the flux density at the edge of the observing band (i.e. at 0.3 or 10 keV).

In the following sections, we use the convention Fν ∝ tανβ, where t is the time post burst, ν is the observing frequency, and α and β are the exponents. Any subscripts are used to indicate the relevant part of the spectrum or frequency band.

Fig.  1 shows the radio light curves from our observing campaign. The flux densities and upper limits are given in Table  1. We detected radio emission at 5 GHz in two of the three observations with e-MERLIN, during epoch one and three (the green crosses in Fig.  1).

Our 10 GHz light curve (blue stars in Fig.  1) from the VLA covers the largest time range, from 12 to 54 d post burst. The 10 GHz behaviour is best described as a shallow power-law decay (⁠|$\alpha _{\rm {10\,\, GHz}} = -0.5\pm 0.1$|⁠ ) from ∼120 |$\mu$| Jy at 12 d to ∼50 |$\mu$| Jy at 54 d. On top of the decaying flux, we detect interobservation variability, which is possibly due to ISS (see Section  3.2). The epoch-to-epoch variability could cause the observed decay rate to deviate significantly from the intrinsic evolution.

The 1.3 GHz MeerKAT light curve (the purple stars in Fig.  1) starts with a very steep rise (⁠|$\alpha _{\rm {1.3\,\, GHz}}\gtrsim 5$|⁠ ) from a 3σ upper limit of 29 |$\mu$| Jy to a detection of 95 |$\mu$| Jy in 6 d. The next three data points show a fairly flat light curve (⁠|$\alpha _{\rm {1.3\,\, GHz}} = 0.1+/-0.2$|⁠ ) to 130 |$\mu$| Jy in the final epoch. The sharpest rise possible for the standard forward shock model is t1.75, which is still far shallower than the observed rise, comes from optically thick synchrotron from a forward shock propagating through a stellar wind environment in the regime where νm < νobs < νSA (Granot & Sari 2002).

Due to the jet’s compactness, radio observations of GRB afterglows are susceptible to scintillation, which can cause significant spectral and temporal variability, especially at early times. We note that scintillation time-scales are also frequency-dependent and that we are sampling variability on different time-scales with the different interferometers due to differing observation lengths. ISS may be the cause of the interobservation variability seen with the VLA. We also search for intraobservation variability within the VLA data, as shown in Fig.  2. We observe some low-level variability, ∼10–20 per cent (using equation 10 from Vaughan et al. 2003), in VLA epochs one and three. VLA epochs two and four show no such variability (see Fig.  2). The flux densities of the last two observations are too low to search for variability.

Intraobservation light curves to show short-term variability in our VLA 10  GHz data set. The first four 10 min epochs are sub-divided into four −2.5 min sub-integrations. We only use the first four observations as these are the brightest four where the source is reliably detected on short time-scales.

We split the two e-MERLIN detections, which are 6 h long each, into four-90 min segments. The left-hand panel of Fig.  3 shows that for the first epoch, radio emission was only detected in two of the four segments. In the final epoch, we only detect radio emission for 90 min out of 6 h.

Short-term variability observed in the first and third e-MERLIN epochs. Each 6-h observation was split into four sub-integrations. In the first observation made 6 d post burst, we detected radio emission in two of the four sub-integrations. In the third observation, made at 28 post burst, we only detected the source in one of the four sub-integrations.

Our MeerKAT data set shows no evidence of intra-observation variability on a time-scale of tens of minutes. The implications of the observed variability are discussed further in Section  3.2.

An optical counterpart to GRB 201216C was detected with the Liverpool telescope and the VLT (red and orange squares in Fig.  1, respectively; Izzo et al. 2020b; Shrestha et al. 2020). Within a couple of hours, the source had faded below detection limits. The r’-band light curve follows a αr = −0.83 ± 0.01 decay from a few minutes post burst (Belkin et al. 2020; Gokuldass et al. 2021). Such a decay rate is consistent with optically thin synchrotron radiation in an ISM environment which gives p = 2.11 ± 0.01 (from α = 3(1 − p)/4; Granot & Sari 2002). An optically thin forward shock in an wind environment would decay more rapidly (for p = 2, α = −1.25), however, with only a handful of detections it is impossible to determine whether the measured decay rate could be due to a combination of a frequency break passing through the optical observing band and a stellar-wind-like environment. Such a combination would result in the optical light curve appearing shallower.

The XRT light curve (black circles in Fig.  1) consists of a significant number of detections early on (t < 1 d) and one further detection around 25 d post burst. We have rebinned the photon-counting mode light curve to reduce any bias towards the earlier detections when fitting a power-law decay to the light curve. Using either 5 or 10 min bins, we obtain a decay rate of αX-ray = −1.9 ± 0.3. The final detection is slightly above the predicted flux density given the above decay rate as shown in Fig.  5.

Various spectra at three separate epochs.

Top panel: X-ray, optical, 10 and 1.3 GHz data for GRB 201216C. Overlaid is a simple single forward shock model. For the optical light curve, we require an optical extinction between 4.3 and 4.6 magnitudes, this is far outside the range of inferred extinction and so we model the light curve with an extinction of 5.3  magnitudes: the lowest value in our inferred range. Bottom panel: the normalized residuals, the ratio of ‘Δ’ (the difference between each observed flux density and the model at that time) to σ (the uncertainty on each measured flux density).

There are many scenarios in which the X-ray light curve agrees with theoretical predictions. (1) Synchrotron radiation from a forward shock above νC in either a homogeneous or wind environment, α = (2−3p)/4, giving p = 3.2 ± 0.4; (2) optically thin forward shock synchrotron emission below νC in a stellar wind environment, α = (1−3p)/4, giving p = 2.9 ± 0.4 (Granot & Sari 2002); (3) optically thin forward shock emission below νC in a homogeneous environment, α = 3(1 − p)/4: p = 3.5 ± 0.4, and (4) the early X-ray light curve is also consistent with a jet break with no significant lateral spreading, with only edge effects considered. In a homogeneous environment, emission above the cooling break should decay as α = −(1 + 3p)/4, where p = 2.2 ± 0.2, and below the cooling break α = −3p/4, where p = 2.5 ± 0.2. In a stellar wind environment, emission below the cooling break should decay as α = −(1 + 3p)/4, where p = 2.2 ± 0.2, and above the cooling break α = −3p/4, where p = 2.5 ± 0.2 (Gao et al. 2013). We can only break the degeneracy between the different potential scenarios by considering the X-ray spectral index measurements and light curves in other wave bands to form a broad-band model.

As well as the long-term evolution, our observations also show evidence of short time-scale (inter- and intra-epoch) variability as a result of ISS.

DISS causes narrow-band fluctuations of the order of unity on a range of time-scales and therefore, can affect radio observations dramatically. Our observations from e-MERLIN show evidence of short time-scale variability, a feature that is inconsistent with the smooth variations expected from GRB afterglows. Here, we explain that the variability observed is a result of small-scale inhomogeneities in the local ISM which causes multipath propagation (DISS) of the radio waves from GRB 201216C (Goodman 1997). For simplicity, the region of the ISM causing the scattering is collapsed into a screen as some distance along the line of sight (Walker 1998).

We placed a lower limit on the intra-observation temporal variability of 30 per cent for our e-MERLIN data set (Vaughan et al. 2003). If the observed variability is due to DISS, we would expect to see narrow-band flux modulations up to one on the time-scale of an hour (Goodman 1997; Walker 1998). Due to signal-to-noise limitations, we cannot search for shorter time-scale variability. We also searched for variability in the spectral domain by dividing the two detection epochs into four sub-bands (centred at 4.8, 4.9, 5.1, and 5.2 GHz), another sign of DISS since DISS is a narrow-band phenomenon. In the first e-MERLIN epoch, radio emission is only detected in the sub-band centred at 4.8 GHz, at 300 ± 30 |$\mu$| Jy. In the bands centred at 4.9, 5.1, and 5.2 GHz, we obtained 3σ upper limits of 156, 144, and 195 |$\mu$| Jy, respectively. A high level significance detection in only a single, narrow frequency band implies that DISS is most likely the origin of the variability at 5 GHz. In the final e-MERLIN epoch, the low flux density of the source means we are unable to detect emission in any of the sub-bands.

Under the assumption that the e-MERLIN variability is caused by DISS, we are able to place constraints on the location of the scattering screen between the Earth and the position of the GRB. If we assume that the screen is located within the Milky Way, by integrating the free electron distribution along the line of sight using the NE2001 model (Cordes & Lazio 2002), we calculate a scattering measure SM−3.5 = 0.69, where SM−3.5 = SM/(⁠|$10^{-3.5}\, \textrm {kpc}~\textrm {m}^{-20/3}$|⁠ ), defined as the characteristic angle that incoming radio waves are scattered by whilst propagating through the ISM and from it infer a transition frequency of 8.8 GHz. Observations below 8.8 GHz are in the regime where it is possible to observe strong scattering, consistent with our conclusion that the variability at 5 GHz is produced by DISS. This scattering measure and transition frequency correspond to a scattering screen at a distance (dscr) of ∼0.9 kpc. We obtain the same value for dscr using the method presented in Goodman ( 1997).

We know that the radio emission observed is still affected by DISS 29 d post burst, based on the short time-scale variability observed in our last e-MERLIN observation (Fig.  3). Therefore, at 29 d post burst the angular size of the jet associated with GRB 201216C must be less than 1μas. At a redshift of 1.1, the distance to GRB 201216C, 1 |$\mu$| as is ∼1 × 1017cm. This size upper limit is consistent with size measurements of other lGRBs made at around 30 d (see fig.  7 of Alexander et al. 2019; Taylor et al. 2004). It is likely that if subsequent observations at 5 GHz were made, they would most likely not have been affected by DISS.

A plot showing the evolution of the break frequencies for both the narrow and wide components from Peng, Königl & Granot ( 2005)’s two component afterglow model, adapted for a stellar-wind environment applied to our data set. The purple dashed, dotted, and solid region denotes the movement of νSA, νm, and νC from the narrow jet, respectively. The dotted-dashed line shows the evolution of νm from the wide component viewed off-axis. The purple shaded regions denote the uncertainties in the location of each frequency break. Overlaid in green, blue, red, and grey are the MeerKAT, VLA, optical, and XRT observing bands, respectively.

A schematic of the geometry of VHE GRB 201216C. The narrow core component is the origin of the VHE, optical, X-ray, and early (< 25 d) radio data. The wider cocoon component produces the late-time radio emission.

According to the thin screen scattering model, at 1.3 GHz (the MeerKAT observing band), we would expect to see variability on time-scales of 10 min with a modulation index of one. We see no intra-observation or narrow-band variability in our MeerKAT data. This indicates that by 29 d post burst, the angular size of the jet has grown larger than the 0.3 |$\mu$| as (0.3 × 1017cm).

where ν0 and ν are the transition and observing frequency, respectively.

The observations at 5 GHz are dominated by the effects of DISS, and due to the sparse cadence, days between each epoch, we are unable to observe the effects of RISS. In the MeerKAT band, the increase in flux density between days 23 and 29 is greater than a factor of three, far higher than the predicted RISS flux modulations of 30 per cent. The observations in which the source detected shows a smooth increase in flux density across the three epochs in which we detect the source. The spacing between each epoch is too large to infer whether the increase in flux density is due to RISS. Therefore, we cannot confidently attribute the MeerKAT flux variations to RISS.

Weak scintillation often affects the data at a level similar to that of the calibration uncertainties (⁠|${\sim}5{-}10{{\ \rm per\ cent}}$|⁠ ); although it can be significantly higher for observing frequencies close to the transition frequency). Our observations at 10 GHz, above the transition frequency, show clear inter-observation variability (the blue stars in Fig.  1), as well as at the ∼10 per cent level on the time-scale of minutes in the two of the first four epochs, see Fig.  2. The flux density of the radio counterpart in the last two VLA epochs are too low to search for intra-observation variability. As a result, we cannot tell if the effects of weak scintillation have faded as the jet grows on the sky.

With a transition frequency of 8.8 GHz, the majority of the VLA observing band falls within the weak scattering regime. However, the VLA’s wide bandwidth means we probably still observe some effects of DISS and RISS. DISS, RISS, and weak ISS are expected to cause variation on time-scales of two to three hours across the VLA band with a modulation index as high as 1 (Walker 1998; Granot & Van der Horst 2014). Such high variability levels are to be expected because our observations are so close to the transition frequency, although the amplitude of the flux modulation is expected to drop rapidly towards high frequencies. Therefore, it is most likely that the variability observed in the VLA band is a combination of DISS and weak scintillation.

The right-most column of Table  1 shows all the in-band radio spectral index measurements calculated using the individual observing bands for epochs where the source was bright enough. Of the three e-MERLIN observations, only the first epoch was bright enough to obtain an in-band spectral index. The scintillation dramatically affects the e-MERLIN spectra as it does the intra-epoch light curves. At 5 d post burst, the only time where the source is bright enough to split the band in two, we measure a 4.8−5.3 GHz spectral index of 7 ± 4. The large uncertainties mean that such a steep result is still compatible with the steepest branch of the synchrotron spectrum in the GRB afterglow scenario, caused by synchrotron self-absorption.

The radio emission at 10 GHz (VLA) is only bright enough in the first three epochs to split the 4 GHz bandwidth into two-2 GHz subbands. In each of these three epochs, the 8–12 GHz spectral index is consistent with being spectrally steep or fairly flat (⁠|$\beta _{\rm {10\,\, GHz}}\ge$| 0). Over the course of the three observations, the VLA in-band spectral index slowly flattens (see Table  1). The wide VLA observing band smears out any narrow-band effects of DISS. In the context of the fireball model, the observations made at 12 and 14 d post burst are too steep to be in the regime where νSA < 10 GHz < νm (β = 1/3). Instead, they are more consistent with 10 GHz < νSA, νm (β = 2), or νm < 10 GHz < νSA (β = 2.5; Granot & Sari 2002).

We also calculated the 1.0–1.7 GHz spectral index for the three MeerKAT observations where we detect radio emission (values are also given in Table  1). At 28 d post burst, we were only able to detect radio emission in the upper half of the band, so we only have a lower limit on the spectral index, here the spectral index is approximately flat. The final two epochs show a spectral index of |$\beta _{\rm {1.3\,\, GHz}}$| <0, consistent with optically thin synchrotron: νm, νSA < 1.3 GHz. The three MeerKAT detections are made after the epochs where the 8–12 GHz spectral indices are calculated, meaning that over the course of the radio campaign, the emission evolves from being optically thick to optically thin.

The Swift-XRT spectra in the range of 0.3–10  keV are each fitted with a power law parametrized by the photon index: Γ, where βX = 1 − Γ. There are no significant variations in Γ over the observing period implying that no break frequency passes through the XRT band. The data taken in photon counting mode result in Γ = 2.0 ± 0.1 (βX = −1.0 ± 0.1). As with the X-ray light curve, βX, emission both below (p = 3.0 ± 0.2) and above (p = 2.0 ± 0.2) νC. The spectral index does not change in the event of a jet break, unlike the light curves.

There are three epochs where we construct broad-band spectra. Fig.  4(a) shows the optical and X-ray data detections (black circle and square), approximately two hours post burst. The blue shaded region denotes the range of possible predicted flux densities under the assumption that νC falls between the optical and X-ray bands at the time of the observations. The range is calculated assuming that βO − X is between βX + 0.5 and βX i.e. νC is at 0.3 keV and in the r’-band, respectively. It is clear from the blue shaded region that the optical flux density is significantly lower than expected from the synthesized synchrotron spectrum. We calculate that βO − X = −0.13 ± 0.02. Given that βX = −1.0 ± 0.2 at this time, if νC falls between the optical and lower end of the X-ray band, we can infer that βO − X should be − 0.5 ± 0.2 (βO − X = βX + 0.5), much steeper than our measured βO − X (Van der Horst et al. 2009). We discuss this classification further in Section  3.4.

Using the detections and upper limits between 20 and 24 d post burst (Fig.  4b), we can construct a broad-band spectrum from 1.3 GHz to 10 keV. Around 20 d, |$\beta _{\rm {10\,\, GHz-X}} = -0.6$|⁠ , which differs from βX = −1.0 ± 0.1 at a 4σ level indicating that it is likely that νC is below the XRT band. We use the value of p (2.0 ± 0.2) derived from the X-ray spectrum and constrain the location of the cooling break at 20 d to be between 8|$\times 10^{15}\, {\rm Hz} \lt \nu _{\textrm {C}} \lt 8 \times 10_{17}\, {\rm Hz}$|⁠ . We also infer the position of the peak of the spectrum and the corresponding flux density: 13 ± 9 GHz and 130 ± 30 |$\mu$| Jy, respectively. The resulting broad-band spectrum from around 20 d post burst can be well described with a series of three power laws, as shown in Fig.  4(b): (1) a low frequency steep spectral component, (2) an optically thin branch where β = −0.5 between νpeak < ν < νC, and (3) the final branch above the cooling break where β = −1.0 ± 0.1. The shaded regions denote the uncertainties (for the optically thin and cooling branches) or variations in possible spectral indices (low frequency optically thick branch).

At 54 d post burst, the broad-band radio spectrum is best described by a single power-law component (Fig.  4c) with a spectral index of βrad = −0.50 ± 0.02, most likely from the optically thin branch of the spectrum below νC where p = 2.00 ± 0.04. This means that by 54 d post burst the peak of the synchrotron spectrum is below 0.9 GHz.

Early time optical observations either placed deep upper limits on any optical emission or obtained very faint detections of the afterglow with respect to the X-ray fluxes, see e.g. Fig.  4(a). Furthermore, the optical spectrum observed by the VLT was very steep, which provided concrete evidence that GRB 201216C is a dark GRB (Vielfaure et al. 2020).

From the optical and X-ray light curves, as well as the broad-band and optical spectra, we can infer that νC is between the optical and X-ray observing bands from 0.05 to ∼1 d post burst. We do not consider the final X-ray data point here because it is a low significance detection. We use this information to place limits on the r’-band flux density if GRB 201216C was not heavily affected by extinction as described in Section  3.3.3. By comparing the inferred and measured optical flux densities, we estimate the extinction to be between 5.3 and 8.6 magnitudes (r’-band); values far in excess of the galactic extinction contribution: AR = 0.12 mag (E(B–V) = 0.05; Schlegel, Finkbeiner & Davis 1998; Fitzpatrick 1999).

Using the empirical relations between extinction and neutral hydrogen column density: NH ≈ 2|$\times 10^{21}\, {\rm cm}^{-2}{\rm A}_{\rm V}$|⁠ , (Predehl & Schmitt 1995; Güver & Özel 2009), we can determine whether or not the line-of-sight hydrogen column density is consistent with the attenuated optical flux densities. For |$N_{\rm H} = 5.07\times 10^{21}\, {\rm cm}^{-2}$|⁠ , from X-ray spectra, we estimate AV ≈ 3 mag (AR ≈ 2 mag). Therefore, an additional source of optical extinction is required. From the observed extinction range, we would expect |$N_{\rm H} = 1-3 \times 10^{22}\, {\rm cm}^{-2}$|⁠ , obtained from the X-ray spectra, which is at least a factor of two higher than the measured NH value.

Given that GRBs occur in regions of high star formation, increased dust in the vicinity of the GRB site is expected, so optically dark GRBs should not be uncommon (Fruchter et al. 2006). Studies of the host galaxies of dark GRBs have shown that the dust distribution is non-uniform, further agreeing with the previous statement that dark GRBs occur in highly obscured regions (Perley et al. 2009). Giant molecular clouds could also be a contributing factor to increased amounts of dust in the vicinity of lGRBs, but would also result in higher measured NH values (Solomon et al. 1987).

The fact that we do not observe such a high NH may also be a result of the AV–NH correlation varying from galaxy-to-galaxy, especially at high redshift (z > 1) where the star formation rate is much higher than in local galaxies. The above calculation assumes a universal AV–NH relation, and so is not necessarily correct for GRB 201216C’s host galaxy. GRB 110709B’s optical darkness was similarly underpredicted by the measured hydrogen column density (Penacchioni et al. 2013; Zauderer et al. 2013). On the other hand, the optical extinction for many GRBs is overestimated, again implying a clear deviation from the Galactic and Magellanic Cloud relations (Perley et al. 2009; Krühler et al. 2011).

Of the four other VHE GRBs detected so far, GRBs 190829A and 190114C have also shown increased optical extinction (Campana et al. 2021; Zhang et al. 2021). Zhang et al. ( 2021) measured an absorption E(B–V) = 0.757 for GRB 190829A and Campana et al. ( 2021) obtained E(B–V) = 0.83 for GRB 190114C. Such high extinctions in the most well-studied VHE GRBs could suggest a potential connection between high density/dusty environments and the presence of VHE emission. Dusty environments could result in strong infrared radiation fields following the reprocessing of optical emission. The infrared radiation could be upscattered to VHE energies in the presence of electrons with sufficiently high Lorentz factors (γe ∼ 106). Further exploration of this idea is outside the scope of this paper and will be considered in future work.

The most simple description of the multifrequency data would be a single forward shock which is produced as the jet decelerates in the circumburst medium. Using our constraints on the positions of the break frequencies and the peak flux from our broad-band spectral considerations (Section  4), combined with well-established analytical afterglow models (e.g. Granot & Sari 2002), we can determine if a single forward shock describes our data well.

In Section  3.3.3, we constrained the location of νC at 20 d post burst to be between 8 × 1015 and 8 × 1017 Hz. Given that we observe no statistically significant breaks in the light curve before day 20, we can assume that until at least day 20, the X-ray emission must originate from above νC. For an X-ray decay of αX = −1.9 ± 0.2, using the binning shown in Fig.  1, we obtain p = 3.2 ± 0.3, which is steeper than p derived from the X-ray spectra (2.0 ± 0.4) at nearly 2σ. These p values are independent of the circumburst medium. Despite the light-curve slope being same, the movement of νC changes depending on the environment. In a stellar wind environment |$\nu _{\textrm {C}}\propto t^{\frac{1}{2}}$| and in a homogeneous medium |$\nu _{\textrm {C}} \propto t^{-\frac{1}{2}}$|⁠ . Given the inference that νC is only just below the XRT band at day 20; if the jet was propagating through a homogeneous environment, we would expect to observe a break in the X-ray light curve due to νC at some time before day 20. Therefore, we conclude that the XRT emission is most likely a result of a stellar-wind environment where νC moves from lower to higher frequencies with time, i.e. towards the XRT band.

The model X-ray light curve is shown as the black line in Fig.  5. From either the X-ray light curve or spectra, it is impossible to determine the location of νC with respect to the 0.3–10 keV band. However, when combining the two with the broad-band spectrum at 20 d, we have good evidence that the jet is propagating through a stellar wind environment. Fig.  6 shows the movement of νC with time according to our stellar wind model (solid purple line). The shaded purple region around the line represents the uncertainty on the location of νC at a given time, derived from the broad-band spectrum constructed around day 20 (Fig.  4b).

The intrinsic brightness of the optical emission is heavily absorbed, but we can assume that the decay rate (t−0.83 ± 0.01) observed is unaffected by the material causing the absorption. For a wind environment, such a decay is too shallow to be produced by optically thin synchrotron in which the r’-band is below νC. The decay, αopt, should follow −2.0 < αopt < −1.3, using an optically thin decay in a stellar wind environment for 2 < p < 3. The shallower decay that we observe could be a result of νm passing through the optical observing band within a few hours of the burst (as previously mentioned in Section  3.1.2). We use this to constrain the position of νm at early times as the purple dotted line in Fig.  6. The red line in Fig.  5 shows the model light curve for a forward shock where νm causes the break around 0.01 d. In order to best fit to the optical data, we consider ∼5.3 magnitudes of optical extinction, the lowest value in the extinction range calculated in Section  3.4. Fig.  5 shows that our forward shock model does not fit the optical data well.

The similar flux densities, to within an order of magnitude, are observed by both XRT before 0.1 d post burst and at the beginning of the VLA observing campaign. We estimate that the peak of the spectrum should be in the optical band around 0.01 d. Fig.  5 shows that at 0.1 d, the most optical and X-ray flux densities are approximately the same. However, when we consider the extinction calculated in Section  3.4, it is clear that |${\rm F}_{\nu ,\, {\rm max}}$| should be significantly higher at this time. The VLA observations, which start at day 20 and show flux densities similar to that measured by XRT 0.04 d post burst, show evidence of the spectral peak passing through the 10 GHz band, meaning that |${\rm F}_{\nu ,\, {\rm max}}$| must have decayed from 0.01 to 20 d post burst to show similar X-ray and radio flux densities. If we assumed an ISM environment, |${\rm F}_{\nu ,\, {\rm max}}$| would be constant with time and would not describe the data at all. Instead, this can be attributed to either a jet break or a stellar-wind environment. A jet break can cause |${\rm F}_{\nu ,\, {\rm max}}$| to decay rapidly as a result of the jet expanding laterally (Fν,max ∝ t−1; Sari, Piran & Halpern 1999). A stellar wind environment means that |${\rm F}_{\nu ,\, {\rm max}}$| will decay with time as |${\rm F}_{\nu ,\, {\rm max}} \propto t^{-0.5}$| when νm > νSA, and |${\rm F}_{\nu ,\, {\rm max}} \propto t^{-\frac{4p+1}{2(p+4)}}$| (t−(0.75−0.93) for 2<p<3) when νSA > νm (Granot & Sari 2002).

The VLA in-band spectral indices of the first three observations are steep (see Table  1). Such values further imply that the peak of the synchrotron spectrum is above 12 GHz until at least 20 d post burst (see the upper panel of Fig.  4). As explained in Section  3.1.1, if the peak of the spectrum was caused by νm then we would expect to measure a spectral index of β = 1/3. The steepness of the spectral indices implies that either νSA and νm or just νSA must be above the VLA observing band. The third VLA epoch is more consistent with a flat spectrum, i.e. β = 1/3 or the peak of the spectrum being at the observing frequency. Combined with the decaying flux density of the final three epochs, we can conclude that the peak of the synchrotron spectrum has passed through the 8–12 GHz band between 20 and 36 d post burst. The most simple scenario for the decay is if only one spectral break is above the VLA band, which moves towards lower frequencies with time, causing the emission above the break to decay with time. In the regime where νm < νSA, |$\nu _{\textrm {SA}} \propto t^{-\frac{3(p+2)}{2(p+4)}}$| (t−(1.0−1.1) for 2<p<3; Granot & Sari 2002). Therefore, we conclude that the peak of the spectrum is caused by νSA.

At 10 GHz, the forward shock model would show a broken power law, with a rise following |$t^{\frac{5}{4}}$| to a peak around day 20 followed by a decay of |$\alpha _{\rm {10\,\, GHz}} = \frac{(1-3p)}{4}$| (t−(1.3−2.0) for 2 < p < 3) as νSA moves through the observing band (Granot & Sari 2002). The theoretical 10 GHz light curve is shown (where p = 2.0) with the blue line in Fig.  5. The optically thick to thin transition provides a reasonable fit to the 10 GHz light curve. We note that the variability due to weak scintillation increases the residuals as shown in the lower panel of Fig.  5.

Fig.  4(c) shows that by 54 d post burst, the 1.3 and 10 GHz light curves are on the same, optically thin branch, of the synchrotron spectrum, with p = 2.00 ± 0.05. The MeerKAT in-band spectral index is also similar at day 41 post burst indicating that by 41 d post burst, 1.3 GHz is above νSA and νm. In order for νSA to be below the MeerKAT band at 41 d post burst, we required νSA ∝ t−3.5, which is significantly faster than the theoretical movement where for p ≈ 2, |$\nu _{\textrm {SA}} \propto t^{-\frac{3(p+2)}{2(p+4)}} = t^{-0.8}$| (Granot & Sari 2002). The measured νSA movement is unphysical when compared to analytical models.

The unphysical movement for νSA complicates the 1.3 GHz model light curve. Using νSA ∝ t−0.8, the model light curve would consist of a single power-law component with |$\alpha _{\rm {1.3\,\, GHz}} = \frac{5}{4}$|⁠ . The light curve would turn over at 200 d as a result of νSA entering the observing band if νSA moved as given in Granot & Sari ( 2002). During the rise, we would expect a spectral index |$\beta _{1.3\, {\rm GHz}} = \frac{5}{2}$|⁠ , not ∼−1 as measured. Furthermore, the peak flux density, according to our single shock model would continue to decay as |${\rm F}_{\nu ,\, {\rm max}} \propto t^{-\frac{4p+1}{2(p+4)}} = t^{-0.75}$| for p = 2, so by the time νSA (the peak of the spectrum) reaches the MeerKAT band, |${\rm F}_{\nu ,\, {\rm max}} \approx 20\,\mu$| Jy, which is a factor of 7 fainter than the observed 1.3 GHz flux density. When compared to the MeerKAT data points (purple circles and downwards facing triangle), the 1.3 GHz model light curve (the purple line) in Fig.  5 show a clear deviation away from the forward shock model. At all times, the model light curve falls far below the observed data. Even if we assume the unphysical movement of νSA, |${\rm F}_{\nu ,\, {\rm max}}$| still decays with time meaning that the modelled 1.3 GHz light curve is predicted to be much fainter than the observed emission.

The late-time change in the evolution of |${\rm F}_{\nu ,\, {\rm max}}$| and νSA could be a result of a change in the circumburst environment. A varying circumburst density distribution, i.e. a deviation from ρ ∝ r−2, where ρ and r are the density and radius for the burst site, could occur as a result of the progenitor star having fluctuating mass-loss rates towards the end of its life. In order to reproduce the observations, we require the circumburst environment to change from ρ ∝ r−2 to ρ ∝ r0: an homogeneous environment. Such a change seems unlikely to reflect the mass-loss history of the progenitor. Furthermore, this cannot explain the MeerKAT light curve or the discrepancy in p derived from the X-ray light curve and spectra.

In conclusion, Fig.  5 shows the results of a single forward shock component model overlaid on the X-ray, optical, 10 and 1.3 GHz light curves. We do not use the e-MERLIN 5 GHz data points in our model as they are heavily affected by DISS. The bottom panel of Fig.  5 shows the normalized residual values for our forward shock model with respect to the data. We are able to reproduce the X-ray and VLA light curves reasonably well. The interobservation variability increases the residuals of the VLA data with respect to the model. However, it is clear from the large residuals for optical and MeerKAT light curves that a single forward shock model is not a good fit. In a single shock scenario, we expect any variation in |${\rm F}_{\nu ,\, {\rm max}}$| to be dictated by the movement of νSA and νm. We acknowledge that our understanding of how |${\rm F}_{\nu ,\, {\rm max}}$| evolves with time early on is poorly constrained because of the faint optical emission. By inferring the range of optical flux densities from the X-ray data, we know that we require a steeper |${\rm F}_{\nu ,\, {\rm max}}$| decay than the afterglow models provide (Granot & Sari 2002). On the other hand, by the time the radio campaign begins, |${\rm F}_{\nu ,\, {\rm max}}$| appears to be constant in time. The rapid decay of |${\rm F}_{\nu ,\, {\rm max}}$| until day ∼20 followed by a transition to a constant |${\rm F}_{\nu ,\, {\rm max}}$| for the rest of the observing campaign is too complex to be attributed to a single jet component.

It is possible that the discrepancies the light curves could be explained with the addition of an extra shock component. Unfortunately, the light curves in any one observing band are too sparse to search for reverse shock emission. For example, optically thin reverse shock emission decays far steeper than the optical light curve. However, we note that two data points are not constraining enough to fully eliminate the possibility of reverse shock contribution. Unfortunately, because our first data point from e-MERLIN at 5 d post burst is dominated by RISS, we do not know the intrinsic flux density value at this time and so cannot determine if there is any reverse shock contribution. The full 5 GHz flux density in the absence of RISS (double/half the observed flux density assuming order of unity variability) is not constraining enough to confirm or reject the presence of reverse shock emission.

The MeerKAT and VLA bands are less affected by ISS and therefore are more appropriate to search for reverse shock emission. We apply a similar methodology to that in Rhodes et al. ( 2020a) in order to determine whether the reverse shock makes a significant contribution to our 1.3 GHz light curve despite not being detectable at 10 GHz. If the reverse shock has faded at 10 GHz by 12 d post burst, then we can determine that the 29 d, whilst the peak of the reverse shock might be in the MeerKAT observing band, Fν,max corresponding to the reverse shock will be significantly fainter than the observed MeerKAT flux densities at this time. The strongest evidence against reverse shock emission in the MeerKAT band is that the sharp rise between 22 and 29 d post burst, inconsistent with the reverse shock scenario.

An additional forward shock component could instead explain the large residuals between the data and the single shock model for the later observations at 1.3 GHz. In order to determine whether an additional forward shock could explain the discrepancies found in Section  4.1, we use the two-component jet model presented in Peng et al. ( 2005) adapted for a stellar wind environment (Chevalier & Li 1999) to better interpret our data, similar to what was applied to GRB 130427A by Van der Horst et al. ( 2014). Structured jets encompass a broad range of geometries and multiple jet components have been invoked in previous GRB afterglow data sets (e.g. Resmi et al. 2005; Racusin et al. 2008). The afterglow from gravitational wave event GW 170817 was inferred to have an ultrarelativistic core surrounded by lower velocity wings (Margutti et al. 2018). In GRB 080319B, the presence of two distinct outflow components was inferred, but the wider component was still more collimated than what we infer for our narrow jet (Racusin et al. 2008). The presence of a wider component is strongly supported by simulations (e.g. Morsony, Lazzati & Begelman 2007), although observations of such a component span a broad range of energetics and opening angles.

For our data on GRB 201216C, we consider a narrow, ultrarelativistic jet launched at a Lorentz factor, Γ0,n > 100, like with the single shock scenario, but with the addition of a wider outflow with Γw,0 < 10 (the subscripts n and w refer to narrow and wide, respectively). It is possible for the wider outflow to be non-relativistic, but we find that a relativistic outflow is more likely given the high luminosity and light-curve behaviour. We can rule out the possibility of supernovae emission as the origin of the late time radio detections: the radio luminosities are an order of magnitude higher than the next most luminous type Ib/c broad-line supernova (Bietenholz et al. 2021).

Each component is uniform within some opening angle θj,n, θj,w. The two components do not interact with each other. In the event of a jet break, we only consider edge effects, not lateral spreading. An on-axis observer would see that the narrow jet dominates at early times. The wider component becomes visible only in systems with favourable geometries and energetics. Fig.  7 shows the geometry of such a system, which undergoes a jet break early on (∼0.05 d) allowing the wider outflow to be visible later on.

In Section  4.1, we noted that the evolution of |${\rm F}_{\nu ,\, {\rm max}}$| is inconsistent with a single forward shock. We can explain a steep decay of Fν,max as well as the X-ray light curve if the narrow component undergoes a jet break before the first X-ray data points at around 0.05 d post burst. As with the single jet scenario, the X-ray light curve follows a single power-law decay. For a jet break scenario in a wind environment, we obtain p = 2.2 ± 0.2, which agrees with p from the X-ray spectra (p = 2.0 ± 0.4) at a 1σ level, slightly better than in the single shock scenario where p = 3.2 ± 0.4 (Section  3.1.3).

At 2 h post burst, we infer from Figs  4(a) and  6 that 3 ≲ |${\rm F}_{\nu ,\, {\rm max}}$| ≲ 15 mJy. By 20 d, |${\rm F}_{\nu ,\, {\rm max}}=$| 0.13 ± 0.03 mJy (Fig.  4b). We find that the best explanation of the steep decay of Fν,max is a result of a combination of the jet break and the wind environment. The combination creates a steeper decay of |${\rm F}_{\nu ,\, {\rm max}}$| compared to a stellar wind-only decay as used in the single shock scenario. The steepening of the light curves due to a jet break is a correction factor of 0.75 (Panaitescu, Mészáros & Rees 1998). Once we consider the optical extinction, we can match the inferred optical flux density from our model with observed optical flux densities. We use the spectral evolution to explain the optical detections here as in Section  4.1: νm passing through the band resulting in an apparent flattening of the optical light curve. The dotted red line in Fig.  8 shows a broken power law where the break is a result of νm passing through the band, in addition to the jet break at 0.05 d. The optical emission is best fit with an extinction range between of 7.4 and 8.0 magnitudes which is within the extinction range inferred from the X-ray observations (5.3–8.6 magnitudes). For the multiple component jet model, we require a larger extinction value to fit the optical light curve well because at earlier times Fν,max is brightest compared to the single shock model. The lower panel of Fig.  8 shows that this optical model fits the data significantly better, compared to in Fig.  5 where in order to get the best fit the assumed extinction was outside the inferred range.

Top panel: X-ray, optical, and radio data. Overlaid are the model light curves corresponding to the two jet component scenario. In this scenario, we require a range of optical extinction between 7.4 and 8.0 magnitudes. Bottom panel: the normalized residuals, the ratio of ‘Δ’ (the difference between each observed flux density and the model at that time) to σ (the uncertainty on each measured flux density).

In terms of the radio data, the narrow jet dominates the early VLA light curve. The narrow jet contribution to the measured 10 GHz emission is shown with the blue dashed line in Fig.  8. The 10 GHz narrow component light curve follows a sharp rise, |$\alpha _{\rm {10\,\, GHz}} = \frac{5}{4}$|⁠ , followed by a decay as |$\alpha _{\rm {10\,\, GHz}} =-\frac{3p}{4} = -1.7$| for p = 2.2. After the peak, the 10 GHz emission decays rapidly giving way to a wider component. Both the rise and the decay are steeper than in the single shock scenario because by the first VLA epoch the narrow jet component has already undergone a jet break. The first two VLA spectra show that the narrow component is self-absorbed. Therefore, the 1.3 GHz narrow jet would peak at around 2 |$\mu$| Jy, and therefore we do not show the 1.3 GHz contribution to the narrow component.

Despite GRB 201216C having a VHE counterpart, our modelling of the afterglow considers only synchrotron emission. Unfortunately, it is not possible to model the SSC emission for GRB 201216C given the sparse sampling in the X-ray energy range and the lack of public VHE data. However, the low value of ϵB we derived from our afterglow modelling of the narrow jet implies that the narrow jet is in a regime where SSC cooling dominates over synchrotron cooling at least at early times when there is a high fraction of electron energy lost due to radiation (Sari & Esin 2001). SSC modelling is outside the scope of this work, however, in future studies, we plan to use new tools such as detailed modelling code by Jacovich, Beniamini & Van der Horst ( 2021) which consider SSC cooling.

Our optical light curve, VLA in-band spectral indices and broad-band spectrum at 20 d, allow us to constrain the locations of νSA, νm, and νC for the narrow jet. We combine analytical models for the movement of the frequency breaks (Fig.  6) along with the decay of |${\rm F}_{\nu ,\, {\rm max}}$|⁠ , and extract physical parameters from this data set: EK,ISO,n, A*, ϵe and ϵB (assuming p = 2.2 from the X-ray data). These ranges of derived values for those parameters are given in Table  2. We also give the opening angle of the narrow component which corresponds to a jet break at 0.05 d. From the jet break, we obtain an opening angle of ∼1–9° (Chevalier & Li 2000) which corresponds to a beaming corrected opening energy between (0.01–20)× 1050 erg.

The physical parameters extracted from our data set using Peng et al. ( 2005)’s two component model. For the wider outflow, we have reduced coverage and therefore assume the same range of values for A* and use fiducial values of ϵe = 0.1 and ϵB = 0.01. We constrain p for the wider outflow from the MeerKAT-VLA spectral index measured at 54 days post burst.

The physical parameters extracted from our data set using Peng et al. ( 2005)’s two component model. For the wider outflow, we have reduced coverage and therefore assume the same range of values for A* and use fiducial values of ϵe = 0.1 and ϵB = 0.01. We constrain p for the wider outflow from the MeerKAT-VLA spectral index measured at 54 days post burst.

The stellar wind environment we infer from this data set is characterized by the parameter A*. Chevalier & Li ( 1999) relates A* to the density profile ρ = Ar−2 where A = |$\dot{M}$| /4πvw = 5 × 1011A* g cm−1. If vw = 1000 km s−1, we obtain a progenitor mass-loss rate (⁠|$\dot{\textrm {M}}$|⁠ ) of (0.6–200)|$\times 10^{-5}\, {\rm M}_{\odot }$| yr−1. Winds of massive stars are heavily dependent on metallicity and GRBs are expected to occur in low metallicity environments. If the metallicity is too high, the mass-loss rate would also be too high resulting an increased loss of angular momentum which inhibits the formation of the GRB. Furthermore, the stellar wind is expected to have a non-spherical distribution with the majority of the material concentrated around the equator. Therefore, at the poles, we expect a less distinct stellar wind profile compared to an equatorial view. The inferred range of mass-loss rate for GRB 201216C from our two-component model is within the range of expected values (Vink & de Koter 2005; Aguilera-Dena et al. 2018).

We can use the |$\dot{\textrm {M}}$| range derived from the afterglow modelling to infer limits on the progenitor mass (Langer 1989; Tramper, Sana & de Koter 2016; Yoon 2017). The upper end of our |$\dot{\textrm {M}}$| range is pushing the boundaries for the progenitor to be a Wolf–Rayet star, independent of metallicity. Low metallicities alongside with high stellar masses would be required to begin to reach such high mass losses, combined with inciting gravity waves (Fuller & Ro 2018) or if the star is reaching the Eddington luminosity (Langer et al. 1994). At the lower end of our |$\dot{\textrm {M}}$| range, the progenitor could be a Wolf–Rayet star of 10 M⊙ at solar metallicity or even 20–25 M⊙ at subsolar metallicity. Since GRBs tend to occur in low metallicity environments, the progenitor mass is likely to be between 12 and 25 M⊙.

The decay of the narrow component allows us to detect the wider outflow. The wider outflow is the origin of the observed 1.3 GHz emission: the MeerKAT light curve shows a sharp rise from 22 d, best described as emission from such a second jet component. Using the behaviour of νm and |${\rm F}_{\nu ,\, {\rm max}}$| for off-axis jets in a stellar wind environment (Chevalier & Li 1999; Peng et al. 2005), we can constrain the time at which the outflow comes into our line of sight (ton) from both the light curve and the spectra. For t< ton, νm ∝ t0, and for t > ton, |$\nu _{\textrm {m}} \propto t^{-\frac{3}{2}}$|⁠ . The dotted-dashed purple line at the bottom of Fig.  6 shows the movement of νm in the wide jet component. Comparison of the spectral index measurements between day 28 (⁠|$\beta _{1.3\, GHz} \gt -0.3$|⁠ ) and 41 (⁠|$\beta _{1.3\, GHz} = -1.1\pm 0.6$|⁠ , see also Fig.  4c) shows that there is some movement of νm between the two epochs. Therefore ton must occur before day 41. We can further constrain the deceleration time by looking at the evolution of |${\rm F}_{\nu ,\, {\rm max}}$|⁠ : for t < ton, |${\rm F}_{\nu ,\, {\rm max}}$| ∝ t3, for t > ton, |${\rm F}_{\nu ,\, {\rm max}}$| |$\propto t^{-\frac{1}{2}}$|⁠ , we place ton at ∼40 d (Peng et al. 2005).

Unlike for the narrow component, we do not have broad-band observations, and as a result we are unable to perform the same detailed modelling as presented in Section  4.2.1. We are able to determine the range of kinetic energies by assuming the same range of A* as for the narrow component and assuming that ϵe and ϵB are 0.1 and 0.01, respectively. We can calculate EK,ISO,w in the wider outflow to be (0.02–50) × 1048 erg. These values are also summarized in Table  2. We find that the isotropic equivalent kinetic energy present in the wide outflow is two and seven orders of magnitude lower than in the narrow component for the same stellar wind profile. We note that the inferred range of kinetic energies for the cocoon is dependent on the assumed values of ϵe and ϵB. With the assumed microphysical parameters, the kinetic energy of the cocoon can be considered as mildly to non-relativistic when compared to other radio transients, similar to radio-detected supernovae (e.g. fig. 5 of Coppejans et al. 2020). However the rapid rise and high luminosity are inconsistent with type Ib/c supernovae (those associated with long GRBs) and so the outflow is more likely to be mildly relativistic.

Fig.  8 shows the wide jet contribution to the 1.3 and 10 GHz light curves as the purple and blue dashed lines, respectively. We also show the total 10 GHz model from both jet components as the solid blue line. We expect the wide component of the jet to make some contribution of the total X-ray flux observed by XRT at the time of our final observation which would make the model closer to the observed flux density. We are unable to quantify the contribution of the wide jet component to the total X-ray flux as we do not know the location of νC with respect to the XRT observing band.

In comparison to our single shock scenario, the normalized residual values are much lower denoting a better fit (shown in the lower panel of Fig.  8). There are still some increased residual values early on in the 10 GHz light curve, although this may be due to weak interstellar scintillation which we cannot model. A much wider outflow makes for a much better fit to the MeerKAT light curve as well as the later VLA data points. We can better quantify whether the more complex, multiple jet component model is a better fit compared to the single forward shock model by performing an F-test. The single jet model has four parameters and the two-component model has only one additional free parameter originating from the wider component, the kinetic energy. The results show that our more complex model is favoured at greater than 4σ significance.

We have presented multiwavelength observations of the fifth VHE GRB. Spectra of the host galaxy confirmed GRB 201216C to be the highest redshift VHE GRB so far, close to the theoretical distance limit beyond which pair-production due to the EBL would prevent the detection of any VHE photons. The faint early optical detections defined this event as a dark GRB where we infer that the optical emission is attenuated by at least 5 magnitudes. Such attenuation is at least two magnitudes greater than that derived from galactic AV–NH relations. Such high extinction could be due to a high dust density in the vicinity of the lGRB, for instance if the stellar progenitor did not travel far from its formation site in the centre of a giant molecular cloud. At radio frequencies, we obtained MeerKAT (1.3 GHz), e-MERLIN (5 GHz), and VLA (10 GHz) observations covering 5 to 55 d after the burst. Our e-MERLIN epochs show evidence of DISS, allowing us to place emitting region size constraints at 29 d post burst of <1 × 1017 cm (Goodman 1997).

We interpret the data set as a whole using two possible scenarios. The first is a single forward shock, but this scenario does not explain the data well because |${\rm F}_{\nu ,\, {\rm max}}$| decays in a non-constant manner. In the first 20 d after the burst, we required a rapid decay however, later on, we needed |${\rm F}_{\nu ,\, {\rm max}}$| to be constant. A varying |${\rm F}_{\nu ,\, {\rm max}}$| decay could be a result of the jet propagating through a highly variable circumburst medium. However, one would also expect the movement of the frequency breaks to vary which is inconsistent with the inferred break evolution. We require νSA to move to lower frequencies at an unphysical rate to explain the late-time broad-band spectra, and even if such a movement were possible, the resulting MeerKAT light curve would still rise too quickly. Both the variation in |${\rm F}_{\nu ,\, {\rm max}}$| and the evolution of νSA are inconsistent with a forward shock model.

Instead, we suggest a second scenario: a jet-cocoon geometry where the earlier emission is dominated by a narrow ultrarelativistic jet which undergoes a jet break at 0.1 d. The later time emission is dominated by a wide-angled, slower moving outflow: a cocoon. We find that the additional component allows for fixing the problems regarding the rapid decay of |${\rm F}_{\nu ,\, {\rm max}}$| early on, produced by the jet break, and then the flattening occurs as the radio observations are now viewing the cocoon as a separate synchrotron component. The cocoon also addresses the unphysical movement of νSA. We find that the jet-cocoon scenario shows the afterglow to be moving through a stellar wind environment of a density similar to that modelled for massive stellar-winds with energies of (0.6–10) × 1052 and (0.02–50)× 1048 erg for the jet and cocoon, respectively. We constrain on the opening angles of the ultrarelativistic jet to be 1−9°. Deeper, more late-time observations are required moving forward in order to better understand and constrain cocoon emission in lGRB events.

The authors thank the referee for their helpful comments. LR acknowledges the support given by the Science and Technology Facilities Council through an STFC studentship. DRA-D is supported by the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of ‘Science and Society’ Action Always strive for excellence - ‘Theodoros Papazoglou’ (Project Number: 01431). This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. e-MERLIN is a National Facility operated by the University of Manchester at Jodrell Bank Observatory on behalf of STFC, part of UK Research and Innovation.

The radio data discussed in this paper are available in Table  1. The data that forms the Swift-XRT light curve is available from the Swift Burst Analyser.

https://github.com/e-merlin/e-MERLIN_CASA_pipeline

https://github.com/ska-sa/tricolour

Ackermann M. et al.  , 2013, ApJS, 209, 11 10.1088/0067-0049/209/1/11

Aguilera-Dena D. R. , Langer N. , Moriya T. J. , Schootemeijer A. , 2018, ApJ, 858, 115 10.3847/1538-4357/aabfc1

Alexander K. D. et al.  , 2017, ApJ, 848, 69 10.3847/1538-4357/aa8a76

Alexander K. D. et al.  , 2019, ApJ, 870, 67 10.3847/1538-4357/aaf19d

Beardmore A. P. et al.  , 2020, GCN Circ., 29061, 1

Belkin S. , Pozanenko A. , Krugov M. , Levkina P. , Klunko E. , Pankov N. , FuN GRB IKI , 2020, GCN Circ., 29210, 1

Bietenholz M. F. , Bartel N. , Argo M. , Dua R. , Ryder S. , Soderberg A. , 2021, ApJ, 908, 75 10.3847/1538-4357/abccd9

Blanch O. et al.  , 2020a, GCN Circ., 28659, 1

Blanch O. et al.  , 2020b, GCN Circ., 29075, 1

Campana S. , Lazzati D. , Perna R. , Grazia Bernardini M. , Nava L. , 2021, A&A, 649, A135 10.1051/0004-6361/202140439

Castro-Tirado A. J. et al.  , 2019, GCN Circ., 23708, 1

Chandra P. et al.  , 2008, ApJ, 683, 924 10.1086/589807

Chandra P. , Frail D. A. , 2012, ApJ, 746, 156 10.1088/0004-637X/746/2/156

Chen W. J. , Urata Y. , Huang K. , Takahashi S. , Petitpas G. , Asada K. , 2020, ApJ, 891, L15 10.3847/2041-8213/ab76d4

Chevalier R. A. , Li Z.-Y. , 1999, ApJ, 520, L29 10.1086/312147

Chevalier R. A. , Li Z.-Y. , 2000, ApJ, 536, 195 10.1086/308914

Coppejans D. L. et al.  , 2020, ApJ, 895, L23 10.3847/2041-8213/ab8cc7

Cordes J. M. , Lazio T. J. W. , 2002, preprint(astro–ph/0207156)

De Colle F. , Kumar P. , Aguilera-Dena D. R. , 2018, ApJ, 863, 32 10.3847/1538-4357/aad04d

Evans P. A. et al.  , 2007, A&A, 469, 379 10.1051/0004-6361:20077530

Evans P. A. et al.  , 2009, MNRAS, 397, 1177 10.1111/j.1365-2966.2009.14913.x

Evans P. A. et al.  , 2010, A&A, 519, A102 10.1051/0004-6361/201014819

Evans P. A. , Swift-XRT Team , 2021, GCN Circ., 29280, 1

Filgas R. et al.  , 2011, A&A, 526, A113 10.1051/0004-6361/201015320

Frail D. A. et al.  , 2000, ApJ, 534, 559 10.1086/308802

Frail D. A. , Kulkarni S. R. , Nicastro L. , Feroci M. , Taylor G. B. , 1997, Nature, 389, 261 10.1038/38451

Fruchter A. S. et al.  , 2006, Nature, 441, 463 10.1038/nature04787

Fuller J. , Ro S. , 2018, MNRAS, 476, 1853 10.1093/mnras/sty369

Gao H. , Lei W.-H. , Zou Y.-C. , Wu X.-F. , Zhang B. , 2013, New Astron. Rev., 57, 141 10.1016/j.newar.2013.10.001

Gokuldass P. , Morris D. , Orange N. , Strausbaugh R. , Cucchiara A. , 2021, GCN Circ., 29674, 1

Goodman J. , 1997, New Astron. Rev., 2, 449 10.1016/S1384-1076(97)00031-6

Granot J. , Sari R. , 2002, ApJ, 568, 820 10.1086/338966

Granot J. , Van der Horst A. J. , 2014, Publ. Astron. Soc. Austr., 31, e008 10.1017/pasa.2013.44

Güver T. , Özel F. , 2009, MNRAS, 400, 2050 10.1111/j.1365-2966.2009.15598.x

H. E. S. S. Collaboration , 2019, GCN Circ., 25566, 1

H. E. S. S. Collaboration 2021, Science, 372, 1081 10.1126/science.abe8560

Heywood I. , 2020, oxkat: Semi-automated imaging of MeerKAT observations, record ascl: 2009.003

Hjorth J. et al.  , 2003, Nature, 423, 847 10.1038/nature01750

Huang X.-L. , Xin L.-P. , Yi S.-X. , Zhong S.-Q. , Qiu Y.-L. , Deng J.-S. , Wei J.-Y. , Liang E.-W. , 2016, ApJ, 833, 100 10.3847/1538-4357/833/1/100

Izzo L. , Malesani D. B. , Zhu Z. P. , Xu D. , de Ugarte Postigo A. , Pursimo T. , 2020a, GCN Circ., 28661, 1

Izzo L. , Malesani D. B. , Kann D. A. , 2020b, GCN Circ., 29066, 1

Jacovich T. E. , Beniamini P. , Van der Horst A. J. , 2021, MNRAS, 504, 528 10.1093/mnras/stab911

Jakobsson P. , Hjorth J. , Fynbo J. P. U. , Watson D. , Pedersen K. , Björnsson G. , Gorosabel J. , 2004, ApJ, 617, L21 10.1086/427089

Jelinek M. et al.  , 2020, GCN Circ., 29070, 1

Kamble A. , Misra K. , Bhattacharya D. , Sagar R. , 2009, MNRAS, 394, 214 10.1111/j.1365-2966.2008.13504.x

Kent B. R. et al.  , 2018, American Astronomical Society Meeting Abstracts #231, 342.14

Kouveliotou C. , Meegan C. A. , Fishman G. J. , Bhat N. P. , Briggs M. S. , Koshut T. M. , Paciesas W. S. , Pendleton G. N. , 1993, ApJ, 413, L101 10.1086/186969

Krühler T. et al.  , 2011, A&A, 534, A108 10.1051/0004-6361/201117428

Kulkarni S. R. et al.  , 1999, ApJ, 522, L97 10.1086/312227

Langer N. , 1989, A&A, 210, 93

Langer N. , Hamann W. R. , Lennon M. , Najarro F. , Pauldrach A. W. A. , Puls J. , 1994, A&A, 290, 819

Lan M.-X. , Wu X.-F. , Dai Z.-G. , 2018, ApJ, 860, 44 10.3847/1538-4357/aac26e

Laskar T. et al.  , 2014, ApJ, 781, 1 10.1088/0004-637X/781/1/1

Laskar T. et al.  , 2016, ApJ, 833, 88 10.3847/1538-4357/833/1/88

Levan A. , Crowther P. , de Grijs R. , Langer N. , Xu D. , Yoon S.-C. , 2016, Space Sci. Rev., 202, 33 10.1007/s11214-016-0312-x

MAGIC Collaboration 2019, Nature, 575, 455 10.1038/s41586-019-1750-x

Marcote B. , Ribó M. , Paredes J. M. , Giroletti M. , Giarratana S. , Ghirlanda G. , 2020, GCN Circ., 29028, 1

Margutti R. et al.  , 2018, ApJ, 856, L18 10.3847/2041-8213/aab2ad

McMullin J. P. , Waters B. , Schiebel D. , Young W. , Golap K. , 2007, in Shaw R. A. , Hill F. , Bell D. J. , eds, ASP Conf. Ser. Vol. 376, Astronomical Data Analysis Software and Systems XVI. Astron. Soc. Pac, San Francisco, p. 127

Moldon J. , 2021, eMCP: e-MERLIN CASA pipeline, record ascl: 2109.006

Morsony B. J. , Lazzati D. , Begelman M. C. , 2007, ApJ, 665, 569 10.1086/519483

Nakar E. , Piran T. , 2017, ApJ, 834, 28 10.3847/1538-4357/834/1/28

Oates S. R. et al.  , 2011, MNRAS, 412, 561 10.1111/j.1365-2966.2010.17928.x

Oates S. R. , Beardmore A. P. , Swift/UVOT Team , 2020, GCN Circ., 29071, 1

Offringa A. R. et al.  , 2014, MNRAS, 444, 606 10.1093/mnras/stu1368

Panaitescu A. , Mészáros P. , Rees M. J. , 1998, ApJ, 503, 314 10.1086/305995

Penacchioni A. V. , Ruffini R. , Bianco C. L. , Izzo L. , Muccino M. , Pisani G. B. , Rueda J. A. , 2013, A&A, 551, A133 10.1051/0004-6361/201220679

Peng F. , Königl A. , Granot J. , 2005, ApJ, 626, 966 10.1086/430045

Perley D. A. et al.  , 2009, AJ, 138, 1690 10.1088/0004-6256/138/6/1690

Predehl P. , Schmitt J. H. M. M. , 1995, A&A, 500, 459

Racusin J. L. et al.  , 2008, Nature, 455, 183 10.1038/nature07270

Racusin J. L. et al.  , 2009, ApJ, 698, 43 10.1088/0004-637X/698/1/43

Ramirez-Ruiz E. , Celotti A. , Rees M. J. , 2002, MNRAS, 337, 1349 10.1046/j.1365-8711.2002.05995.x

Rees M. J. , Meszaros P. , 1992, MNRAS, 258, 41 10.1093/mnras/258.1.41P

Resmi L. et al.  , 2005, A&A, 440, 477 10.1051/0004-6361:20041642

Rhodes L. et al.  , 2020a, MNRAS, 496, 3326 10.1093/mnras/staa1715

Rhodes L. , Fender R. , Bray J. , Williams D. R. A. , 2020b, GCN Circ., 28945, 1

Sari R. , Esin A. A. , 2001, ApJ, 548, 787 10.1086/319003

Sari R. , Piran T. , Narayan R. , 1998, ApJ, 497, L17 10.1086/311269

Sari R. , Piran T. , Halpern J. P. , 1999, ApJ, 519, L17 10.1086/312109

Schlegel D. J. , Finkbeiner D. P. , Davis M. , 1998, ApJ, 500, 525 10.1086/305772

Shrestha M. , Melandri A. , Smith R. , Steele I. A. , Kobayashi S. , Mundell C. , Gomboc A. , Guidorzi C. , 2020, GCN Circ., 29085, 1

Solomon P. M. , Rivolo A. R. , Barrett J. , Yahil A. , 1987, ApJ, 319, 730 10.1086/165493

Starling R. L. C. , Wijers R. A. M. J. , Hughes M. A. , Tanvir N. R. , Vreeswijk P. M. , Rol E. , Salamanca I. , 2005, MNRAS, 360, 305 10.1111/j.1365-2966.2005.09042.x

Taylor G. B. , Frail D. A. , Berger E. , Kulkarni S. R. , 2004, ApJ, 609, L1 10.1086/422554

Tramper F. , Sana H. , de Koter A. , 2016, ApJ, 833, 133 10.3847/1538-4357/833/2/133

Valeev A. F. et al.  , 2019, GCN Circ., 25565, 1

Van der Horst A. J. et al.  , 2008, A&A, 480, 35 10.1051/0004-6361:20078051

Van der Horst A. J. et al.  , 2014, MNRAS, 444, 3151 10.1093/mnras/stu1664

Van der Horst A. J. , Kouveliotou C. , Gehrels N. , Rol E. , Wijers R. A. M. J. , Cannizzo J. K. , Racusin J. , Burrows D. N. , 2009, ApJ, 699, 1087 10.1088/0004-637X/699/2/1087

Vaughan S. , Edelson R. , Warwick R. S. , Uttley P. , 2003, MNRAS, 345, 1271 10.1046/j.1365-2966.2003.07042.x

Vielfaure J. B. et al.  , 2020, GCN Circ., 29077, 1

Vink J. S. , de Koter A. , 2005, A&A, 442, 587 10.1051/0004-6361:20052862

Vreeswijk P. M. et al.  , 2018, GCN Circ., 22996, 1

Walker M. A. , 1998, MNRAS, 294, 307 10.1046/j.1365-8711.1998.01238.x

Wijers R. A. M. J. , Galama T. J. , 1999, ApJ, 523, 177 10.1086/307705

Yoon S.-C. , 2017, MNRAS, 470, 3970 10.1093/mnras/stx1496

Zauderer B. A. et al.  , 2013, ApJ, 767, 161 10.1088/0004-637X/767/2/161

Zhang W. , Woosley S. E. , Heger A. , 2004, ApJ, 608, 365 10.1086/386300

Zhang L.-L. , Ren J. , Huang X.-L. , Liang Y.-F. , Lin D.-B. , Liang E.-W. , 2021, ApJ, 917, 95 10.3847/1538-4357/ac0c7f

Oxford University Press is a department of the University of Oxford. It furthers the University's objective of excellence in research, scholarship, and education by publishing worldwide

Sign In or Create an Account

This PDF is available to Subscribers Only

For full access to this pdf, sign in to an existing account, or purchase an annual subscription.