The late-time afterglow evolution of long gamma-ray bursts GRB 160625B and GRB 160509A

TL;DR

This study analyzes late-time afterglow observations of two long gamma-ray bursts, revealing complex jet behaviors, potential magnetar engines, and challenges to standard jet models through multi-wavelength data and modeling.

Contribution

It provides detailed post-jet-break observations and modeling of GRB afterglows, highlighting deviations from simple models and suggesting complex jet structures and magnetar central engines.

Findings

Post-break decline inconsistent with simple jet edge models.

Radio light curves suggest highly chromatic jet breaks.

Magnetar engine models require extreme parameters.

Abstract

We present post-jet-break \textit{HST}, VLA and \textit{Chandra} observations of the afterglow of the long $\gamma$-ray bursts GRB 160625B (between 69 and 209 days) and GRB 160509A (between 35 and 80 days). We calculate the post-jet-break decline rates of the light curves, and find the afterglow of GRB 160625B inconsistent with a simple $t^{-3/4}$ steepening over the break, expected from the geometric effect of the jet edge entering our line of sight. However, the favored optical post-break decline ($f_{\nu} \propto t^{-1.96 \pm 0.07}$) is also inconsistent with the $f_{\nu} \propto t^{-p}$ decline (where $p \approx 2.3$ from the pre-break light curve), which is expected from exponential lateral expansion of the jet; perhaps suggesting lateral expansion that only affects a fraction of the jet. The post-break decline of GRB 160509A is consistent with both the $t^{-3/4}$ steepening and with $f_{\nu} \propto t^{-p}$. We also use {\sc boxfit} to fit afterglow models to both light curves and find both to be energetically consistent with a millisecond magnetar central engine, although the magnetar parameters need to be extreme (i.e. $E \sim 3 \times 10^{52}$ erg). Finally, the late-time radio light curves of both afterglows are not reproduced well by {\sc boxfit} and are inconsistent with predictions from the standard jet model; instead both are well represented by a single power law decline (roughly $f_{\nu} \propto t^{-1}$) with no breaks. This requires a highly chromatic jet break ($t_{j,\mathrm{radio}} > 10 \times t_{j,\mathrm{optical}}$) and possibly a two-component jet for both bursts.

Figures

Tables

Table 1. Table 1: Log of our late-time HST /WFC3 observations of GRB 160625B.

Phase	MJD	$t_{\exp }$	F606W	corrected F606W
(d)		(s)	(mag)	(mag)
71.5	57636.4	2400	$25.38\pm 0.03$	$25.36\pm 0.04$
140.2	57705.1	4800	$26.76\pm 0.06$	$26.67\pm 0.07$
498.3	58063.2	4800	…	…
503.6	58068.5	4800	…	…

Table 2. Table 2: Log of our late-time HST /WFC3 and GTC/CIRCE observations of GRB 160509A.

Phase	MJD	$t_{\exp ,\mathrm{F110W}}$	F110W	$t_{\exp ,\mathrm{F160W}}$	F160W	$t_{\exp ,H}$	$H$
(d)		(s)	(mag)	(s)	(mag)	(s)	(mag)
5.8	57523.2	…	…	…	…	3060	$20.50\pm 0.17$
24.8	57542.2	…	…	…	…	2100	$\ge 21.9$
35.3	57552.7	2697	$27.11\pm 0.10$	2797	$26.07\pm 0.07$	…	…
422.1	57939.5	2697	…	2797	…	…	…

Table 3. Table 3: Log of our late-time Chandra /ACIS-S observations of GRB 160625B and GRB 160509A.

Phase	MJD	$t_{\exp }$	$f_{\nu }$(5 keV)
(d)		(ks)	(erg s-1 cm-2 keV-1)
160625B
69.8	57634.7	19.80	($1.47\pm 0.29)\times 10{}^{-15}$
142.3	57707.3	45.84	…
$144.3\pm 2.2$777Combination of the 142.3 and 146.2 d epochs.	$57709.3\pm 2.2$	69.56	($3.21\pm 0.79)\times 10{}^{-16}$
146.2	57711.2	23.72	…
160509A
42.1	57559.5	24.75	($1.38\pm 0.25)\times 10{}^{-15}$

Table 4. Table 4: Log of our VLA radio observations of GRB 160625B and GRB 160509A. The GRB 160625B points until 31.3 d were also reported in Troja et al. ( 2017 ) , but without the calibration uncertainty.

Phase	MJD	$\nu$	$f_{\nu }$	Configuration
(d)		(GHz)	($\mu$Jy)
160625B
4.5	57569.4	4.8	$104\pm 16$	B
4.5	57569.4	7.4	$454\pm 27$	B
4.5	57569.4	19	$278\pm 35$	B
4.5	57569.4	25	$204\pm 36$	B
13.4	57578.3	4.8	$377\pm 25$	B
13.4	57578.3	7.4	$310\pm 21$	B
13.4	57578.3	22	$163\pm 20$	B
31.3	57596.2	7.4	$113\pm 16$	B
31.3	57596.2	22	$88\pm 19$	B
58.3	57623.2	6.1	$75\pm 11$	B
58.3	57623.2	22	$52\pm 13$	B
209.0	57773.9	6.1	$16\pm 5$	A
160509A
23.9	57541.3	6.0	$80\pm 8$	B
23.9	57541.3	9.0	$71\pm 7$	B
36.9	57554.3	5.0	$50\pm 7$	B
36.9	57554.3	6.9	$52\pm 7$	B
36.9	57554.3	8.5	$41\pm 6$	B
36.9	57554.3	9.5	$29\pm 6$	B
79.9	57597.3	6.0	$27\pm 6$	B
79.9	57597.3	9.0	$25\pm 5$	B

Table 5. Table 5: Parameters of the best smooth broken power law fits to the GRB 160625B light curves and of the single power law (SPL) fits to the early and late decline, ignoring the bump(s) between 8.5 and 26.5 d.

Parameter	$\omega =3$	$\omega =10$	SPL
$t_{j,\mathrm{F606W}}$	$24\pm 3$ d	$17\pm 2$ d	$17\pm 4$ d
${\alpha }_{1,\mathrm{F606W}}$	$-0.96\pm 0.01$	$-0.96\pm 0.01$	$-0.97\pm 0.01$
${\alpha }_{2,\mathrm{F606W}}$	$-2.27\pm 0.13$	$-1.96\pm 0.07$	$-1.94\pm 0.13$
Reduced ${\chi }^2$	5.5	4.4	1.8
$t_{j,X}$	$27\pm 5$ d	$22\pm 4$ d	$22\pm 5$ d
${\alpha }_{1,X}$	$-1.24\pm 0.02$	$-1.24\pm 0.02$	$-1.25\pm 0.03$
${\alpha }_{2,X}$	$-2.40\pm 0.19$	$-2.23\pm 0.15$	$-2.20\pm 0.13$
Reduced ${\chi }^2$	0.91	0.81	0.84

Table 6. Table 6: Best-fit physical parameters of the best boxfit fits to GRB 160625B at three different values of the participation fraction ξ 𝜉 \xi .

Parameter	$\xi =1$	$\xi =0.1$	$\xi =0.01$
$p$	2.30	2.05	2.05
$E_{K,\mathrm{iso}}$ (erg)	$1.8\times {10}^{54}$	$1.4\times {10}^{54}$	$1.3\times {10}^{55}$
${ϵ}_e$	0.13	0.25	0.024
${ϵ}_B$	$0.030$	$3.0\times {10}^{-4}$	$5.8\times {10}^{-5}$
$n$ (cm-3)	$1.1\times {10}^{-5}$	0.18	0.96
${\theta }_j$ (rad)	0.059	0.14	0.13
${\theta }_j$ (deg)	$3.4$	7.8	7.2
${\theta }_{obs}$ (rad)	$0.012$	$1.1\times {10}^{-3}$	$1.1\times {10}^{-3}$
${\theta }_{obs}$ (deg)	$0.69$	0.07	0.06
$E_{\mathrm{tot}}$ (erg)	$8.3\times {10}^{51}$	$4.1\times {10}^{52}$	$1.3\times {10}^{53}$
$\eta$	0.62	0.68	0.19
${\chi }^2/d.o.f$	8.6	4.6	4.5

Table 7. Table 7: Parameters of the best smooth broken power law fits to the GRB 160509A X-ray light curve.

Parameter	$\omega =3$	$\omega =10$
$t_{j,X}$	$3.2\pm 0.9$ d	$3.7\pm 0.8$ d
${\alpha }_{1,X}$	$-1.06\pm 0.10$	$-1.20\pm 0.06$
${\alpha }_{2,X}$	$-1.98\pm 0.10$	$-1.96\pm 0.09$
Reduced ${\chi }^2$	0.84	0.85

Table 8. Table 8: Best-fit physical parameters of the best boxfit fits to GRB 160509A at three different values of the participation fraction ξ 𝜉 \xi .

