A Reverse Shock in GRB 181201A

The Swift X-ray Telescope (XRT; Burrows et al. 2005) began observing GRB 181201A 0.16 days after the burst and detected a fading X-ray afterglow at R.A. = 21^h17^m1120, decl. = −12°37' 509 (J2000), with an uncertainty radius of 14 (90% containment).²⁰ XRT continued observing the afterglow for 20 days in photon-counting mode. The GRB entered a Sun constraint for Swift at ≈21 days after the burst and was dropped from the Swift schedule. We obtained 4.4 ks of Swift ToO observations on 2019 April 7 (≈127 days), when it was visible again. The X-ray afterglow was weakly detected. The data were automatically calibrated using the XRT pipeline and incorporated into the online light curve. We use a photon index of Γ_X ≈ 1.77 and an unabsorbed counts-to-flux conversion rate from the Swift website of 4.75 × 10⁻¹¹ erg ${\mathrm{cm}}^{-2}$ ct⁻¹ to convert the 0.3–10 keV count rate light curve²¹ to flux density at 1 keV for subsequent analysis.

We derived aperture photometry from pipeline-processed UVOT images downloaded from the Swift data center²² with the uvotproduct (v2.4) software. We used a circular source extraction region with a 35 radius, followed by an aperture correction based on the standard 5'' aperture for UVOT photometry (Poole et al. 2008). We present our UVOT measurements in Table 1.

Four MASTER-net observatories (MASTER-Kislovodsk, MASTER-IAC, MASTER-SAAO, and MASTER-OAFA) participated in rapid-response and follow-up observations of GRB 181201A, beginning with 10 s exposures at 02:40:30 UT on 2019 December 1, ≈150 s after the burst and 10 s after the INTEGRAL alert in the BVRI and Clear bands²³ (Podesta et al. 2018). The MASTER autodetection system (Lipunov et al. 2010, 2019) discovered a bright optical afterglow at R.A. = 21^h17^m1120, decl. = −12°37' 514 (J2000) with unfiltered magnitude ≈13.2 mag (Vega). MASTER follow-up observations continued until ≈3 days after the burst. We carried out astrometry and aperture photometry (using a ≈25 aperture radius) for observations taken at each telescope separately, calibrated to six nearby SDSS reference stars. We verified our photometry by measuring the brightness of comparison stars with similar brightness to the afterglow, and incorporated the derived systematic calibration uncertainty into the photometric uncertainty (Troja et al. 2017; Lipunov et al. 2019).

We observed the afterglow of GRB 181201A with the Gamma-Ray Burst Optical/Near-Infrared Detector (GROND) on the MPI/ESO 2.2 m telescope in La Silla in Chile beginning 0.91 days after the burst in g'r'i'z'JHK filters. We calibrated the data with tools and methods standard for GROND observations and detailed in Krühler et al. (2008). We performed aperture photometry calibrated against PanSTARRS field stars in the g'r'i'z' filters and against the 2MASS catalog in JHK.

We observed the afterglow with the Ohio State Multi-Object Spectrograph (OSMOS; Martini et al. 2011) on the 2.4 m Hiltner telescope at MDM Observatory on 2018 December 4, with 120 s exposures each in the BVRI bands. We obtained two later epochs of imaging using the Templeton detector on the 1.3 m McGraw-Hill telescope on the nights of 2018 December 9 and 10, obtaining $4\times 300\,{\rm{s}}$ each in r' and i' bands. Observations on December 10 were taken at a high airmass of 1.91 and under poor seeing. We performed aperture photometry with a 2'' aperture on the Hiltner and McGraw-Hill images, referenced to SDSS.

We obtained Liverpool Telescope (LT; Steele et al. 2004) imaging with the infrared/optical camera (IO:O) at ≈2.7 days, comprising 3 × 120 s exposures in g'r'i' bands. We performed aperture photometry using the Astropy Photutils package (Bradley et al. 2016) calibrated to SDSS. We also collected all optical photometry for this event published in GCN circulars. We present our combined optical and NIR data set in Table 2.

We obtained five epochs of ALMA Band 3 (3 mm) observations of GRB 181201A spanning 0.88 days to 29.8 days after the burst through program 2018.1.01454.T (PI: Laskar). The observations utilized two 4 GHz wide basebands centered at 91.5 and 103.5 GHz, respectively, in the first epoch, and at 90.4 and 102.4 GHz in subsequent epochs. One of J2258−2758, J2148+0657, and J2000−1748 was used for bandpass and flux density calibration, while J2131−1207 was used for complex gain calibration in all epochs. We analyzed the data using the Common Astronomy Software Application (CASA; McMullin et al. 2007), and the afterglow was well detected in all epochs. We performed one round of phase-only self-calibration in the first two epochs, where the deconvolution residuals indicated the presence of phase errors. We list the results of our ALMA observations in Table 3.

We observed the afterglow using the VLA starting 2.9 days after the burst through programs 18A-088 and 18B-407 (PI: Laskar). We detected and tracked the flux density of the afterglow from 1.2 GHz to 37 GHz over multiple epochs until ≈164 days after the burst. We used 3C 48 as the flux and bandpass calibrator and J2131−1207 as the complex gain calibrator. Where phase errors appeared to be present after deconvolution, we performed phase-only self-calibration in epochs with sufficient signal-to-noise. The effect of self-calibration impacted the measured flux density by ≲15% and only at the highest frequencies (≳15 GHz). We carried out data reduction using CASA and list the results of our VLA observations in Table 3.

We observed the afterglow using the upgraded Giant Meterwave Radio Telescope (uGMRT) through program 35_065 (PI: Laskar) starting 12.5 days and 13.5 days after the burst in Bands 4 (center frequency 550 MHz) and 5 (center frequency 1450 MHz), respectively. The observations utilized the full 400 MHz bandwidth of the uGMRT Wide-band Backend (GWB) system. We used 3C 48 for bandpass and flux density calibration, and J2131−1207 for complex gain calibration. We carried out data reduction using CASA and list the results of our uGMRT observations in Table 3.

One of the most striking features of GRB 181201A revealed by our VLA and ALMA observations is a spectral bump in the centimeter- to millimeter-band SED at 3.9 days (Figure 1). The feature spans 3–24 GHz and can be fit with a smoothly broken power law, with low-frequency index²⁴ β₁ = 1.43 ± 0.23, high-frequency index β₂ = −0.35 ± 0.03, break frequency ν_peak = 5.5 ± 0.4 GHz, and peak flux density F_ν,peak = 1.48 ±0.06 mJy. The steeply rising spectrum below the spectral peak suggests the underlying emission is self-absorbed.²⁵ Fitting the VLA centimeter-band data above 24 GHz and the ALMA observations at the same epoch with a single power-law model yields ${\beta }_{24-\mathrm{ALMA}}=0.09\pm 0.02$.

Zoom In Zoom Out Reset image size

The radio and millimeter SED at 3.9 days clearly exhibits two distinct components. The simplest explanation for multiple components in the radio is to attribute one to an RS and the other to an FS. Because the RS emission is expected to peak at lower frequencies than the FS and the RS spectrum is expected to cut off steeply above the RS cooling frequency (Kobayashi 2000), the lower-peaking component is more naturally associated with RS emission.