Parameter	$\xi =1$	$\xi =0.1$	$\xi =0.01$
$p$	2.29	$2.13$	2.05
$E_{K,\mathrm{iso}}$ (erg)	$8.5\times {10}^{53}$	$3.8\times {10}^{53}$	$3.8\times {10}^{55}$
${ϵ}_e$	0.19	0.45	$5.7\times {10}^{-3}$
${ϵ}_B$	0.015	$1.7\times {10}^{-5}$	$5.8\times {10}^{-4}$
$n$ (cm-3)	$2.1\times {10}^{-5}$	$18.1$	$6.1\times {10}^{-3}$
${\theta }_j$ (rad)	0.046	0.20	0.045
${\theta }_j$ (deg)	2.6	11.5	2.6
${\theta }_{obs}$ (rad)	0.026	0.12	0.027
${\theta }_{obs}$ (deg)	1.5	7.0	1.5
$E_{\mathrm{tot}}$ (erg)	$1.7\times {10}^{51}$	$2.5\times {10}^{52}$	$3.9\times {10}^{52}$
$\eta$	0.50	0.69	0.02
${\chi }^2/d.o.f$	1.8	1.9	1.8

Equations

$$f_{\nu}=f_{\nu,0}\Big{[}\Big{(}\frac{t}{t_{j}}\Big{)}^{-\omega\alpha_{1}}+\Big{(}\frac{t}{t_{j}}\Big{)}^{-\omega\alpha_{2}}\Big{]}^{-\frac{1}{\omega}}$$

$$E_{\mathrm{rot}}\approx 2\times 10^{52}~{}\mathrm{erg}~{}\frac{M}{1.4~{}M_{\odot}}\Big{(}\frac{R}{10~{}\mathrm{km}}\Big{)}^{2}\Big{(}\frac{P_{0}}{1~{}\mathrm{ms}}\Big{)}^{-2}~{},$$

Full text

gamma-ray burst: general — gamma-ray burst: individual (GRB 160625B; GRB 160509A) — relativistic processes

††facilities: HST(WFC3), Chandra(ACIS-S), VLA, GTC(CIRCE)††software: drizzlepac (Hack et al., 2013); pyraf (Science Software Branch at STScI, 2012); iraf (Tody, 1986); ciao (Fruscione et al., 2006); casa (McMullin et al., 2007); isis (Alard, 2000; Alard & Lupton, 1998); boxfit (van Eerten et al., 2012)

1 Introduction

Gamma-ray bursts (GRBs) are among the most luminous transient events in the universe. Through their association with broad-lined type Ic supernovae (e.g. Iwamoto et al., 1998; Woosley & Bloom, 2006; Hjorth & Bloom, 2012), long GRBs (LGRBs; duration of the prompt $\gamma$-ray emission more than 2 s) have been established as the terminal core-collapse explosions of massive stars at cosmological distances (e.g. Paczynski, 1986; Woosley, 1993; MacFadyen & Woosley, 1999), where an ultra-relativistic jet is launched and breaks out of the stellar envelope, generating the initial prompt emission of $\gamma$ rays through an as yet unclear mechanism (for a review on GRB physics, see e.g. Piran, 2004; Kumar & Zhang, 2015). The central engine responsible for launching the jet and powering the emission may be either accretion onto a black hole formed in the core collapse (Woosley, 1993) or rotational energy released through the spin-down of a nascent magnetar (e.g. Bucciantini et al., 2008, 2009). The prompt emission of a GRB is followed by an afterglow from X-ray to radio frequencies – synchrotron emission from an external shock created by the interaction between the circumburst medium (CBM) and the highly collimated and relativistically beamed jet (e.g. Paczynski & Rhoads, 1993; Sari et al., 1998; Piran, 2004). The flux density of the afterglow declines as a power law of the form $f_{\nu}\propto t^{\alpha}$.

As the jet interacts with the CBM, it decelerates and the relativistic beaming effect diminishes over time (on the order of days or weeks after a long GRB; e.g. Racusin et al., 2009). This results in an achromatic jet break in the afterglow light curve when the relativistic beaming angle ($\Gamma^{-1}$, where $\Gamma$ is the bulk Lorentz factor in the jet) becomes comparable to the opening angle of the jet (Rhoads, 1999; Sari et al., 1999), with a steeper power-law decline after the break. The post-break decline is affected by a geometric ’edge effect’, in contrast to the situation pre-break where the observer only sees a fraction of the jet front and hence behaviour consistent with an isotropic fireball model. This phenomenon is believed to steepen the decline slope $\alpha$ by $-3/4$ over the break assuming a constant-density CBM, or by $-1/2$ in the case of a wind-like CBM (e.g. Mészáros & Rees, 1999; Panaitescu & Mészáros, 1999; Kumar & Zhang, 2015). Another effect is that, around the same time as this happens, transverse sound waves become able to cross the jet and lateral expansion starts, exponentially decelerating the shock wave. Theoretically the post-break slope in this scenario is expected to be equal to $-p$ (e.g. Sari et al., 1999), where $p$ is the index of the electron Lorentz factor distribution ($N(\gamma)\propto\gamma^{-p}$), typically estimated to be between 2 and 3. There is, however, evidence from numerical simulations that the lateral expansion is unimportant until a later stage – at least unless the jet is very narrow, $\theta_{j}\lesssim 3$ deg (Lyutikov, 2012; Granot & Piran, 2012). At even later times, the jet is expected to be better described as a non-relativistic fireball in the Sedov-von Neumann-Taylor regime, resulting in a somewhat flatter decline (e.g. Frail et al., 2000; van der Horst et al., 2008).

Simulations of relativistic shocks have resulted in values around $p\approx 2.2$ (e.g. Bednarz & Ostrowski, 1998; Gallant et al., 1999; Kirk et al., 2000). In the X-rays, the pre-break light curve tends to follow a decline around $t^{-1.2}$ (albeit with some variation; e.g. Piran, 2004; Zhang et al., 2006); thus both of these effects result in a roughly similar post-break decline (i.e. $\sim t^{-2}$, though with high uncertainties due to the the fast decline and the resulting faintness; often there are not enough data to distinguish between $t^{-1.9}$ and $t^{-2.2}$). Thus determining the exact scenario observationally requires late-time observations of the rapidly declining afterglows to constrain this slope.

The Large Area Telescope (LAT) on the Fermi Gamma-ray Space Telescope has detected a number of GRBs at relatively high energies (MeV to GeV) since the launch of Fermi in 2008. These are often among the most energetic GRBs, consistent with the Amati correlation between isotropic-equivalent energy $E_{\mathrm{iso}}$ and the peak of the energy spectrum (Amati et al., 2002), and can haveisotropic-equivalent energies on the order of $10^{54}$ erg (Cenko et al., 2011). Some of these most energetic bursts do not exhibit the expected jet breaks, suggesting larger opening angles than expected and making them even more energetic intrinsically (De Pasquale et al., 2016; Gompertz & Fruchter, 2017). With beaming-corrected energies on the order of $10^{52}$ erg, magnetar spin-down models struggle to produce the required power (Cenko et al., 2011). Thus examining the late-time evolution of the LAT bursts can shed light on the physics of the most energetic GRBs.

In this paper, we present results from our late-time Hubble Space Telescope (HST), Karl G. Jansky Very Large Array (VLA) and Chandra X-ray Observatory imaging observations of the afterglows of two LAT bursts, GRB 160625B and GRB 160509A. GRB 160625B was discovered by the Gamma-ray Burst Monitor (GBM) on Fermi on 2016 June 25 at 22:40:16.28 UT (MJD 57564.9; Dirirsa et al., 2016) and detected by the LAT as well. Xu et al. (2016) determined its redshift to be $z=1.406$. It was one of the most energetic $\gamma$-ray bursts ever observed with $E_{iso}\sim 3\times 10^{54}$ erg (Wang et al., 2017; Zhang et al., 2018), and a well-studied object with a multi-frequency follow-up that revealed signs of a reverse shock within the jet (Alexander et al., 2017). The jet break time was unusually long, around 20 days, as expected from unusually bright GRBs (the median time is $\sim 1$ d, with more energetic bursts having longer break times; see Racusin et al., 2009). GRB 160509A was detected by GBM and LAT on 2016 May 9 at 08:59:04.36 UT (MJD 57517.4; Roberts et al., 2016; Longo et al., 2016a, b) at a redshift of $z=1.17$ (Tanvir et al., 2016). With $E_{iso}=8.6\pm 1.1\times 10^{53}$ erg, this was another luminous burst that exhibited signs of a reverse shock as well (Laskar et al., 2016).

Our observations of GRB 160625B make its follow-up one of the longest post-jet-break optical and X-ray follow-ups of a GRB afterglow111The post-break light curve of GRB 060729 (Grupe et al., 2010) and GRB 170817A (Hajela et al., 2019) has been followed up longer, while GRB 130427A was followed for $\sim 1000$ days (De Pasquale et al., 2016), but exhibited no jet break., thus providing one of the best estimates of the post-break decline in these bands so far, while for GRB 160509A no prior estimates of the infrared/optical post-break decline could be made due to the very sparse light curve.