The ALMA millimeter-band light curve can be fit with a single power law, with decay rate α_ALMA = −0.71 ± 0.01 (Figure 2). From the radio SED at 3.9 days, we know that ${\nu }_{{\rm{m}},{\rm{f}}}\gt {\nu }_{\mathrm{ALMA}}$ at this time. In this spectral regime, the expected rise rate for FS emission is t^1/2 in the ISM environment and t⁰ in the wind environment, both of which are incompatible with the observed millimeter-band light curve. Thus, the millimeter-band light curve at <4 days originates from a separate component. The shallowest decline rate for an RS light curve at $\nu \gtrsim {\nu }_{{\rm{m}},{\rm{r}}}$ is α ≈ −2.1 for a relativistic RS and α − 1.5 for a nonrelativistic RS, in each case corresponding to a wind environment with p ≈ 2. Both of these are steeper than the observed millimeter-band light curve. However, the sum of the contributions from the two emission components (RS and FS) may yield the observed shallower decay. We investigate this further in Section 5.2.

Zoom In Zoom Out Reset image size

For the remainder of this article, we associate the low- and high-frequency components with emission from the RS and the FS, respectively, and consider other scenarios in Section 5.5. In this framework, our observations of GRB 181201A provide the first instance where a single radio SED can be instantaneously and cleanly decomposed into RS and FS components, with the RS dominating at the lower frequency bands and the underlying FS emission clearly visible as an optically thin, rising spectrum at higher frequencies. We note that a more detailed characterization of these components, including the determination of the precise locations of their spectral break frequencies, requires a multiband model fit, which we undertake in Section 4.

The identification of the FS component now allows us to constrain the circumburst density profile. We note that there is no evidence for a jet break until ≈21 days, as the X-ray light curve declines as a single power law from ≈0.15 days to ≈21 days (α_X = −1.22 ± 0.01; Figure 3). The rising radio spectrum to ≈90 GHz at 3.9 days indicates that ${\nu }_{{\rm{m}},{\rm{f}}}\gtrsim 90\,\mathrm{GHz}$ at this time. The relatively shallow spectral index of β ≈ 0.1 compared to the expected index of β = 1/3 may result from the smoothness and potential proximity of the ${\nu }_{{\rm{m}},{\rm{f}}}$ break (Granot & Sari 2002).

Zoom In Zoom Out Reset image size

Because ${\nu }_{{\rm{m}},{\rm{f}}}\gtrsim 90\,\mathrm{GHz}$ at 3.9 days, the FS peak flux density ${F}_{\nu ,{\rm{m}},{\rm{f}}}\gtrsim {F}_{3\mathrm{mm}}(3.9\,\mathrm{days})\approx 1.2$ mJy. In the ISM environment, this peak flux density, F_p, remains fixed and must appear at a lower centimeter-band frequency, ν_cm, at a later time given by

$$\begin{eqnarray}&&{t}_{{\rm{p}}}={\left(\displaystyle \frac{{\nu }_{\mathrm{ALMA}}}{{\nu }_{\mathrm{cm}}}\right)}^{2/3}{\left(\displaystyle \frac{{F}_{{\rm{p}}}}{{F}_{\mathrm{ALMA}}}\right)}^{2}.\end{eqnarray} \tag{ 1 }$$

This implies a flux density of F_p ≳ 1.2 mJy at t_p ≳ 13.2 days at 16 GHz, which is a factor of 2 brighter than the observed 16 GHz light curve at the corresponding time. Indeed, all centimeter-band light curves are inconsistent with this extrapolation, and thus inconsistent with an ISM-like environment.

We now show that the optical to X-ray data are also inconsistent with an ISM environment and indicate a wind-like circumburst medium. The GROND SED at 0.91 days corrected for a Galactic extinction of A_V = 0.17 mag (Schlafly & Finkbeiner 2011) can be fit with a single power law, with index β_opt = −0.49 ± 0.04 (Figure 4). On the other hand, the spectral index between the GROND K-band and the X-rays is steeper, ${\beta }_{\mathrm{NIR},{\rm{X}}}=-0.75\pm 0.02$, indicating that a spectral break lies between the optical/NIR and X-rays. The most natural explanation for these observations is ${\nu }_{{\rm{m}},{\rm{f}}}\lt {\nu }_{\mathrm{opt}}\,\lt {\nu }_{{\rm{c}},{\rm{f}}}\lt {\nu }_{{\rm{X}}}$. Fitting the optical and X-ray SED at 0.9 days, 2.9 days, and 4.9 days with broken power-law models²⁶ fixing the low-frequency index to β_opt = −0.5 and the high-frequency index to β_X = −1.0, we find a break frequency ν_b =(5.7 ± 1.5) × 10¹⁵ Hz, (8.1 ± 2.0) × 10¹⁵ Hz, and (9.1 ±3.3) × 10¹⁵ Hz, respectively (Figure 4). The rising break frequency is roughly consistent with the evolution of ${\nu }_{{\rm{c}},{\rm{f}}}\propto {t}^{1/2}$ in a wind-like environment. We caution that the inferred location of the break is only indicative, as it is degenerate with the spectral slopes and break smoothness.

Zoom In Zoom Out Reset image size

For ${\nu }_{{\rm{m}},{\rm{f}}}\lt {\nu }_{\mathrm{opt}}\lt {\nu }_{{\rm{c}}}$, we expect a spectral index of β = (1−p)/2. Combined with the observed GROND SED at 0.9 days, this implies p ≈ 2, which yields a light-curve decay rate of α_opt ≈ −1.25 for the wind environment and α_opt ≈ −0.75 for the ISM environment. The observed decay rate in the Swift/UVOT uvw2-band and in the ground-based optical r'-band observations at ≈0.2–8.7 days is α_UV = −1.21 ± 0.02, consistent with the wind case and inconsistent with an ISM environment (Figure 3).

Above ${\nu }_{{\rm{c}},{\rm{f}}}$, we expect a decay rate of α = (2–3p)/4 ≈ −1. A power-law fit to the X-ray light curve yields a decay rate of α_X = −1.212 ± 0.013; however, the light curve after ≈6 days is consistently above the single power-law fit. Adding a temporal break at t_b ≈ 0.4 days yields a better fit with α_X,1 = −1.7 ± 0.1 and α_X,2 = −0.96 ± 0.02. The steep early decay may imply that the cooling break is not far below the X-ray band, and we verify this with our full multiwavelength model in Section 4.1. On the other hand, the late-time excess is unusual, and we discuss this further in Section 5.3.

The measured X-ray photon index of ${{\rm{\Gamma }}}_{{\rm{X}}}=1.77\pm 0.06$ implies β_X = −0.77 ± 0.06, which is not consistent with the β = −p/2 ≈ −1 expected for p ≈ 2. We note that the Klein–Nishina (KN) correction may yield a spectral index of β = −3(p − 1)/4 ≈ −0.75 in some regimes (Nakar et al. 2009), and it is possible that this may impact the X-ray spectrum but not the light curve.²⁷ We defer a detailed discussion of the KN corrections to future work.