Our observations and data reduction process are described in Section 2. Our analysis and results are presented in Section 3. In Section 4, we discuss the implications of our findings, and finally present our conclusions in Section 5. All magnitudes are in the AB magnitude system (Oke & Gunn, 1983) and all error bars correspond to $1\sigma$ confidence intervals. We use the cosmological parameters $H_{0}=69.6$ km s*-1* Mpc*-1*, $\Omega_{m}=0.286$ and $\Omega_{\Lambda}=0.714$ (Bennett et al., 2014).

2 Observations and data reduction

Late-time imaging observations of GRB 160625B were performed using HST/WFC3 and the F606W filter on 2016 September 5 (71.5 d) and 2016 November 13 (140.2 d). A template image of the host galaxy was created by combining images obtained with the same setup on 2017 November 6 (498.3 d) and 11 (503.6 d). At this time the contribution of the afterglow itself was a factor of $\sim 13$ fainter than at 140 d, assuming a $f_{\nu}\propto t^{\alpha}$ decline where $\alpha=-2$. Imaging of GRB 160509A in the $H$ band was performed using the Canarias InfraRed Camera Experiment (CIRCE; Eikenberry et al., 2018) instrument on Gran Telescopio Canarias (GTC) on 2016 May 15 (5.8 d) and 2016 June 3 (24.8 d). Late-time imaging of GRB 160509A was done using HST/WFC3 and the F110W and F160W filters on 2016 June 13 (35.3 d); template images of the host galaxy in these filters were obtained on 2017 July 5 (422.1 d), when, assuming $\alpha=-2$, the afterglow was a factor of 143 fainter. Our HST observations of both bursts were executed as part of program GO 14353 (PI Fruchter), and these data are available at10.17909/t9-yvpg-xb33 (catalog 10.17909/t9-yvpg-xb33) (GRB 160625B) and10.17909/t9-11cx-cv41 (catalog 10.17909/t9-11cx-cv41) (GRB 160509A).

Basic reduction and flux calibration of the HST images was performed by the HST calwf3 pipeline. The calibrated images were corrected for distortion, drizzled (Fruchter & Hook, 2002) and aligned to a common world coordinate system using the astrodrizzle, tweakreg and tweakback tasks in the drizzlepac222http://drizzlepac.stsci.edu/ package in pyraf333http://www.stsci.edu/institute/software_hardware/pyraf. The two epochs of GRB 160625B in November 2017 were combined into one template image. Subtraction of the template images and aperture photometry of the afterglows were done using iraf444iraf is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation.. Basic reduction of the GTC/CIRCE data was done using standard iraf tasks. The HST F160W template image was subtracted from the CIRCE images using the isis 2.2 package (Alard & Lupton, 1998; Alard, 2000). Flux calibration was done using field stars in the Two-Micron All Sky Survey (2MASS) catalog555http://www.ipac.caltech.edu/2mass/ (Skrutskie et al., 2006), and aperture photometry was performed using standard iraf tasks. At 24.8 d, we were unable to detect the afterglow and only obtained a ($3\sigma$) limit of $H\geq 21.9$ mag.

The measured magnitudes of GRB 160625B were corrected for over-subtraction caused by the continued presence of a faint afterglow in the template image. Assuming a post-jet-break decline of $\alpha=-2.0\pm 0.2$ (obtained from a single-power-law fit to uncorrected $>25$ d data, with errors rounded up to be conservative), the afterglow flux present in the template image was estimated to be $2.0\pm 1.0$ per cent of the flux at 71.5 d or $7.5\pm 2.6$ per cent of the flux at 140.2 d, and thus the images at these epochs were over-subtracted by approximately these amounts. The magnitudes were adjusted for this; the errors of the corrected magnitudes include an estimate of the uncertainty of the over-subtraction. The magnitudes of GRB 160509A were not corrected, as the contribution of the afterglow in the template image was only estimated to be 0.7 per cent of the 35.3 d brightness. The log of optical observations and measured and corrected magnitudes of GRB 160625B are presented in Table 1, while Table 2 contains the near-infrared observations of GRB 160509A. Figure 1 shows our F606W band images and the resulting template subtractions of GRB 160625B, while Figure 2 shows the F160W image and subtraction of GRB 160509A.