To summarize, the radio SED at 3.9 days and the millimeter-band light curve indicate the presence of an emission component separate from that responsible for the X-ray and optical emission. We interpret this low-frequency component as RS emission. Furthermore, we have described several observations that favor a wind-like circumburst environment: (i) the high millimeter-band flux density at 4 days that does not appear in the lower-frequency radio light curve, (ii) a slowly rising break between the optical and X-rays, and (iii) the optical light curve, which is consistent with ${\nu }_{{\rm{m}},{\rm{f}}}\lt {\nu }_{\mathrm{opt}}\lt {\nu }_{{\rm{c}}}$ in a wind environment. We therefore focus on the wind environment in the remainder of this article.

The first two MASTER observations at 156 and 209 s after the burst (156 s and 209 s, respectively) reveal a rapid brightening, α ≈ 2.8. Such brightenings have only been seen in a few other optical afterglows (Molinari et al. 2007; Liang et al. 2010; Melandri et al. 2010) and have been interpreted as the onset of the afterglow, where the increase in flux density corresponds to an ongoing transfer of energy from the ejecta to the FS, ultimately resulting in the deceleration of the ejecta (Rees & Mészáros 1992; Mészáros & Rees 1993; Kobayashi et al. 1999; Sari & Piran 1999b; Kobayashi & Zhang 2007). Unfortunately, our data do not capture the time of the peak in the optical bands. As the FS is not fully set up during this onset period, we do not include these points in our multiwavelength modeling (Section 4). We discuss these data further in Section 5.1.

Guided by the basic considerations discussed above, we fit all available X-ray to radio observations with a model comprising an FS and RS produced by a relativistic jet propagating in a wind-like circumburst environment. The free parameters for the FS model are p, ${E}_{{\rm{K}},\mathrm{iso},52}$, A_*, ${\epsilon }_{{\rm{e}}}$, ${\epsilon }_{{\rm{B}}}$, and the extinction ${A}_{{\rm{V}}}$. Given the relatively shallow decline rate of the millimeter-band light curve (α ≈ −0.72) corresponding to the RS component, we focus on the nonrelativistic RS model, which allows for a shallower light-curve decay than a relativistic RS. The free parameters for the RS model²⁸ are the break frequencies ${\nu }_{{\rm{a}},{\rm{r}}}$, ${\nu }_{{\rm{m}},{\rm{r}}}$, ${\nu }_{{\rm{c}},{\rm{r}}}$ and the RS peak flux, ${F}_{\nu ,{\rm{m}},{\rm{r}}}$ all at a reference time ${t}_{\mathrm{ref}}\equiv 0.01$ days selected to be ≈few × T₉₀ (and hence comparable to ${t}_{\mathrm{dec}}$). Additionally, we allow the parameter g, corresponding to the evolution of the RS Lorentz factor with radius after the deceleration time (Γ ∝ R^−g), to vary (Mészáros & Rees 1999; Kobayashi & Sari 2000). Finally, we allow for a constant contribution to the light curves from an underlying host galaxy in bands where late-time flattening is evident (optical griz bands and the centimeter bands from 2.7 to 24.5 GHz). All free parameters, including RS and FS components and the underlying host galaxy contamination, are fit simultaneously. Further details of the modeling procedure, including the likelihood function employed, are described elsewhere (Laskar et al. 2013, 2014, 2016).

We performed a Markov Chain Monte Carlo (MCMC) exploration of the multidimensional parameter space using emcee (Foreman-Mackey et al. 2013). We ran 512 walkers for 10,000 steps and discarded a few (≲0.15%) initial nonconverged steps with low log likelihood as burn-in. Our highest likelihood model with p ≈ 2.1, ${\epsilon }_{{\rm{e}}}\approx 0.37$ ${\epsilon }_{{\rm{B}}}\approx 9.6\times {10}^{-3}$, E_K,iso ≈ 2.2 × 10⁵³ erg, ${A}_{* }\approx 1.9\times {10}^{-2}$, and negligible extinction, ${A}_{{\rm{V}}}\lesssim 0.02$ mag, together with an RS model with g ≈ 2.7, ${\nu }_{{\rm{a}},{\rm{r}}}\approx 4.2\,\mathrm{GHz}$, ${\nu }_{{\rm{c}},{\rm{r}}}\approx 90\,\mathrm{GHz}$, and F_ν,a,r ≈ 0.6 mJy at 3.9 days fits the X-ray to radio light curves and radio SEDs well (Figures 5 and 6). The peak in the centimeter-band spectrum at 3.9 days corresponds to ${\nu }_{{\rm{a}},{\rm{r}}}$. The sum of the RS self-absorbed spectrum and the FS optically thin ν^1/3 spectrum below ${F}_{\nu ,{\rm{a}},{\rm{r}}}$ (Figure 6) explains the intermediate β ≈ 1.4 spectral index below the centimeter-band peak (Section 3.1). We present a selection of SEDs spanning radio to X-rays in Figure 7, and all light curves together in Figure 8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For this model, the FS cooling frequency ${\nu }_{{\rm{c}},{\rm{f}}}\approx 2.6\,\times {10}^{16}$ Hz at 1 day lies between the optical and X-rays, as expected from the basic considerations in Section 3.2. The proximity of ${\nu }_{{\rm{c}},{\rm{f}}}$ to the X-ray band (${\nu }_{{\rm{c}},{\rm{f}}}\approx 0.3\,\mathrm{keV}$ at 8 days) explains the moderately steep decline in the X-ray light curve, compared to the expectation for ${\nu }_{{\rm{X}}}\gtrsim {\nu }_{{\rm{c}},{\rm{f}}}$. Furthermore, ${\nu }_{{\rm{a}},{\rm{f}}}\approx 60\,\mathrm{MHz}$ is below the VLA and GMRT frequencies and is therefore not constrained, resulting in degeneracies between the model parameters (Figure 9). The FS characteristic frequency at 4 days, ${\nu }_{{\rm{m}},{\rm{f}}}\approx 6\times {10}^{11}$ Hz, is above the ALMA band, as expected from the radio/millimeter SED at that time (Section 3.1). We list the properties of the best-fit model and summarize the marginalized posterior density function in Table 4.

Zoom In Zoom Out Reset image size

Table 4.  FS and RS Model Parameters

Parameter	Best Fit	MCMCa
Forward Shock
p	2.11	2.11 ± 0.01
$\mathrm{log}{\epsilon }_{{\rm{e}}}$	−0.43	$-{0.39}_{-0.17}^{+0.12}$
$\mathrm{log}{\epsilon }_{{\rm{B}}}$	−2.01	$-{2.20}_{-0.78}^{+0.42}$
$\mathrm{log}{A}_{* }$	−1.72	$-{1.66}_{-0.13}^{+0.23}$
$\mathrm{log}({E}_{{\rm{K}},\mathrm{iso}}/$ erg)	53.3	${53.41}_{-0.16}^{+0.13}$
$\mathrm{log}({A}_{{\rm{V}}}/$ mag)	−1.83	$-{1.84}_{-0.18}^{+0.14}$
$\mathrm{log}{\nu }_{{\rm{a}},{\rm{f}}}$ (Hz)	7.74b	⋯
$\mathrm{log}{\nu }_{{\rm{m}},{\rm{f}}}$ (Hz)	12.7	⋯
$\mathrm{log}{\nu }_{{\rm{c}},{\rm{f}}}$ (Hz)	16.4	⋯
$\mathrm{log}{F}_{\nu ,{\rm{m}},{\rm{f}}}$ (mJy)	0.62	⋯
Reverse Shock
$\mathrm{log}({\nu }_{{\rm{a}},{\rm{r}}}/$ Hz)	12.1	12.11 ± 0.06
$\mathrm{log}({\nu }_{{\rm{m}},{\rm{r}}}/$ Hz)	8.6	${8.64}_{-0.45}^{+1.01}$
$\mathrm{log}({\nu }_{{\rm{c}},{\rm{r}}}/$ Hz)	14.7	${14.65}_{-0.22}^{+0.25}$
$\mathrm{log}({F}_{\nu ,{\rm{m}},{\rm{r}}}/$ mJy)	5.2	${5.16}_{-0.56}^{+0.26}$
g	2.72	${2.74}_{-0.31}^{+0.38}$

Notes. Frequencies and flux densities are calculated at 1 day (FS) and 10⁻² days (RS).

^aSummary statistics from the marginalized posterior density distributions, with median and ±34.1% quantiles (corresponding to ±1σ for Gaussian distributions; Figure 9). ^bThis frequency is not directly constrained by the data.

Download table as:  ASCIITypeset image

We now use the MCMC results to derive the initial Lorentz factor of the jet. The jet's initial Lorentz factor can be determined from joint RS and FS observations by solving the following equations self-consistently (Zhang et al. 2003; Laskar et al. 2018a) at the deceleration time, ${t}_{\mathrm{dec}}$,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\nu }_{{\rm{m}},{\rm{r}}}}{{\nu }_{{\rm{m}},{\rm{f}}}} & \sim & {{\rm{\Gamma }}}_{0}^{-2}{R}_{{\rm{B}}},\\ \displaystyle \frac{{\nu }_{{\rm{c}},{\rm{r}}}}{{\nu }_{{\rm{c}},{\rm{f}}}} & \sim & {R}_{{\rm{B}}}^{-3},\\ \displaystyle \frac{{F}_{\nu ,{\rm{m}},{\rm{r}}}}{{F}_{\nu ,{\rm{m}},{\rm{f}}}} & \sim & {{\rm{\Gamma }}}_{0}{R}_{{\rm{B}}},\end{array}\end{eqnarray} \tag{ 2 }$$

where ${{\rm{\Gamma }}}_{0}$ is the initial Lorentz factor of the GRB jet, ${R}_{{\rm{B}}}\equiv {({\epsilon }_{{\rm{B}},{\rm{r}}}/{\epsilon }_{{\rm{B}},{\rm{f}}})}^{1/2}$ is the relative magnetization of the ejecta, and we have assumed similar values of ${\epsilon }_{{\rm{e}}}$ and the Compton Y parameter in the two regions. The effective magnetizations of the shocked ejecta and shock external medium may well be very different, and the introduction of ${R}_{{\rm{B}}}$ allows us to investigate such variations caused, for instance, by the presence of magnetic fields produced in the initial GRB and advected in the outflow (Usov 1992; Mészáros & Rees 1997; Zhang et al. 2003).

This system of equations has three unknowns (${t}_{\mathrm{dec}}$, ${R}_{{\rm{B}}}$, and Γ₀) and can be inverted exactly to yield

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{dec}} & = & {t}_{\mathrm{ref}}{\left[\left(\displaystyle \frac{{\nu }_{{\rm{c}},{\rm{f}}}^{{\prime} }}{{\nu }_{{\rm{c}},{\rm{r}}}^{{\prime} }}\right){\left(\displaystyle \frac{{F}_{\nu ,{\rm{m}},{\rm{f}}}^{{\prime} }}{{F}_{\nu ,{\rm{m}},{\rm{r}}}^{{\prime} }}\right)}^{2}\left(\displaystyle \frac{{\nu }_{{\rm{m}},{\rm{f}}}^{{\prime} }}{{\nu }_{{\rm{m}},{\rm{r}}}^{{\prime} }}\right)\right]}^{\tau },\\ \tau & \equiv & {\left[\displaystyle \frac{15}{2(4-k)}-\displaystyle \frac{8(g+3)}{7(2g+1)}\right]}^{-1}\end{array}\end{eqnarray} \tag{ 3 }$$

where the primed quantities refer to values at the reference time of ${t}_{\mathrm{ref}}=0.01$ days used in the MCMC analysis. For our best-fit model, this yields ${t}_{\mathrm{dec}}\approx 5.3\times {10}^{-3}$ days or ≈460 s; the full MCMC analysis gives $\mathrm{log}({t}_{\mathrm{dec}}/\mathrm{days})=-{2.27}_{-0.14}^{+0.13}$. We note that this is longer than the burst duration, T₉₀ ≈ 172 s (Section 2), consistent with the thin shell case (Zhang et al. 2003).

Given this value of ${t}_{\mathrm{dec}}$, which satisfies Equation (2), we can in principle compute the initial Lorentz factor and magnetization,

$$\begin{eqnarray}\begin{array}{rcl}{R}_{{\rm{B}}} & = & {\left(\displaystyle \frac{{\nu }_{{\rm{c}},{\rm{f}}}}{{\nu }_{{\rm{c}},{\rm{r}}}}\right)}^{1/3}\\ {{\rm{\Gamma }}}_{0} & = & \displaystyle \frac{{F}_{\nu ,{\rm{m}},{\rm{r}}}}{{F}_{\nu ,{\rm{m}},{\rm{f}}}}{R}_{{\rm{B}}}^{-1}.\end{array}\end{eqnarray} \tag{ 4 }$$