Late-time X-ray imaging of both GRBs was performed using Chandra/ACIS-S in VFAINT mode (proposal ID 17500753, PI Fruchter). GRB 160625B was observed on 2016 September 3 (69.8 d), 2016 November 15 (142.3 d) and 2016 November 19 (146.2 d). The latter two epochs were combined to obtain the flux at $144.3\pm 2.2$ d, as the flux of the afterglow was not expected to vary significantly over a few days at this time. GRB 160509A was observed on 2016 June 20 (42.1 d). Reprocessing of the Chandra level 1 data was performed using the chandra_repro script within the ciao v. 4.9 software (caldb v. 4.7.7; Fruscione et al., 2006), and aperture photometry was done using iraf. The web-based Portable Interactive Multi-Mission Simulator (pimms666https://heasarc.gsfc.nasa.gov/docs/software/tools/pimms.html) was used to convert count rates in the 0.3 – 10 keV range to unabsorbed flux densities at 5 keV. For GRB 160625B, we used a Galactic neutral hydrogen column density $N_{\mathrm{H,MW}}=9.76\times 10^{20}$ cm*-2* (Willingale et al., 2013), a photon index of $\Gamma_{X}=1.86$ and an intrinsic absorption of $N_{\mathrm{H,int}}=2.1\times 10^{21}$ cm*-2* as derived by Alexander et al. (2017). These parameters are also consistent with the initial analysis by Melandri et al. (2016). For GRB 160509A, we used a Galactic neutral hydrogen column density $N_{\mathrm{H,MW}}=2.12\times 10^{21}$ cm*-2* (Willingale et al., 2013), a photon index of $\Gamma_{X}=2.07$ and an intrinsic absorption of $N_{\mathrm{H,int}}=1.52\times 10^{22}$ cm*-2*, following Laskar et al. (2016). $\Gamma_{X}$ is assumed to be constant over the light curve break. The log of X-ray observations and derived flux densities is presented in Table 3.

GRB 160625B was observed in the radio using the VLA in the $C$, $K$, $X$ and/or $Ku$ bands at five epochs between 2016 March 30 (4.5 d) and 2017 January 20 (209.0 d), and GRB 160509A in the $C$ and $X$ bands on 2016 June 2 (23.9 d), 2016 June 15 (36.9 d) and 2016 July 28 (79.9 d) (program IDs S81171 and SH0753, PI Cenko and Fruchter respectively). The observations were done in the B configuration, apart from the last GRB 160625B point where configuration A was used. The log of our observations is presented in Table 4. The data were reduced using the Common Astronomy Software Applications package (CASA; McMullin et al., 2007)888https://casa.nrao.edu. Calibration was carried out using the standard VLA calibration pipeline provided in CASA. For GRB 160625B we used J2049+1003 as our complex gain calibrator and 3C48 as our flux and bandpass calibrator. For GRB 160509A we used J2005+7752 as our complex gain calibrator and 3C48 as our flux and bandpass calibrator. After calibration, the data were manually inspected for radio-frequency interference flagging. Imaging was carried out using the clean algorithm in interactive mode in CASA. Flux densities reported in Table 4 correspond to peak flux densities measured in a circular region centered on the GRB position, with radius comparable to the nominal full width half maximum of the VLA synthesized beam in the appropriate configuration and frequency band. The reported errors include the VLA calibration uncertainty, which is assumed to be 5 per cent below 18 GHz and 10 per cent above it999(https://science.nrao.edu/facilities/vla/docs/manuals/oss/performance/fdscale).

3 Analysis

3.1 GRB 160625B

As our HST observations took place after the jet break, we combined our data set with earlier ground-based observations. Both Alexander et al. (2017) and Troja et al. (2017) have published SDSS $r^{\prime}$ band light curves of GRB 160625B. However, there is a slight ($\sim 0.1$ mag) systematic offset between these data, so in our light curve fits we have only used the Troja et al. (2017) data set, which has a larger number of data points and which was directly tied to the PanSTARRS magnitude system. Magnitudes of GRB 160625B in the $r^{\prime}$ band were converted to flux density at the central wavelength of the F606W filter (5947 Å) assuming a spectral slope of $f_{\nu}\propto\nu^{-0.68}$ between the characteristic synchrotron frequency $\nu_{m}$ and the cooling frequency $\nu_{c}$ (Alexander et al., 2017). As the optical spectrum with $\beta=-0.68\pm 0.07$ is consistent with the expected index of $\beta=-0.65$ when $p=2.3$ (also consistent with the light curve; see Section 4.2.1), host extinction is assumed to be negligible. Optical fluxes have been corrected for Galactic reddening, $E(B-V)=0.1107$ mag (Schlafly & Finkbeiner, 2011), assuming the Cardelli et al. (1989) extinction law. In the X-ray, we combined our Chandra data with the GRB 160625B light curve from the Swift/XRT lightcurve repository101010http://www.swift.ac.uk/xrt_curves/ (Evans et al., 2007, 2009), converted to 5 keV flux densities using pimms as described in Section 2.

We then fitted a smooth broken power law of the form

[TABLE]

to the light curve, where $t_{j}$ is the jet break time, $\alpha_{1}$ is the pre-break power-law slope, $\alpha_{2}$ the post-break slope, and $\omega$ describes the sharpness of the break. We fitted this function to both the optical and the X-ray curve using two values, 3 and 10, for $\omega$ (a value of 3 was found consistent with most GRB observations by Liang et al., 2007, but some events were found to require a sharper break with $\omega=10$). The results of the fit parameters are presented in Table 5. The pre-break decline $\alpha_{1}$ does not depend on the choice of $\omega$; we found $\alpha_{1,\mathrm{F606W}}=-0.96\pm 0.01$ and $\alpha_{1,\mathrm{X}}=-1.24\pm 0.02$ in both cases. The best fit to the post-break decline was $\alpha_{2,\mathrm{F606W}}=-2.27\pm 0.13$ and $\alpha_{2,\mathrm{X}}=-2.40\pm 0.19$ assuming $\omega=3$, and $\alpha_{2,\mathrm{F606W}}=-1.96\pm 0.07$ and $\alpha_{2,\mathrm{X}}=-2.23\pm 0.15$ when $\omega=10$. The optical and X-ray light curves and our best fits in both cases are shown in Figure 3.

We also fitted the decline using a single power law before 8.5 d and another after 26.5 d, ignoring the points in the vicinity of the break itself. The $r^{\prime}$ band light curve contains at least one smooth ’bump’ feature, possibly two depending on $t_{j}$ (we discuss the nature of the bump in Section 4.1). These may disturb the optical broken power-law fits; the reduced $\chi^{2}$ values of these fits are rather high, although the small errors also contribute to this. The result is $\alpha_{2,\mathrm{F606W}}=-1.94\pm 0.13$, nearly exactly coinciding with the $\omega=10$ case but with a $\sim 2.5\sigma$ difference to $\omega=3$. Repeating this in the X-ray results in $\alpha_{2,\mathrm{X}}=-2.20\pm 0.13$, which is also almost identical to the $\omega=10$ case. A simultaneous single power law fit to both post-break light curves results in $\alpha_{2}=-2.01\pm 0.09$.

Assuming an achromatic break, we determined $t_{j}$ by taking the weighted average of $t_{j,\mathrm{F606W}}$ and $t_{j,\mathrm{X}}$. In the $\omega=10$ case, the result is $t_{j}=19\pm 2$ d. Assuming an instantaneous break (corresponding to $\omega=\infty$) between the single power law fits, the resulting jet break times are consistent, $t_{j,\mathrm{F606W}}=17\pm 4$ d and $t_{j,\mathrm{X}}=22\pm 5$ d, and the weighted average $t_{j}=19\pm 3$ d. In the $\omega=3$ case, we obtained $t_{j}=25\pm 3$ d.

For the radio light-curve of GRB 160625B, we combined flux measurements from Alexander et al. (2017) and Troja et al. (2017) with our own data. At 58.3 d and 209.0 d we have observations at 6.1 GHz; we therefore obtained flux densities at 6.1 GHz by power-law interpolation between 5 and 7.1 GHz literature values at 22.5 and 48.4 d. We also scaled the 7.4 GHz flux at 31.34 d assuming the same power law as at 22.5 d. Points earlier than 22.5 d were ignored in the analysis of the late afterglow due to the influence of the reverse shock (Alexander et al., 2017). The resulting best fit for the late-time light curve is $\alpha_{\mathrm{6.1GHz}}=-1.08\pm 0.11$ as shown in Figure 3.

Additionally, we used the boxfit v.2 afterglow fitting code (van Eerten et al., 2012), based on the Afterglow Library111111http://cosmo.nyu.edu/afterglowlibrary/index.html, to fit the light curve. The library of models itself was constructed using the relativistic hydrodynamics code ram (Zhang & MacFadyen, 2006). boxfit then uses a downhill simplex method with simulated annealing to find the best fit, interpolating between these models. We omitted the pre-break radio points due to the influence of the reverse shock in the early light curve, and all the radio points below 5 GHz due to possible strong Milky Way scintillation (Alexander et al., 2017). We also included the ultraviolet to near-infrared frequency data from Troja et al. (2017). We assumed an ISM-like CBM (the light curve rules out a wind-type CBM; see Section 4.2.1) and performed the fit with three different values of the participation fraction $\xi$, i.e. the fraction of electrons accelerated by the shock into a non-thermal power-law distribution. Simulations indicate this value can be as low as 0.01 (Sironi & Spitkovsky, 2011; Sironi et al., 2013; Warren et al., 2018); we used fixed values of 1 (commonly assumed in the literature), 0.1 and 0.01. All other model parameters were allowed to vary within the full range allowed by boxfit. The resulting best-fit parameters are summarized in Table 6. Taking the isotropic-equivalent $\gamma$-ray energy $E_{\mathrm{iso}}=3.0\times 10^{54}$ erg (with the fluence from Svinkin et al., 2016), we also calculate the geometry-corrected total energy and the efficiency $\eta=E_{\mathrm{iso}}/(E_{\mathrm{K,iso}}+E_{\mathrm{iso}})$ for the conversion of kinetic energy to $\gamma$-rays. These fits, however, fail to reproduce the measured power law slope of $\alpha_{\mathrm{6.1GHz}}=-1.08\pm 0.11$, instead predicting a break in the radio light curve around $\sim 100$ d (associated with the passage of $\nu_{m}$ through this band). See Figure 4 for our best boxfit light curve fits. For clarity, we plot the $U$, F606W and $H$ bands, covering the optical/infrared behavior from early to late times, but omit the other optical/infrared bands, which exhibit very similar behavior (see Troja et al., 2017). While the late-time 6.1 GHz light curve can be reproduced slightly better at low $\xi$ values, the fit at higher frequencies or earlier times is still somewhat worse; we show 22 GHz as an example.

As some optical and X-ray observations are nearly contemporaneous, we can construct the spectral energy distribution (SED) of GRB 160625B. Figure 5 shows the SED at four epochs around or after the break, along with spectra produced by boxfit at these epochs. The power-law slope of the SED, $\beta$, between the optical ($r$) and X-ray (5 keV) bands, steepens slightly over time, from $-0.79\pm 0.02$ between 3 and 10 d to $-0.86\pm 0.04$ at 141 d. This is steeper than $-0.65$, expected from $p\approx 2.3$ implied by the early optical and X-ray light curves (see Section 4.2.1) for $\nu<\nu_{c}$, but shallower than $-1.15$, which is expected for $\nu>\nu_{c}$. Alexander et al. (2017) obtain an early X-ray spectral slope similar to this, $\beta_{X}=-0.86^{+0.09}_{-0.10}$, and explain this as $\nu_{c}$ being located just below the X-ray band. However, according to the UKSSDC Swift Burst Analyser121212http://www.swift.ac.uk/burst_analyser/00020667/ the X-ray photon index $\Gamma_{X}$ (and thus the spectral slope in X-ray) does not significantly evolve over the first 30 d but stays around $\sim 1.8$, after which the spectrum seems to flatten to $\Gamma_{X}\sim 1.1$. This feature may not be real, though, as the Burst Analyzer light curve deviates much more from a clean power law when this is used in flux calculation – thus we assume a constant $\Gamma_{X}$131313The post-break X-ray slope would not change by changing $\Gamma_{X}$ at the latest Swift points, as Chandra points would be affected equally – but $t_{j,\mathrm{X}}$ could be delayed.. If $\nu_{c}$ was initially just below X-ray and changed as $\nu_{c}\propto t^{-1/2}$, one would expect the spectrum to instead steepen over time to its $\nu\gg\nu_{c}$ value. We discuss this evolution further in Section 4.2.1.

3.2 GRB 160509A

It was noted in Laskar et al. (2016) that the host galaxy of GRB 160509A contributes substantially to the optical and infrared photometry, and that the event occurred behind a significant amount of extinction in the host galaxy. In order to estimate the host galaxy extinction along the line of sight to the GRB, we removed the foreground Galactic reddening of $E(B-V)=0.2519$ mag (Schlafly & Finkbeiner, 2011) using the Cardelli et al. (1989) law, and assumed a $f_{\nu}\propto\nu^{\beta}$ SED, where $\beta=-0.6$ (consistent with $\nu<\nu_{c}$ and $p\approx 2.2$, determined based on the X-ray spectrum and light curve by Laskar et al., 2016). For the host, we assume the Pei (1992) extinction law for the Small Magellanic Cloud (SMC), as both Kann et al. (2006) and Schady et al. (2012) found the extinction curve in the SMC consistent with their samples. We fitted the observed optical-infrared SED simultaneously at two epochs, corrected using this extinction curve, to find the required extinction correction to match $\beta=-0.6$. The GRB flux in the $g^{\prime}$ band at 1 d was estimated by subtracting the observed flux at 28 d ($g^{\prime}=25.39\pm 0.12$; Laskar et al., 2016) from the flux at 1.0 d ($g^{\prime}=25.03\pm 0.15$; Cenko et al., 2016). The host is assumed to dominate at 28 d due to the flatness of the light curve even after the X-ray break. In the $J$ band, we subtracted the flux of the host galaxy measured in the HST F110W filter (using a 1 arcsec aperture) from the flux at 1.2 d ($J\approx 19.7$; Tanvir et al., 2016). The $r^{\prime}$ band was not included in the SED, as the late and early fluxes are consistent within $1\sigma$ (Cenko et al., 2016; Laskar et al., 2016). Our F110W and F160W observations at 35.3 d made up the other epoch to be fitted simultaneously. The resulting host extinction is $A_{V}=2.8\pm 0.1$ mag in the rest-frame (this is somewhat lower than the result obtained by Laskar et al., 2016, using an afterglow model fit where the host flux was a free parameter). Using the Pei (1992) law, the extinction correction in F160W (approximately $i$-band in the rest frame) is thus 1.5 mag. In the Milky Way, the adopted $N_{\mathrm{H,int}}=1.52\times 10^{22}$ cm-2 would correspond to $A_{V}\approx 6.9$ mag (Güver & Özel, 2009), suggesting a low $A_{V}/N_{H}$ ratio for Milky Way standards but higher than that of most GRB hosts. This ratio is consistent with the $A_{V}$ vs. $N_{H}/A_{V}$ relation in Krühler et al. (2011). As in the case of GRB 160625B, we combined our Chandra data of GRB 160509A with the data from the Swift/XRT light curve repository converted to 5 keV flux densities.

The CIRCE $H$-band fluxes were converted to the narrower F160W filter assuming $\beta=-0.6$. The F160W and X-ray data and our power-law fits are presented in Figure 6, and the parameters of the fits are listed in Table 7. For our power law fits we ignore the data points before $\sim 0.5$ d ($4\times 10^{4}$ s), as the early X-ray light curve may contain a plateau and/or a flare; see Figure 6. In this case the smooth- and sharp-break scenarios give similar results: the best fit for the post-break decline for $\omega=3$ is $\alpha_{2,\mathrm{X}}=-1.98\pm 0.10$ and for $\omega=10$, $\alpha_{2,\mathrm{X}}=-1.96\pm 0.09$. The jet-break times, $3.2\pm 0.9$ d and $3.7\pm 0.8$ d, respectively, are consistent with each other as well.

In the radio, we obtained the fluxes at 6 and 9 GHz at the epochs earlier than 79.9 d by power-law interpolation between observed fluxes – our measurements at 36.9 d and those published in Laskar et al. (2016) at earlier times. We then fitted a single power law to the points where the reverse shock should no longer dominate the radio flux (i.e. $\geq 10$ days; Laskar et al., 2016). The resulting decline slopes are $\alpha_{\mathrm{6GHz}}=-0.91\pm 0.11$ and $\alpha_{\mathrm{9GHz}}=-0.92\pm 0.13$. Since the reverse shock may still be contributing a non-negligible fraction of the flux at 10 d, we also performed the fit without this epoch. The results are consistent but less constraining: $\alpha_{\mathrm{6GHz}}=-1.07\pm 0.18$ and $\alpha_{\mathrm{9GHz}}=-0.92\pm 0.21$. The slopes at other frequencies between 5 and 16 GHz, fitted from 10 to 20 d, are all consistent with these, ranging from $-0.80\pm 0.10$ (7.4 GHz) to $-1.02\pm 0.04$ (8.5 GHz). In F160W and/or $H$, we only have two points and an upper limit; therefore we simply measure the decline assuming a single power law. As the first point at 5.8 d is after the jet break time we obtained from the X-ray fit, there should be no significant deviation from a single power law. The measured decline is $\alpha_{2,\mathrm{F160W}}=-2.09\pm 0.10$, consistent within $1\sigma$ with the X-ray decline.

Using boxfit, we again fitted the light curve at three different values of $\xi$: 1, 0.1 and 0.01. As with the power-law fits, the X-ray points before 0.6 d were ignored, since boxfit cannot accommodate continuous energy injection. Radio points with a significant reverse shock contribution were also ignored (i.e. $<10$ d; at frequencies $<5$ GHz also 10.03 d; see Laskar et al., 2016). We ran boxfit with the boosted-frame wind-like CBM model (with both strong and medium boost) and a lab-frame model with ISM-like CBM, as the lack of optical data makes it difficult to distinguish between different CBM profiles (although the ISM scenario is tentatively favored by Laskar et al., 2016). However, as shown in Figure 7, our fits in a wind CBM do not reproduce the jet break clearly detected in the X-ray light curve. Even with the parameters in Laskar et al. (2016), the break only appears at $\sim 100$ d and the X-ray fit is much worse than with an ISM-type CBM. Thus the analytical model and boxfit seem to disagree on how the jet behaves in a wind-type CBM, and we concentrate on the ISM scenario. The best ISM fits are shown in Figure 8; Figure 9 shows the SED at three post-break epochs along with specra produced by boxfit at these epochs. Our resulting best-fit parameters are summarized in Table 8. These fits (including the wind fits) again fail to match the observed shape of the radio light curve, although the amplitude of the flux can be reproduced at some epochs.

4 Discussion

4.1 The shape of the break

In the X-ray, we find little difference in the reduced $\chi^{2}$ values of the fits between a sharp and a smooth break for GRB160625B. In the optical, however, fixing $\omega=3$ results in a visible and significant residual of $4.2\sigma$ at 140.2 d, while fixing $\omega=10$ results in a residual of $1.5\sigma$. The reduced $\chi^{2}$ of the latter fit is also slightly smaller. In the optical light curve, one can see either one slight bump or two, depending on the break time. These deviations from a perfect power law may disturb the fit and cause the high $\chi^{2}$ values, which suggests that one should also try only using the post-break points. Simply fitting a single power law to the points after 26.5 d results in consistency with the $\omega=10$ case. We thus conclude that while both values of $\omega$ remain plausible, a sharp break with $\omega=10$ is more likely. A sharp break also implies a small viewing angle $\theta_{obs}$ (Ryan et al., 2015), which is compatible with the boxfit results for this burst.

The post-jet-break decline of GRB 160625B has been previously estimated to be $f_{\nu}\propto t^{\alpha_{2}}$, where generally $\alpha_{2}\sim-2.3$ and its error roughly 0.5 (Alexander et al., 2017; Fraija et al., 2017; Lü et al., 2017). These estimates are largely consistent with both sharp and smooth breaks (and with our results listed in Table 5). However, all of these results are based on observations no later than $\sim 50$ d from the burst ($\sim 2.5\times t_{j}$, compared to our latest observations at $\sim 7\times t_{j}$), and their post-break fluxes mostly include relatively large uncertainties. In addition, Troja et al. (2017) obtained a more precise post-break slope of $\alpha_{2}=-2.57\pm 0.04$ and Strausbaugh et al. (2018) obtained $\alpha_{2,\mathrm{optical}}\approx 1.6$ and $\alpha_{2,\mathrm{X}}=-2.06\pm 0.22$, but their optical slope is inconsistent with our later-time optical data in both cases.

Troja et al. (2017) placed their estimate of the jet break at 14 d, during the ’bump’ in the light curve between $\sim 8$ d and $\sim 16$ d. Using the same data, Strausbaugh et al. (2018) suggested a break at 12.6 d at the peak of the bump, which they took as brightening of the jet toward its edges. However, our later-time data require a later break and a steeper $\alpha_{2}$, leading us to suggest the bump may still be due to angular brightness differences or perhaps the result of density fluctuations in the CBM, but not necessarily a sign of a bright edge – and seemingly not simultaneous with a true jet break. The bump is not seen in the X-rays, which is also consistent with a density fluctuation, as the flux above $\nu_{c}$ is insensitive to ambient density (Kumar, 2000). Strausbaugh et al. (2018) also suggest that a slowly changing spectral slope in the optical bands indicates a gradual cooling transition instead of a $\nu_{c}$ break in the spectrum, and that the optical spectrum eventually becomes consistent with $\beta\sim-1.1$, i.e. the slope above $\nu_{c}$, which would disfavor a CBM density fluctuation because of this insensitivity. We, however, measure $\beta=-0.86\pm 0.04$ between F606W and 5 keV at 141 d, suggesting that $\nu_{c}$ is still above optical frequencies but below X-ray at this time. Thus we cannot rule out either scenario for the bump, but we can place the jet break at an epoch after the bump.

In the case of GRB 160509A, the $\chi^{2}$ values of the fits with different $\omega$ are close to equal and the post-break slopes are in agreement. A higher $\theta_{obs}$ results in a softer break (Ryan et al., 2015), so in this case, considering that $\theta_{obs}\sim 0.6\theta_{j}$ (from boxfit), one would expect the break to be softer than for GRB 160625B where $\theta_{obs}$ is much smaller or close to zero. One can attempt to resolve this by finding inconsistencies in estimates of $p$ based on the pre-break light curve and spectrum. The X-ray spectrum, with a slope of $\beta=-1.07\pm 0.04$, is consistent with $p\approx 2.2$ and with $\nu_{c}$ being below the X-ray band (Laskar et al., 2016). As a result, we can use $\alpha=(2-3p)/4$ independent of the CBM distribution (Granot & Sari, 2002); in the case of $\omega=3$ we obtain $p=2.08\pm 0.14$ and for $\omega=10$, $p=2.27\pm 0.08$. While the former is closer to the measured post-break decline, both values are consistent with 2.2.

4.2 Physical implications

4.2.1 GRB 160625B

Based on the well-constrained pre-break light curve of the afterglow of GRB 160625B, one can estimate the electron energy distribution index $p$: below the cooling frequency $\nu_{c}$, in the case of a wind-like CBM, $\alpha_{\mathrm{wind}}=(1-3p)/4$, while for a constant-density CBM similar to the interstellar medium (ISM), $\alpha_{\mathrm{ISM}}=3(1-p)/4$ (Granot & Sari, 2002). Thus, in the optical, one obtains $p=1.63\pm 0.02$ in the wind case and $p=2.29\pm 0.02$ in the ISM case. Above $\nu_{c}$, in both cases $\alpha=(2-3p)/4$. Comparing the optical and X-ray spectra and fluxes Alexander et al. (2017) argue that $\nu_{c}$ lies below the X-ray frequencies after $\sim 1.2\times 10^{4}$ s, and thus the early X-ray light curve gives us $p=2.29\pm 0.06$. This is also consistent with the spectrum below the X-ray frequencies (Alexander et al., 2017), and thus, as the $p$ values in the wind scenario are mutually inconsistent, an ISM-like density profile is favored. Fraija et al. (2017) infer a transition from wind-like to ISM-like CBM at $\sim 8000$ s.

When only taking into account the relativistic visible-edge effect (Mészáros & Rees, 1999), the slope of the decline is expected to steepen in the jet break by a factor of $t^{-3/4}$ in a constant-density CBM. In the $\omega=10$ case, the difference between the pre- and post-break power laws is $\Delta\alpha_{\mathrm{F606W}}=-1.00\pm 0.08$ in the optical and $\Delta\alpha_{\mathrm{X}}=-0.99\pm 0.16$ in the X-ray. Thus a $t^{-3/4}$ factor can be ruled out in the optical at a $>3\sigma$ level (although in the X-ray, only at a $\sim 1.5\sigma$ level). The difference is larger in the $\omega=3$ case ($>4\sigma$ and $>2\sigma$ respectively), and therefore a simple edge effect is inconsistent with our observations regardless of whether the break is sharp or smooth (the $t^{-1/2}$ factor from a wind-like CBM is, of course, even less plausible).

If one assumes a smooth break ($\omega=3$), both the optical and X-ray post-break decline rates are consistent with the form $f_{\nu}\propto t^{-p}$, for $p\approx 2.3$, as expected from exponential lateral expansion (Rhoads, 1999; Sari et al., 1999). At first glance, the favored sharp-break scenario seems to make GRB 160625B inconsistent with a $f_{\nu}\propto t^{-p}$ decline in the optical band (the X-ray slope is still consistent with it) and would seem to require another physical mechanism. One explanation could be that the true jet break is due to a combination of the visible-edge effect and more limited lateral expansion. The steepening in both bands is a factor of $t^{-1}$, steeper than the $t^{-3/4}$ expected from the edge effect (Mészáros & Rees, 1999; Panaitescu & Mészáros, 1999), and the resulting $\alpha_{2}$ values are only consistent within $2\sigma$, while the full exponential lateral expansion scenario described by Rhoads (1999) should result in identical slopes. In some numerical simulations, lateral expansion has been found to initially involve only the outer layer of the jet carrying a fraction of its energy, and the bulk of the material remains unaffected for some time (van Eerten & MacFadyen, 2012), while the results of Rhoads (1999) require the assumption that the entire jet expands at the speed of sound. On the other hand, completely ignoring the lateral expansion was found to result in insufficient steepening across the jet break. This scenario seems consistent with our results.

A complication was noted by Gompertz et al. (2018), who find that using different synchrotron relations to estimate $p$ (such as using the spectral index or the pre- or post-break decline) typically results in different estimates, with an intrinsic scatter on the value of $p$ of $0.25\pm 0.04$ (we will denote this as $\sigma_{p}$). They argue this is probably caused by emission from GRB afterglows not behaving exactly as the rather simplified analytical models predict141414We note that the inconsistency between $p$ values derived from the optical and X-ray pre-break slopes assuming a wind-type CBM is $>2\sigma_{p}$, so an ISM-like density profile is still favored.. Taking this scatter into account, both $\alpha_{2,\mathrm{F606W}}$ and $\alpha_{2,\mathrm{X}}$ in the $\omega=10$ case (or simply using only the $>26.5$ d points and a single power law) are in fact consistent within $\approx 1\sigma_{p}$ with $f_{\nu}\propto t^{-p}$. Thus lateral expansion at the speed of sound can still account for the observed late-time decline. Using closure relations for both the light curve and the spectrum, Gompertz et al. (2018) found a best fit of $p=2.06\pm 0.13$ for GRB 160625B, which is consistent with our results in both bands within $\sigma_{p}$. In any case, for this burst some form of lateral expansion is required, and the edge effect alone is insufficient.

We can also attempt to use the results from boxfit to determine if the magnetar spin-down power source is consistent with the GRB. The rotational energy that can be extracted from a millisecond magnetar is (Lü & Zhang, 2014; Kumar & Zhang, 2015)

[TABLE]

where $M$ is the mass, $R$ the radius and $P_{0}$ the initial spin period of the newborn magnetar. Metzger et al. (2015) placed a limit of $\sim 1\times 10^{53}$ erg on the maximum energy of a newborn magnetar in extreme circumstances (in terms of mass and spin period). Therefore the energy requirements of all the fits from boxfit may technically be achievable with the magnetar model, but with the (more realistic) low $\xi$ values the required energy approaches or exceeds even this maximum limit. The exceptionally high $E_{\mathrm{iso}}$ can be due to a relatively narrow jet and a lower explosion energy instead, but this requires a high $\xi$ that is inconsistent with simulations by Sironi & Spitkovsky (2011) and Warren et al. (2018) – the best fit at $\xi=1$ also results in an extremely low density more typical to intergalactic environments. We do point out a caveat that the parameters of the best fits show a non-monotonic dependence on $\xi$, with notable degeneracy between parameters.

We have attempted to use boxfit to estimate errors for the best-fit parameters as well. However, as a result of what seems to be a bug in the error estimation routine of boxfit (G. Ryan and H. van Eerten, private communication), some of the errors are clearly incorrect and, therefore, we have not included errors in our Table 6. This mostly manifests as error limits that either do not include the best fit or where the best-fit value of a parameter is always the lower limit as well151515In other cases, such as the $\xi=0.1$ case of GRB 160625B, the errors are seemingly reasonable ($p=2.05\pm 0.01$, $E_{\mathrm{K,iso}}=1.4_{-1.3}^{+1.2}\times 10^{54}$ erg, $\epsilon_{e}=0.25_{-0.13}^{+0.10}$, $\epsilon_{B}=3.0_{-2.0}^{+106.3}\times 10^{-4}$, $n=0.18_{-0.15}^{+0.58}$ cm*-3*, $\theta_{j}=0.14\pm 0.03$ rad and $\theta_{obs}=1.1_{-1.1}^{+5.9}\times 10^{-3}$ rad) and the relative ranges of each parameter comparable to those found by Alexander et al. (2017). These values thus give an indication of how well each parameter is constrained.. We also note that, as the shape of the radio light curve is not well reproduced in any of our fits, error limits could be misleading in any case. As a consistency check for the rest of the code, we have run boxfit using the Alexander et al. (2017) forward shock parameters, which are similar to our $\xi=1$ results. The output light curves and spectra are similar to the analytical ones and reproduce the early behavior of the afterglow well, although post-break fluxes are somewhat under-predicted.

We also note that the 6.1 GHz light curve of GRB 160625B is not successfully reproduced by boxfit, and the jet model struggles to explain the late slope of $\alpha_{\mathrm{6.1GHz}}=-1.08\pm 0.11$ and the lack of an observed jet break. At low $\xi$ values, the boxfit fit is somewhat better, but only if one ignores the 22.5 d point, where a low-frequency scattering event by an intervening screen, suggested by Alexander et al. (2017), may contribute to the flux. The radio SED shows a peak centered at 3 GHz between 12 and 22 d, which then disappears. Even so, the fit at 22 GHz is worse at all $\xi$ values. At 48 d, the radio SED is consistent with being entirely flat, which is only plausible in the standard model around a very smooth $\nu_{m}$ break. While the low $\xi$ fits do place the $\nu_{m}$ passage at roughly this time, the peak in the boxfit spectrum is too sharp, and in earlier spectra the lowest frequencies must then be brightened by a factor of ten or so by the proposed scattering. The shape may instead be altered by another emission source contributing to the spectrum (see below).

Theoretically expected post-break values in the slow-cooling scenario ($\nu_{m}<\nu_{c}$) are $-p$ or $-1/3$, depending on which side of $\nu_{m}$ the band is located (Rhoads, 1999). As the jet break is a geometric effect, we should see it in every band, but this is not the case: we can set a limit of $t_{j,\mathrm{6.1GHz}}\gtrsim 10\times t_{j,\mathrm{F606W}}$. The possibilities given by the standard jet model that are consistent with the slope are:

• •

Post-break, $\nu_{c}<\nu_{m}$, i.e. fast-cooling: $\alpha_{\mathrm{6.1GHz}}$ is consistent with the expected decline of $\alpha_{2}=-1$. However, the measured $\alpha_{1,\mathrm{F606W}}=-0.96\pm 0.01$ does not match the pre-break decline expected at any frequency in this scenario.

• •

Pre-break, $\nu_{m}<6.1~{}\mathrm{GHz}<\nu_{c}$: $\alpha_{\mathrm{6.1GHz}}$ is consistent with $p=2.4$ and $\alpha=3(1-p)/4=-1.05$ (Granot & Sari, 2002). However, the spectral index between radio and optical is $-0.35\pm 0.03$ at 22 d and $-0.49\pm 0.01$ at 140 d, which is intermediate between the indices expected above and below $\nu_{m}$ (respectively, $(1-p)/2\approx-0.65$ and $1/3$) and thus implies that $\nu_{m}>6.1$ GHz at 140 d, or that $p\approx 2.0$.

• •

A transition to a non-relativistic flow, $\nu_{m}<6.1~{}\mathrm{GHz}<\nu_{c}$: the expected slope is $(21-15p)/10$ (van der Horst, 2007), resulting in $p=2.12\pm 0.08$, which is consistent with our estimate within $\sigma_{p}$. However, such a transition is not seen in the optical or X-ray bands.

The LGRB population has been observed to be comprised of a radio-quiet and a radio-loud population, where the radio-quiet GRBs are incompatible with a simple sensitivity effect and indicate an actual deficit in radio flux compared to theory (Hancock et al., 2013). Lloyd-Ronning et al. (2019) further argued that the two populations originate in different progenitor scenarios. This deficit in radio flux implies some mechanism that suppresses the expected synchrotron emission at radio frequencies. Since our findings indicate that the radio light curve of GRB 160625B (and GRB 160509A; see below) is incompatible with the higher frequencies, the source of the radio emission that we do see may not be the same as that of the optical and X-ray synchrotron emission. This seems to suggest that even in (at least) some radio-loud GRBs, the same mechanism may be in effect. Furthermore, if the radio emission is generated by another source, this source is not active in the radio-quiet GRBs for some reason. We have run the boxfit fitting code with $\xi=1$ and all radio fluxes divided by ten to investigate if the standard model allows suppression of the radio flux simply through adjusting the parameters. The resulting best fit over-predicts all radio fluxes by at least a factor of a few at all times. This implies a caveat that, at least in some cases, including another, dominant radio source without an additional suppression mechanism may over-predict the radio flux. Another caveat with this is that, unless the second component is coupled to the ’main’ source, getting a total radio flux compatible with one component may require fine-tuning. If such a mechanism is widespread, one would expect some GRBs to have radio fluxes unambiguously too high for the standard model, which, to our knowledge, has not been seen.

One explanation for the ’extra’ radio source, with its lack of a jet break and the requirement of $6.1~{}\mathrm{GHz}>\nu_{m}$, could be a two-component jet, where a narrow jet core is surrounded by a cocoon with a lower Lorentz factor (Berger et al., 2003; Peng et al., 2005), resulting in a different source with different physical parameters dominating the radio emission, and thus a different break time and $\nu_{m}$. This does not result in a deficit in radio synchrotron flux, only an inconsistency between the light curve shape and the standard model. For an on-axis or slightly off-axis burst ($\theta_{obs}<\theta_{j,\mathrm{narrow}}$), the wider component would not contribute significantly to the optical light curve if its kinetic energy is lower than that of the narrow component (Peng et al., 2005). This may also affect the required energy, but without robust modeling it is difficult to say whether the consistency with a magnetar energy source would change.

Strausbaugh et al. (2018) suggested a scenario where a very smooth cooling transition (i.e. not a normal spectral break) is moving through the optical and infrared frequencies, starting at early times, and the optical spectrum becomes consistent with $\nu>\nu_{c}$ by $\sim 50$ d. This would indicate a unique cooling behavior inconsistent with the standard expectations. The observed lack of evolution of the Swift spectra until 30 d implies that the X-ray spectral slope $\beta_{X}$ is not the result of a $\nu_{c}$ break right below the X-ray frequencies, as this would require the spectrum to soften over time to its slope at $\nu\gg\nu_{c}$. Furthermore the optical-to-X-ray index is observed to gradually steepen and eventually become similar to $\beta_{X}$. This is qualitatively consistent with the reddening of the optical spectrum noted by Strausbaugh et al. (2018). In addition, $\beta_{X}$ indicates a different $p$ than the X-ray light curve; this agrees with the implication of Gompertz et al. (2018) that some physics is missing or simplified in the relevant closure relations. Another possible explanation is that a Klein-Nishina correction (Nakar et al., 2009) is needed above $\nu_{c}$; this can result in $\beta=3(1-p)/4$, which would imply $p\approx 2.1$. This harder spectrum is expected to dominate when the $\epsilon_{e}/\epsilon_{B}$ ratio is high, which would fit the low-$\xi$ boxfit results.

4.2.2 GRB 160509A

In the case of GRB 160509A, the change in X-ray decay slope across the break, $\Delta\alpha_{\mathrm{X}}=-0.75\pm 0.11$ for a sharp break and $\Delta\alpha_{\mathrm{X}}=-0.92\pm 0.15$ for a smooth break. Thus we cannot exclude the $t^{-3/4}$ factor expected from the edge effect alone in an ISM-like medium. The $t^{-1/2}$ factor expected in the case of a wind medium is inconsistent with the observations at a $2.3\sigma$ or $3\sigma$ level, depending on $\omega$. However, when considering the intrinsic $p$ scatter of $\sigma_{p}=0.25$ (Gompertz et al., 2018), $\alpha_{2,{\mathrm{X}}}$ is also consistent with a $f_{\nu}\propto t^{-p}$ decline. Thus we cannot say conclusively whether lateral expansion is important in the jet of GRB 160509A, but it does not seem necessary. In the IR, the measured slope of $\alpha_{2,\mathrm{F160W}}=2.09\pm 0.10$ is marginally consistent ($1.1\sigma$) with $p\approx 2.2$, but a lack of pre-break data prevents us from determining $\Delta\alpha_{\mathrm{F160W}}$.

The decline of the afterglow in the radio after 10 d is about $f\propto t^{-0.9}$ at both 6 and 9 GHz (and consistent with this at other frequencies where fewer points are available). This is again inconsistent with the expected post-jet-break slope of $-p$ or $-1/3$ in the slow-cooling case, respectively above and below the characteristic synchrotron frequency $\nu_{m}$ (Rhoads, 1999). As with GRB 160625B, we list the possibilities consistent with this decline, allowed by standard jet theory:

• •

Post-break, $\nu_{c}<\nu_{m}$, i.e. fast-cooling: $\alpha=-1$ is expected and consistent with $\alpha_{\mathrm{radio}}$, but this scenario is incompatible with the measured IR-to-X-ray spectral index $-0.74\pm 0.09$ at 35 d, as the expected index is $-p/2\approx-1.1$ (a photon index consistent with this is indeed seen in the X-ray at earlier times according to the UKSSDC Swift Burst Analyser161616http://www.swift.ac.uk/burst_analyser/00020607/ – $<\beta_{X}>=1.06\pm 0.04$ between 1 and 10 d – indicating that $\nu_{c}$ is still between X-ray and optical frequencies and $\nu_{c}>\nu_{m}$ at 35 d).

• •

Pre-break, $\nu_{m}<6~{}\mathrm{GHz}<\nu_{c}$, ISM-like CBM: $\alpha_{\mathrm{radio}}$ is consistent with the expected decline ($\alpha=3(1-p)/4=-0.9$ assuming $p=2.2$; Granot & Sari, 2002), but the observed spectral index of $-0.40\pm 0.01$ between F160W and 9 GHz at 35 d implies $\nu_{m}>9$ GHz.

• •

A transition to a non-relativistic flow, $\nu_{m}<6~{}\mathrm{GHz}<9~{}\mathrm{GHz}<\nu_{c}$: the expected slope is $(21-15p)/10$, resulting in $p=2.01\pm 0.08$ – again consistent with our estimate within $\sigma_{p}$. However, such a transition is not seen in the X-ray light curve, which continues to evolve consistently with a relativistic flow.

The best boxfit fit at $\xi=1$ places a smooth, and thus off-axis, jet-break at a later time, around 35 d in all bands, in which case the radio light curve would include contamination from the reverse shock at early times, changing the decline slope (see Figure 8). This is because boxfit attempts to fit a model with a late break to $f_{\nu}\propto t^{-p}$ in order to match the radio light curve, which has no observed break. It is incompatible with the broken power-law fit with $t_{j}\sim 3.5$ d, though, and at lower, more realistic values of $\xi$ the break is placed at an earlier time. This scenario is therefore not supported. Instead, for GRB 160509A we can place a lower limit of $t_{j,\mathrm{radio}}\gtrsim 20\times t_{j,\mathrm{X}}$ based on the broken power-law fit. The situation in the radio frequencies is thus qualitatively very similar to that of GRB 160625B, and the same mechanisms may well be in effect.

We note that Kangas & Fruchter (2019) are, in fact, able to get a plausible fit to the GRB 160509A radio light curve using an analytical fit based on the standard model, but only if the light curve smoothly turns over to a $t^{-p}$ decline immediately after the last radio epoch, which is suspicious as their sample contains several GRBs with no unambiguously observed radio breaks, and many cases where the standard model does not fit the radio light curve. We also note that as Laskar et al. (2016) showed, the radio SED seems to remain roughly flat after the reverse shock influence on the light curve fades ($\sim 20$ d), which might again be caused by another emission component. As boxfit also disagrees with this analytical model, one or the other is in doubt. The issue will be addressed in more detail in the upcoming revised version of that paper.

A boxfit simulation using the FS parameters of the Laskar et al. (2016) analytical model agrees fairly well with the X-ray data and reproduces the rough magnitude of the radio light curve but not its shape (assuming some RS contribution not accounted for by boxfit), but over-predicts the IR flux by a factor of about 10. Their IR light curve does not include host subtraction, and they fit for extinction as another free parameter in their model. Our host subtraction allows us to estimate the extinction and true IR fluxes independently, and in light of this the Laskar et al. (2016) model becomes incompatible with the IR data. Thus our boxfit results provide a better reproduction of the light curve in the IR. However, again, the fit parameters show a non-monotonic dependence on $\xi$. As with GRB 160625B above, boxfit was clearly unable to produce meaningful error bars for the parameters in some cases, and these are not included in Table 8171717Again, the ranges of each parameter at $\xi=0.1$, which are large but not obviously incorrect, may provide some indication of how well each parameter is constrained: ($p=2.13_{-0.01}^{+0.02}$, $E_{\mathrm{K,iso}}=3.8_{-3.4}^{+24.8}\times 10^{53}$ erg, $\epsilon_{e}=0.45_{-0.20}^{+0.31}$, $\epsilon_{B}=1.7_{-0.7}^{+2.7}\times 10^{-5}$, $n=18_{-18}^{+1530}$ cm*-3*, $\theta_{j}=0.20_{-0.16}^{+0.18}$ rad and $\theta_{obs}=0.12_{-0.12}^{+0.21}$ rad). – and, as the radio light curve is again problematic for the fit, would be misleading in any case.

Keeping in mind the caveats associated with our best boxfit fits, we can use them to estimate the energy requirements. The geometry-corrected jet energy $1.8\times 10^{51}$ erg at $\xi=1$ is well below the maximum rotational energy of a millisecond magnetar (see Section 4.2.1) Once again, we deem the lower $\xi$ values more realistic based on simulations (Sironi & Spitkovsky, 2011; Warren et al., 2018) and the fact that $\xi=1$ results in an extremely low density. Low $\xi$ values require energies around $\sim 3\times 10^{52}$ erg, which again strains the magnetar spin-down model but does not rule it out. Thus GRB 160509A also seems compatible with a magnetar power source.

For both GRBs considered here (Tables 6 and 8), but especially for GRB 160509A, the efficiency $\eta$ of the prompt $\gamma$-ray emission depends on the value of $\xi$ used in the fitting, but not monotonically: with $\xi=0.01$ one obtains a much lower value for $\eta$ than otherwise. In both cases, the difference in $\chi^{2}$ between the $\xi=0.1$ and $\xi=0.01$ fits is minimal, and in the case of GRB 160509A, so is the difference between $\xi=1$ and $\xi=0.01$; thus, we cannot reliably distinguish between these scenarios. In the literature, it is commonly assumed that $\xi=1$, and high values of $\eta$ are obtained: for example, Lloyd-Ronning & Zhang (2004) find values as high as $\eta\sim 1$ depending on $E_{iso}$, and mostly $\eta\gtrsim 0.3$. Such a high efficiency is used as a criterion for successful models of prompt emission: e.g. the internal shock mechanism tends to result in $\eta\lesssim 0.1$ (Kumar & Zhang, 2015, and references therein). Our results may indicate that, if very low values of $\xi$ are more realistic (Warren et al., 2018), one should not dismiss models based on low efficiency.

5 Conclusions

We have presented our late-time optical, radio and X-ray observations of the afterglows of GRB 160625B and GRB 160509A. We have fitted broken power law functions to the data, combined with light curves from the literature, to constrain the jet break time and the post-jet-break decline, and used the numerical afterglow fitting software boxfit (van Eerten et al., 2012) to constrain the physical parameters and energetics of the two bursts. Our conclusions are as follows:

Regardless of the sharpness of the GRB 160625B jet break, we find that the effect of the jet edges becoming visible as the jet decelerates is alone insufficient to explain the post-jet-break light curves. A full lateral expansion break onto a $t^{-p}$ decline is also inconsistent with the favored sharp break. The light curve behavior seems qualitatively consistent with the edge effect combined with only a fraction of the jet expanding at the speed of sound (van Eerten & MacFadyen, 2012). It is also possible that an intrinsic scatter in the electron Lorentz factor distribution index $p$ exists, the result of simplified synchrotron theory and closure relations that do not necessarily reflect the true complexity of the emission region (Gompertz et al., 2018). This scenario combined with lateral expansion is also consistent with our results. For GRB160509A we are unable to exclude any of the considered scenarios due to the scarcity of the available data.

Based on the best fits from boxfit, the geometry-corrected energy requirements of both GRBs are consistent with a magnetar spin-down energy source – albeit only in extreme cases when the ’participation fraction’ (fraction of electrons accelerated into a non-thermal distribution) is fixed at $\xi=0.1$ or $\xi=0.01$, requiring energies of $\sim 3\times 10^{52}$ or even $\sim 10^{53}$ erg. As simulations have shown these lower fractions to be more realistic (e.g. Warren et al., 2018), it seems that magnetar spin-down alone struggles to produce the required energies unless the nascent magnetar has extreme properties (Metzger et al., 2015).

However, neither boxfit nor analytical relations from standard jet theory (e.g. Rhoads, 1999; Granot & Sari, 2002) can provide a good fit to the radio data of either GRB, which are consistent with a single power law that requires the jet break to occur much later in radio than in the other bands. Both GRBs also show an almost flat radio SED at relatively late times (tens of days; see Laskar et al., 2016; Alexander et al., 2017). The higher frequencies do conform to expectations from the jet model, though. This might be the result of a multi-component jet, but that would require the wide component of the jet to dominate the light curve, and simultaneously suppressed flux from the narrow component. A similar behavior (a radio decline described by a single power law with $\alpha=-1.19\pm 0.06$ until $\sim 60$ d) was recently reported for GRB 171010A by Bright et al. (2019). We explore this problem further in a companion paper (Kangas & Fruchter, 2019), and find that these GRBs are not exceptional in this regard.

We thank the anonymous referee for comments that improved the paper. We also thank Hendrik van Eerten, Geoffrey Ryan, Alexander van der Horst, Paz Beniamini and Chryssa Kouveliotou for helpful discussions. A.Co. acknowledges support from the National Science Foundation via CAREER award #1455090. A.dU.P. acknowledges support from a Ramon y Cajal fellowship (RyC-2012-09975), from the Spanish research project AYA2017-89384-P, and from the State Agency for Research of the Spanish MCIU through the ”Center of Excellence Severo Ochoa” award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). A.P. wishes to acknowledge support by the European Research Council via the ERC consolidating grant #773062 (acronym O.M.J.). A.Cu. acknowledges the support of NASA MIRO grant NNX15AP95A. Based on observations made with the NASA/ESA Hubble Space Telescope (programme GO 14353, PI Fruchter), obtained through the data archive at the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. Support for this work was also provided by the National Aeronautics and Space Administration through Chandra Award Number 17500753, PI Fruchter, issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester as well as observations made with the Gran Telescopio Canarias (GTC), installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. Development of CIRCE was supported by the University of Florida and the National Science Foundation (grant AST-0352664), in collaboration with IUCAA. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Development of the boxfit code was supported in part by NASA through grant NNX10AF62G issued through the Astrophysics Theory Program and by the NSF through grant AST-1009863. Simulations for boxfit version 2 have been carried out in part on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP) of the research cooperation ”Excellence Cluster Universe” in Garching, Germany.