However, we find that whereas ${\nu }_{{\rm{c}},{\rm{r}}}$ and ${\nu }_{{\rm{c}},{\rm{f}}}$ are reasonably well constrained, ${F}_{\nu ,{\rm{m}},{\rm{r}}}$ is completely degenerate with ${\nu }_{{\rm{m}},{\rm{r}}}$. This is because the RS peak in the radio band corresponds to ${\nu }_{{\rm{a}},{\rm{r}}}$, and hence the data do not constrain ${\nu }_{{\rm{m}},{\rm{r}}}$. In the observed regime of ${\nu }_{{\rm{m}},{\rm{r}}}\lt {\nu }_{{\rm{a}},{\rm{r}}}$, the peak flux density of the RS component, ${F}_{\nu ,{\rm{a}},{\rm{r}}}\propto {F}_{\nu ,{\rm{m}},{\rm{r}}}{({\nu }_{{\rm{a}},{\rm{r}}}/{\nu }_{{\rm{m}},{\rm{r}}})}^{(1-p)/2}$ results in the degeneracy, ${F}_{\nu ,{\rm{m}},{\rm{r}}}\propto {\nu }_{{\rm{m}},{\rm{r}}}^{(1-p)/2}\propto {\nu }_{{\rm{m}},{\rm{r}}}^{-0.5}$ for p ≈ 2, which exactly explains the observed correlation between these parameters (Figure 9). Our MCMC samples span four orders of magnitude in ${\nu }_{{\rm{m}},{\rm{r}}}$ and correspondingly about two orders of magnitude in ${F}_{\nu ,{\rm{m}},{\rm{r}}}$. The sampling is clipped at the lower end in ${\nu }_{{\rm{m}},{\rm{r}}}$ by the prior on ${\nu }_{{\rm{m}},{\rm{r}}}\gt {10}^{8}\,\mathrm{Hz}$ (corresponding to ${F}_{\nu ,{\rm{m}},{\rm{r}}}\lesssim {10}^{8}$ mJy) and at the lower end in ${F}_{\nu ,{\rm{m}},{\rm{r}}}$ by the actual observed peak flux density of the RS component (corresponding to ${\nu }_{{\rm{m}},{\rm{r}}}\lesssim {10}^{12}$ Hz). This degeneracy between ${\nu }_{{\rm{m}},{\rm{r}}}$ and ${F}_{\nu ,{\rm{m}},{\rm{r}}}$ does not affect the calculation of ${t}_{\mathrm{dec}}$ (Equation (3)), because ${\nu }_{{\rm{m}},{\rm{r}}}$ and ${F}_{\nu ,{\rm{m}},{\rm{r}}}$ enter in the combination ${F}_{\nu ,{\rm{m}},{\rm{r}}}{\nu }_{{\rm{m}},{\rm{r}}}^{2}\propto {\nu }_{{\rm{m}},{\rm{r}}}^{2-p}$, which is only weakly dependent on ${\nu }_{{\rm{m}},{\rm{r}}}$ for p ≈ 2 (i.e., the degenerate terms cancel out). Therefore, we can constrain ${t}_{\mathrm{dec}}$ and ${R}_{{\rm{B}}}$ through this system of equations, but not ${{\rm{\Gamma }}}_{0}$. We find ${R}_{{\rm{B}}}\approx 0.65$ for the best-fit solution, and ${R}_{{\rm{B}}}={0.64}_{-0.08}^{+0.07}$ from the MCMC.

To infer ${{\rm{\Gamma }}}_{0}$ in this case, we turn to an energy balance argument (Zhang et al. 2003). At ${t}_{\mathrm{dec}}$, the energy in the swept-up material (E_sw) is similar to the rest energy of the ejecta (which is the same as the explosion kinetic energy, E),

$$\begin{eqnarray}&&{E}_{\mathrm{sw}}={{\rm{\Gamma }}}_{0}^{2}{M}_{\mathrm{sw}}{c}^{2}\approx E,\end{eqnarray} \tag{ 5 }$$

where the swept-up mass at radius, R,

$$\begin{eqnarray}&&{M}_{\mathrm{sw}}=\displaystyle \frac{4\pi {{AR}}^{3-k}}{3-k}\end{eqnarray} \tag{ 6 }$$

Combining Equations (5) and (6), the deceleration radius

$$\begin{eqnarray}&&R={\left[\displaystyle \frac{(3-k)E}{4\pi A{{\rm{\Gamma }}}_{0}^{2}{c}^{2}}\right]}^{1/(3-k)}={{\rm{\Gamma }}}_{0}^{-2/(3-k)}{l}_{\mathrm{Sedov}},\end{eqnarray} \tag{ 7 }$$

where ${l}_{\mathrm{Sedov}}$ is the Sedov length, which depends only on the energy and density, parameters that are known from the FS modeling. For a point source moving along the line of sight, this radius corresponds to the deceleration time given by $R=2{{\rm{\Gamma }}}_{0}^{2}{{ct}}_{\mathrm{dec}}/(1+z)$, which gives the initial Lorentz factor

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{0}\approx {\left[\displaystyle \frac{{l}_{\mathrm{Sedov}}}{c}\left(\displaystyle \frac{2{t}_{\mathrm{dec}}}{1+z}\right)\right]}^{\tfrac{3-k}{2(4-k)}}.\end{eqnarray} \tag{ 8 }$$

We note that this is a simplified treatment, and a more rigorous analysis requires integrating over the full fluid profile and using the specific enthalpy of the shocked region while calculating the swept-up mass. However, for our purposes, this approximate expression will suffice to estimate ${{\rm{\Gamma }}}_{0}$. We compute ${{\rm{\Gamma }}}_{0}$ for all MCMC samples using Equation (8), resulting in ${{\rm{\Gamma }}}_{0}\,={103}_{-8}^{+10}$ (Figure 10).

Zoom In Zoom Out Reset image size

We now return to the optical rise observed at ≲2.4 × 10⁻³ days with MASTER (Figure 5). These observations were performed at $\lesssim {t}_{\mathrm{dec}}$ in our best-fit RS+FS model (Section 4 and Table 4), consistent with the afterglow onset scenario. The spectral ordering at ${t}_{\mathrm{dec}}$ is ${\nu }_{{\rm{m}},{\rm{r}}}\lt {\nu }_{{\rm{a}},{\rm{r}}}\lt {\nu }_{\mathrm{opt}}\approx {\nu }_{{\rm{c}},{\rm{r}}}$ and ${\nu }_{{\rm{a}},{\rm{f}}}\lt {\nu }_{{\rm{c}},{\rm{f}}}\lt {\nu }_{\mathrm{opt}}\lt {\nu }_{{\rm{m}},{\rm{f}}}$. In this model, the RS light curve at $\lesssim {t}_{\mathrm{dec}}$ is expected to rise at the rate α = (p − 1)/2 ≈ 0.5, while the FS light curve should rise as α = 1/2 (Zou et al. 2005). However, the observed rise of α ≈ 2.8 is not consistent with either prediction. The proximity of ${\nu }_{{\rm{c}},{\rm{r}}}$ to ${\nu }_{\mathrm{opt}}$ at the deceleration time may suppress the RS flux. On the other hand, the FS model also overpredicts the observed flux. One possible solution is the rapid injection of energy into the FS during this period. We find that an energy increase at the rate E ∝ t^1.8 until ${t}_{\mathrm{dec}}$ roughly accounts for the optical rise²⁹ (Figure 11). Such an injection could arise, for instance, by the presence of energy at Lorentz factors $\lesssim {{\rm{\Gamma }}}_{0}$ owing to velocity stratification in the ejecta (Sari & Mészáros 2000), which has previously been inferred in other GRB afterglows (Laskar et al. 2015, 2018a, 2018c). An alternate scenario could be a pair-dominated optical flash, as suggested for GRB 130427A (Vurm et al. 2014). However, the paucity of data around the optical peak precludes a more detailed analysis.

Zoom In Zoom Out Reset image size

The 3 mm model light curve underpredicts two of the ALMA detections at ≈1.9 and ≈3.9 days by ≈20% (Figure 12). In particular, the observed light curve exhibits a shallower decline than afforded by the current model, where the transition between RS dominant to FS dominant at ≈1.5 days results in an inflection in the model light curve, which is not seen in the data. One possibility for a shallower RS decline rate is for the RS to be refreshed by slower ejecta catching up with the FS (Sari & Mészáros 2000). This scenario requires a distribution of Lorentz factors in the jet, which has indeed been inferred in other events (Laskar et al. 2015, 2018c, 2018a) and remains feasible here. We note that a similar discrepancy was present in the only other ALMA light curve of a GRB available at this date, that of GRB 161219B (Laskar et al. 2018a). In that case, the model underpredicts the data at ≈3 days, corresponding to a similar RS–FS transition. These similarities may indicate that the RS evolution is more complex than previously thought, or that the millimeter-band emission is not dominated by the RS at all, and is in fact a separate component. Further ALMA 3 mm observations of GRBs, together with detailed multiwavelength modeling, may help shed light on this issue.

Zoom In Zoom Out Reset image size

The observed 1 keV light curve exhibits a late-time shallowing (Figure 3) and can therefore not be fit perfectly in our model (Figure 5), where only steepenings are expected (for instance, following a jet break; Sari et al. 1999). As a result, our best-fit model underpredicts the X-ray flux density at 127 days. Above ${\nu }_{{\rm{m}},{\rm{f}}}$, the afterglow light curves are indeed expected to become shallower at late times when the FS becomes nonrelativistic at ${t}_{\mathrm{NR}}$ (Waxman et al. 1998). One simple definition of ${t}_{\mathrm{NR}}$ corresponds to the time when Γβ ≈ 1 (i.e., ${\rm{\Gamma }}\approx \sqrt{2}$). For the best-fit model, these parameters correspond to ${t}_{\mathrm{NR}}\approx 220$ yr, which cannot explain the X-ray excess at ≈127 days. Another possibility is that the late-time excess corresponds to a contaminating source within the XRT point-spread function, which has been observed at least once before; however, this possibility is quite rare and unlikely. One way to test this is to obtain observations of the afterglow with Chandra; however, at the time of writing, the target was within an extended Sun constraint for Chandra and not visible.

We note that late-time excess X-ray flux has been uncovered in other events (Fong et al. 2014; Margutti et al. 2015; Laskar et al. 2018b). Suggested explanations have included dust echoes, inverse-Compton scattering, late-time radiation from the spin down of a magnetar central engine, and interaction of the FS with the shocked stellar wind at the edge of the stellar wind bubble (Shao & Dai 2007; Chandra et al. 2008; Liang et al. 2008; Mimica & Giannios 2011; Gat et al. 2013; Fong et al. 2014; Margutti et al. 2015). For reference, the shock radius at 127 days is ≈30 pc for our best-fit parameters, which is significantly larger than the mean distance between O stars in a typical massive stellar cluster (Figer et al. 1999).

The 1–2 GHz GMRT light curves appear to exhibit an unexpected brightening of the afterglow between 30 and 48 days after the burst. The corresponding radio maps were dynamic range limited, with noise dominated by image deconvolution errors. Thus, it is possible that our reported uncertainties underestimate the true variation of the non-Gaussian image noise. In that scenario, the light curve may be consistent with being flat and therefore host dominated. We caution that while a similar brightening is evident in the GMRT 750 MHz observations, they appear to be absent at 620 MHz. The two bands were observed together as part of a broader frequency coverage with the upgraded GMRT (GWB) system, and no systematic errors are expected in the light curve produced from different parts of the observing band.

On the other hand, we note that late-time light-curve flattenings or rebrightenings have been observed in several radio afterglows (Berger et al. 2003a; Panaitescu & Kumar 2004; Alexander et al. 2017; Laskar et al. 2018a, 2018b). If real, the rebrightening in this event might indicate (i) late-time energy injection due to slower ejecta shells³⁰ (Sari & Mészáros 2000; Kathirgamaraju et al. 2016), which has been inferred in several instances (Laskar et al. 2015, 2018a, 2018c), or (ii) interaction of the FS with the wind termination shock, a scenario discussed above in the context of the X-ray light-curve flattening (Section 5.3). Any contribution from the emerging supernova is expected to peak on ≳10 yr timescales and therefore an unlikely explanation for this excess (Barniol Duran & Giannios 2015). The emergence of a counterjet, such as has been claimed in centimeter-band observations of GRB 980703 at ≳500 days after the burst (Li & Song 2004; Perley et al. 2017) should occur at $\gtrsim {t}_{\mathrm{NR}}\approx 220\,\mathrm{yr}$ for this event (Section 5.3) and is also disfavored.

We note that our derived value of g ≈ 2.7, while higher than the theoretically expected upper limit for a wind environment (Zou et al. 2005), is comparable to g ≈ 2.3 in GRB 140304A (Laskar et al. 2015). Previous studies have found even higher values of g ≈ 5 in the case of the RS emission from GRB 130427A (Laskar et al. 2013; Perley et al. 2014). Higher values of g result in a more slowly evolving spectrum and are preferred by the fit due to the relatively slowly evolving centimeter/millimeter emission observed in these events, including in the event studied here. Owing to this mismatch with the theoretical prediction for g, as well as some of the other issues described above, it is natural to ask whether the RS+FS model is the only one that can explain the observations. It is possible that the emission observed here and in similar events in the past, which has been ascribed to the RS, is in fact from a different emission mechanism.

One possible alternate explanation is to invoke a two-component jet, one in which the optical and X-ray emission arises from a narrower, faster jet than that producing the centimeter-band observations (Berger et al. 2003c; Peng et al. 2005; Oates et al. 2007; Racusin et al. 2008; Holland et al. 2012). In such a case, we expect the radio peak frequency and peak flux density to evolve according to the standard FS prescription as ${\nu }_{{\rm{m}}}\propto {t}^{-3/2}$ and ${F}_{\nu ,{\rm{m}}}\propto {t}^{-1/2}$, respectively. From the SED at 3.9 days (Section 3.1 and Figure 1), this would imply an R-band peak at ≈2.2 × 10⁻³ days with a flux density of ≈60 mJy. While this would overpredict the MASTER observations at ≈2.4 × 10⁻² days, if the jet producing this emission component were still decelerating at this time, then the flux could be significantly lower. Furthermore, we find that this component would not contribute significantly to the remainder of the optical light curve. For instance, at the time of the earliest UVOT observation at ≈0.18 days, the predicted flux density of ≈0.2 mJy is below the observed flux density of ≈0.6 mJy. Whereas it is not possible for us to confirm or rule out two-component jet models, we note that decoupling the centimeter/millimeter emission from the optical/X-ray emission reduces the predictive power of such a model.

Another alternate scenario to explain the radio spectral bump at 3.9 days is the presence of a population of non-accelerated ("thermal") electrons, which are not accelerated by the passage of the FS into a relativistic power-law distribution, but instead form a Maxwellian distribution at lower energies (Eichler & Waxman 2005). However, a disjoint spectral component such as that observed here would likely require a very specific electron distribution, one in which the typical Lorentz factor of the thermal electrons is much lower than the minimum Lorentz factor of the shock-accelerated electrons, the so-called "cold electron model" (Ressler & Laskar 2017).

The emission signature from such a population is expected to be quite broad. We demonstrate this by setting the electron participation fraction, usually assumed to be unity, to f_NT ≈ 0.6, and computing the resultant spectrum as described in Ressler & Laskar (2017). We use a cold electron model, with the ratio of electron temperature (T_e) to gas temperature (T_g) set at ${\eta }_{{\rm{e}}}\equiv {T}_{{\rm{e}}}/{T}_{{\rm{g}}}\approx {10}^{-2}$. The resultant spectrum exhibits a broad peak in the centimeter band and does not match the observations as well as the RS model (Figure 13). While increasing f_NT increases the sharpness of the peak (relative to the underlying FS emission), doing so also increases the relative brightness of the thermal bump beyond the data, and so no combination of f_NT and η_e appears to match the centimeter-band data. A more detailed analysis requires a complete parameter search over f_NT and ${\eta }_{{\rm{e}}}$ that also includes the FS parameters. Such a study is currently not possible due to the computational intensiveness of the thermal electron calculation, and we defer it to future work.

Zoom In Zoom Out Reset image size

Our modeling reveals degeneracies between the FS model parameters (Figure 9), which cannot be explained by the unknown value of ${\nu }_{{\rm{a}},{\rm{f}}}$ alone. For instance, the correlation contours between ${A}_{* }$ and ${E}_{{\rm{K}},\mathrm{iso}}$ lie parallel to the degenerate line ${E}_{{\rm{K}},\mathrm{iso}}\propto {A}_{* }^{1.4}$. However, the expected relation between these parameters when ${\nu }_{{\rm{a}}}$ is unknown is ${E}_{{\rm{K}},\mathrm{iso}}\propto {A}_{* }^{-0.2}$ (Granot & Sari 2002; Laskar et al. 2014). We note that the Compton Y parameter varies by about two orders of magnitude along the axis of the observed degeneracy, and thus the change in the physical parameters is compensated for by varying Y.

Our prescription for calculating Compton Y relies on a simplistic global approximation (Sari & Esin 2001; Laskar et al. 2015). However, we note that in low-density environments, the KN suppression limits the maximum value of Y (Nakar et al. 2009). We derive an expression for this maximum value of Y corresponding to the Lorentz factor of cooling electrons in a wind environment by substituting Equation (18) from Panaitescu & Kumar (2000) for the effective density as a function of observer time in Equation (64) from Nakar et al. (2009) for the maximum value of Y when KN corrections are included to obtain

$$\begin{eqnarray}\begin{array}{rcl}{Y}_{\max }({\gamma }_{{\rm{c}}}) & = & 10{e}^{\tfrac{27{\left(3.2-p\right)}^{2}-17}{(p+2)(4-p)}}{e}^{\left(\tfrac{8-p}{2(p+2)}\mathrm{ln}0.73\right)}\\ & & \times \,{\epsilon }_{{\rm{e}},-1}^{\tfrac{5(p-1)}{p+2}}{E}_{53}^{\tfrac{p-3}{p+2}}{A}_{* }^{\tfrac{8-p}{p+2}}{t}_{\mathrm{days}}^{-\tfrac{(2p-1)}{(p+2)}}.\end{array}\end{eqnarray} \tag{ 9 }$$

For our best-fit parameters (Table 4), this evaluates to ${Y}_{\max }({\gamma }_{{\rm{c}}})\approx 2{t}_{\mathrm{days}}^{-0.8}$, which is lower than some of the values in the degenerate models discussed. Thus, a future, more accurate calculation incorporating KN corrections may resolve this source of degeneracy.

The data do not reveal a jet break (Rhoads 1999; Sari et al. 1999), and hence the degree of collimation remains unknown for this burst. As of the time of writing, we are continuing VLA observations of this event and expect that late-time observations may allow for beaming-independent calorimetry provided the contribution from the underlying host is negligible (Frail et al. 2000; Berger et al. 2003c, 2004; Frail et al. 2005; Shivvers & Berger 2011; De Colle et al. 2012). With that said, we note that our modeling does require some host contribution at 2.7–24.5 GHz at the ≈50 μJy level (Figure 5), corresponding to a star formation rate of ≈2 M_⊙ yr⁻¹ (Yun & Carilli 2002; Berger et al. 2003a); hence, it is not clear whether such calorimetry will be feasible for this burst. For a jet break time of ${t}_{\mathrm{jet}}\gtrsim 163$ days, we can only set a lower limit to the jet half-opening angle of ${\theta }_{\mathrm{jet}}\gtrsim 6^\circ$, corresponding to a minimum beaming-corrected energy of 1.4 × 10⁵² erg.

The radiative efficiency, on the other hand, is independent of the beaming correction. For our best-fit model, we find ${\eta }_{\mathrm{rad}}\approx 0.35$. This value corresponds to a slightly lower efficiency than for other GRBs at similar ${E}_{\gamma ,\mathrm{iso}}$ as derived from X-ray afterglow observations (Lloyd-Ronning & Zhang 2004). However, the observed distribution of ${\eta }_{\mathrm{rad}}$ for Swift GRBs spans a broad range (Berger et al. 2003b; Nousek et al. 2006; Zhang et al. 2007; Beniamini et al. 2015), and this value of ${\eta }_{\mathrm{rad}}$ is unremarkable.

There is a strong observed correlation between multiwavelength detections of RS emission with low-density circumburst environments, typically ${n}_{0}$, ${A}_{* }\lesssim {10}^{-2}$ (Laskar et al. 2018a). We note that the circumburst density for our best-fit model, ${A}_{* }\approx 0.02$, is also similarly low. We have previously remarked on the possibility that this is an observational selection effect, because low-density environments are more likely to result in a slow-cooling RS with long-lasting emission (Laskar et al. 2013, 2016, 2018a, 2019; Kopac et al. 2015; Alexander et al. 2017). In the case of GRB 161219B, we used constraints on the RS cooling frequency to determine a maximal particle density for the ISM environment under which RS emission was expected and showed that the observed density satisfied the constraint in that case. We repeat that argument here, now for the wind medium.

Combining Equation (2) with the expression for FS cooling frequency from Granot & Sari (2002),

$$\begin{eqnarray}\begin{array}{rcl}{\nu }_{{\rm{c}},{\rm{f}}} & = & 4.40\times {10}^{10}(3.45-p){e}^{0.45p}{\left(1+z\right)}^{-3/2}\\ & & \times \,{\epsilon }_{{\rm{B}}}^{-3/2}{A}_{* }^{-2}{E}_{{\rm{K}},\mathrm{iso},52}^{1/2}{\left(1+{Y}_{{\rm{f}}}\right)}^{-2}{t}_{{\rm{d}}}^{1/2}\mathrm{Hz}\\ & \approx & {10}^{13}{A}_{* }^{-2}\mathrm{Hz}\\ & \sim & {R}_{{\rm{B}}}^{3}{\nu }_{{\rm{c}},{\rm{r}}}\end{array}\end{eqnarray} \tag{ 10 }$$

for the best-fit parameters at t ≈ 1 days yields

$$\begin{eqnarray}\begin{array}{rcl}{A}_{* } & \approx & {\left(\displaystyle \frac{{\nu }_{{\rm{c}},{\rm{r}}}}{{10}^{13}\mathrm{Hz}}\right)}^{-1/2}{R}_{{\rm{B}}}^{-3/2}\\ & \lesssim & 0.24{R}_{{\rm{B}}}^{-3/2}{\mathrm{cm}}^{-3},\end{array}\end{eqnarray} \tag{ 11 }$$

which is indeed satisfied for GRB 181201A where ${A}_{* }\,\approx 2\times {10}^{-2}$ (corresponding to a mass-loss rate $\dot{M}\approx 2\,\times {10}^{-7}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ for a wind speed of 10³ km s⁻¹). We note that a similar low wind parameter of ${A}_{* }\approx 3\times {10}^{-3}$ was inferred in the case of GRB 130427A, which also exhibited a strong RS signature (Laskar et al. 2013; Perley et al. 2014). Our observation of RS emission in GRB 181201A therefore supports the hypothesis that RS emission is more easily detectable in low-density environments, in both the ISM and wind scenarios.

During our analysis above, we have assumed a nonrelativistic RS. We now check this assumption by calculating the Lorentz factor of the RS itself. During the shock crossing, the Lorentz factor of the shocked fluid in the cold ejecta frame is given by

$$\begin{eqnarray}&&\bar{\gamma }\approx \displaystyle \frac{1}{2}\left(\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{FS}}/\sqrt{2}}{{{\rm{\Gamma }}}_{0}}+\displaystyle \frac{{{\rm{\Gamma }}}_{0}}{{{\rm{\Gamma }}}_{\mathrm{FS}}/\sqrt{2}}\right),\end{eqnarray} \tag{ 12 }$$

where Γ_FS is the FS Lorentz factor, assumed to be highly relativistic. The RS Lorentz factor is then given by the shock jump condition (Blandford & McKee 1976)

$$\begin{eqnarray}&&{\bar{{\rm{\Gamma }}}}_{\mathrm{RS}}=\displaystyle \frac{(\bar{\gamma }+1){\left[\hat{\varsigma }\left(\bar{\gamma }-1\right)+1\right]}^{2}}{\hat{\varsigma }(2-\hat{\varsigma })(\bar{\gamma }-1)+2},\end{eqnarray} \tag{ 13 }$$

where we take the adiabatic index $\hat{\varsigma }=4/3$, appropriate for a relativistic fluid (Kobayashi & Sari 2000). We compute these quantities from our MCMC results and find ${\bar{{\rm{\Gamma }}}}_{\mathrm{RS}}={1.97}_{0.07}^{+0.08}$, corresponding to a shock speed of ${\bar{\beta }}_{\mathrm{RS}}=0.74\pm 0.02$ and four-velocity ${\bar{{\rm{\Gamma }}}}_{\mathrm{RS}}{\bar{\beta }}_{\mathrm{RS}}=1.7\pm 0.1$. Thus, the RS is at most mildly relativistic and definitely not ultrarelativistic (${\bar{{\rm{\Gamma }}}}_{\mathrm{RS}}{\bar{\beta }}_{\mathrm{RS}}\,\gg 1$). Whereas the expressions we have used for the RS evolution in the modeling are strictly valid only in the asymptotic regime of a Newtonian RS, they are expected to be approximately correct for cold shells where the postshock evolution deviates from the Blandford–McKee solution. A more detailed analysis requires numerical calculation of the RS hydrodynamics and radiation simultaneously with the MCMC fitting, which is beyond the scope of this work.

The relative magnetization of ${R}_{{\rm{B}}}\approx 0.6$ observed in this event is similar to ${R}_{{\rm{B}}}\approx 1$ for GRB 161219B (Laskar et al. 2018a) and ${R}_{{\rm{B}}}\approx 0.6$ for GRB 140304A (Laskar et al. 2018c). On the other hand, these values are all smaller than $1\lesssim {R}_{{\rm{B}}}\lesssim 5$ for GRB 130427A (Laskar et al. 2013), ${R}_{{\rm{B}}}\approx 8$ for GRB 160509A (Laskar et al. 2016), ${R}_{{\rm{B}}}\sim 1\mbox{--}100$ for GRB 160509A (Alexander et al. 2017), and ${R}_{{\rm{B}}}\approx 1$–100 derived for GRBs with optical flashes (Gomboc et al. 2008, 2009; Japelj et al. 2014). Interpreting these values requires accounting for the variation in ${\epsilon }_{{\rm{B}},{\rm{f}}}$ between events, which we perform next.

The plasma magnetization in the ejecta, σ, defined as the ratio of the magnetic energy density to the rest-mass energy density, corresponds to $\sigma \approx {\epsilon }_{{\rm{B}},{\rm{r}}}{\bar{{\rm{\Gamma }}}}_{\mathrm{RS}}$. For the nonrelativistic RSs observed in all of these cases, this implies $\sigma \approx {\epsilon }_{{\rm{B}},{\rm{r}}}={R}_{{\rm{B}}}^{2}{\epsilon }_{{\rm{B}},{\rm{f}}}$. A detailed statistical comparison of σ between GRBs requires full posterior density functions for ${R}_{{\rm{B}}}$ in each case, which are not currently available. However, we can use the typical parameters in each case for a rough comparison. For our best-fit parameters, we have ${\epsilon }_{{\rm{B}},{\rm{r}}}\approx 2.6\times {10}^{-3}$ for GRB 181201A. For GRB 160625B, the large feasible range of ${R}_{{\rm{B}}}\sim 1$–100 makes comparisons less meaningful. For GRB 160509A, ${R}_{{\rm{B}}}^{2}{\epsilon }_{{\rm{B}},{\rm{f}}}\gt 1;$ however, the RS in that case was long lived and continuously refreshed, thus suggesting additional physics not captured by the simple calculation of σ. For the other bursts with radio RS emission, we find ${\epsilon }_{{\rm{B}},{\rm{r}}}\approx 6\times {10}^{-2}$ (161219B), ${\epsilon }_{{\rm{B}},{\rm{r}}}\approx 2\,\times {10}^{-2}$ (140304A), and ${\epsilon }_{{\rm{B}},{\rm{r}}}\approx 0.2-1$ (130427A). Thus, for all GRBs with detected and well-constrained RS emission, σ ranges between 10⁻³ and ≈1, indicating that the ejecta are weakly magnetized. However, we caution that this outcome is largely built into the formalism for nonrelativistic RSs, for which $\sigma \approx {\epsilon }_{{\rm{B}},{\rm{r}}}$. The multiwavelength detection and characterization of a relativistic RS might enable us to explore a broader regime in σ.

In summary, the detection of the RS signature in this GRB argues against the scenario where the jet is Poynting flux dominated, as in that case, no RS emission is expected (Lyutikov & Camilo Jaramillo 2017). The mildly relativistic RS observed in this case argues against propagation in a highly magnetized, relativistic outflow, such as in a magnetar wind (Metzger et al. 2011). This event brings the number of GRBs with radio RS detections to seven (Kulkarni et al. 1999; Sari & Piran 1999a; Laskar et al. 2013, 2016, 2018a, 2018c, 2019; Perley et al. 2014; Alexander et al. 2017).