The Late-time Afterglow Evolution of Long Gamma-Ray Bursts GRB 160625B and GRB 160509A

Tuomas Kangas; Andrew S. Fruchter; S. Bradley Cenko; Alessandra Corsi; Antonio de Ugarte Postigo; Asaf Pe’er; Stuart N. Vogel; Antonino Cucchiara; Benjamin Gompertz; John Graham; Andrew Levan; Kuntal Misra; Daniel A. Perley; Judith Racusin; Nial Tanvir

doi:10.3847/1538-4357/ab8799

1. Introduction

Gamma-ray bursts (GRBs) are among the most luminous transient events in the universe. Through their association with broad-lined SNe Ic (e.g., Iwamoto et al. 1998; Woosley & Bloom 2006; Hjorth & Bloom 2012), long GRBs (LGRBs; duration of the prompt γ-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 γ-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_ν ∝ t^α.

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 an LGRB; e.g., Racusin et al. 2009). This results in an achromatic jet break in the afterglow light curve when the relativistic beaming angle (Γ⁻¹, where Γ 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, behavior consistent with an isotropic fireball model. This phenomenon is believed to steepen the decline slope α 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(γ) ∝ γ^−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, θ_j ≲ 3° (Granot & Piran 2012; Lyutikov 2012). At even later times, the jet is expected to be better described as a nonrelativistic 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 ≈ 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., ∼t⁻², 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–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_iso and the peak of the energy spectrum (Amati et al. 2002), and can haveisotropic-equivalent energies on the order of 10⁵⁴ 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⁵² 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 57,564.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 γ-ray bursts ever observed with E_iso ∼ 3 × 10⁵⁴ erg (Wang et al. 2017; Zhang et al. 2018) and a well-studied object with a multifrequency 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 ∼1 day, 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 57,517.4; Longo et al. 2016a, 2016b; Roberts et al. 2016) at a redshift of z = 1.17 (Tanvir et al. 2016). With E_iso = 8.6 ± 1.1 × 10⁵³ 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 afterglow,¹⁵ 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, we present our conclusions in Section 5. All magnitudes are in the AB magnitude system (Oke & Gunn 1983), and all error bars correspond to 1σ confidence intervals. We use the cosmological parameters H₀ = 69.6 km s⁻¹ Mpc⁻¹, Ω_m = 0.286, and Ω_Λ = 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 days) and 2016 November 13 (140.2 days). A template image of the host galaxy was created by combining images obtained with the same setup on 2017 November 6 (498.3 days) and 11 (503.6 days). At this time, the contribution of the afterglow itself was a factor of ∼13 fainter than at 140 days, assuming a f_ν ∝ t^α decline where α = −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 days) and 2016 June 3 (24.8 days). Late-time imaging of GRB 160509A was done using HST/WFC3 and the F110W and F160W filters on 2016 June 13 (35.3 days); template images of the host galaxy in these filters were obtained on 2017 July 5 (422.1 days), when, assuming α = −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 at 10.17909/t9-yvpg-xb33 (GRB 160625B) and 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 drizzlepac¹⁶ package in pyraf.¹⁷ The two epochs of GRB 160625B in 2017 November were combined into one template image. Subtraction of the template images and aperture photometry of the afterglows were done using iraf.¹⁸ 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 catalog¹⁹ (Skrutskie et al. 2006), and aperture photometry was performed using standard iraf tasks. At 24.8 days, we were unable to detect the afterglow and only obtained a (3σ) limit of H ≥ 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 α = −2.0 ± 0.2 (obtained from a single power-law fit to uncorrected >25 days data, with errors rounded up to be conservative), the afterglow flux present in the template image was estimated to be 2.0 ± 1.0% of the flux at 71.5 days or 7.5 ± 2.6% of the flux at 140.2 days, 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% of the 35.3 days 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.

**Figure 1.** Afterglow and host galaxy of GRB 160625B in the F606W band. Panel (a): the afterglow and the host galaxy at 71.5 days; panel (b): 140.2 days; panel (c): the combined template at ∼500 days; panel (d): the template-subtracted image at 71.5 days; and panel (e): the subtraction at 140.2 days. North is up and east is to the left in all panels. The black north–south line corresponds to one arcsecond. The afterglow location is indicated with red tick-marks.
Download figure:
Standard image High-resolution image

**Figure 2.** Afterglow and host galaxy of GRB 160509A in the F160W band. Panel (a): the afterglow and the host galaxy at 35.3 days; panel (b): the template at 422.1 days; and panel (c): the template-subtracted image at 35.3 days. North is up and east is to the left in all panels. The black north–south line corresponds to one arcsecond. The afterglow location is indicated with red tick-marks. The afterglow is very weak compared to the host galaxy, making a template subtraction crucial for this target.
Download figure:
Standard image High-resolution image

Table 1. Log of Our Late-time HST/WFC3 Observations of GRB 160625B

Phase	MJD	t_exp	F606W	Corrected F606W
(days)		(s)	(mag)	(mag)
71.5	57636.4	2400	25.38 ± 0.03	25.36 ± 0.04
140.2	57705.1	4800	26.76 ± 0.06	26.67 ± 0.07
498.3	58063.2	4800	⋯	⋯
503.6	58068.5	4800	⋯	⋯

Download table as: ASCII Typeset image

Table 2. Log of Our Late-time HST/WFC3 and GTC/CIRCE Observations of GRB 160509A

Phase	MJD	${t}_{\exp ,{\rm{F}}110{\rm{W}}}$	F110W	${t}_{\exp ,{\rm{F}}160{\rm{W}}}$	F160W	${t}_{\exp ,H}$	H
(days)		(s)	(mag)	(s)	(mag)	(s)	(mag)
5.8	57523.2	⋯	⋯	⋯	⋯	3060	20.50 ± 0.17
24.8	57542.2	⋯	⋯	⋯	⋯	2100	≥21.9
35.3	57552.7	2697	27.11 ± 0.10	2797	26.07 ± 0.07	⋯	⋯
422.1	57939.5	2697	⋯	2797	⋯	⋯	⋯

Download table as: ASCII Typeset image

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 days), 2016 November 15 (142.3 days), and 2016 November 19 (146.2 days). The latter two epochs were combined to obtain the flux at 144.3 ± 2.2 days, 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 days). 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 (pimms²⁰ ) 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}_{{\rm{H}},\mathrm{MW}}=9.76\times {10}^{20}$ cm⁻² (Willingale et al. 2013), a photon index of Γ_X = 1.86, and an intrinsic absorption of ${N}_{{\rm{H}},\mathrm{int}}=2.1\,\times {10}^{21}$ cm⁻² 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}_{{\rm{H}},\mathrm{MW}}=2.12\times {10}^{21}$ cm⁻² (Willingale et al. 2013), a photon index of Γ_X = 2.07 and an intrinsic absorption of ${N}_{{\rm{H}},\mathrm{int}}=1.52\times {10}^{22}$ cm⁻², following Laskar et al. (2016). Γ_X is assumed to be constant over the light-curve break. The log of the X-ray observations and derived flux densities is presented in Table 3.

Table 3. Log of Our Late-time Chandra/ACIS-S Observations of GRB 160625B and GRB 160509A

Phase	MJD	t_exp	f_ν(5 keV)
(days)		(ks)	(erg s⁻¹ cm⁻² keV⁻¹)
160625B

69.8	57634.7	19.80	(1.47 ± 0.29) × 10⁻¹⁵
142.3	57707.3	45.84	⋯
144.3 ± 2.2^a	57709.3 ± 2.2	69.56	(3.21 ± 0.79) × 10⁻¹⁶
146.2	57711.2	23.72	⋯

160509A

42.1	57559.5	24.75	(1.38 ± 0.25) × 10⁻¹⁵

Note.

^aCombination of the 142.3 and 146.2 days epochs.

Download table as: ASCII Typeset image

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 days) and 2017 January 20 (209.0 days), and GRB 160509A in the C and X bands on 2016 June 2 (23.9 days), 2016 June 15 (36.9 days), and 2016 July 28 (79.9 days) (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).²¹ 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 FWHM 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% below 18 GHz and 10% above it.²²

Table 4. Log of Our VLA Radio Observations of GRB 160625B and GRB 160509A

Phase	MJD	ν	f_ν	Configuration
(days)		(GHz)	(μJy)
160625B

4.5	57569.4	4.8	104 ± 16	B
4.5	57569.4	7.4	454 ± 27	B
4.5	57569.4	19	278 ± 35	B
4.5	57569.4	25	204 ± 36	B
13.4	57578.3	4.8	377 ± 25	B
13.4	57578.3	7.4	310 ± 21	B
13.4	57578.3	22	163 ± 20	B
31.3	57596.2	7.4	113 ± 16	B
31.3	57596.2	22	88 ± 19	B
58.3	57623.2	6.1	75 ± 11	B
58.3	57623.2	22	52 ± 13	B
209.0	57773.9	6.1	16 ± 5	A

160509A

23.9	57541.3	6.0	80 ± 8	B
23.9	57541.3	9.0	71 ± 7	B
36.9	57554.3	5.0	50 ± 7	B
36.9	57554.3	6.9	52 ± 7	B
36.9	57554.3	8.5	41 ± 6	B
36.9	57554.3	9.5	29 ± 6	B
79.9	57597.3	6.0	27 ± 6	B
79.9	57597.3	9.0	25 ± 5	B

Note. The GRB 160625B points until 31.3 days were also reported in Troja et al. (2017) but without the calibration uncertainty.

Download table as: ASCII Typeset image

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 Sloan Digital Sky Survey $r^{\prime}$ band light curves of GRB 160625B. However, there is a slight (∼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_ν ∝ ν^−0.68 between the characteristic synchrotron frequency ν_m and the cooling frequency ν_c (Alexander et al. 2017). As the optical spectrum with β = −0.68 ± 0.07 is consistent with the expected index of β = −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 light-curve repository²³ (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

$\begin{eqnarray}&&{f}_{\nu }={f}_{\nu ,0}{\left[{\left(\displaystyle \frac{t}{{t}_{j}}\right)}^{-\omega {\alpha }_{1}}+{\left(\displaystyle \frac{t}{{t}_{j}}\right)}^{-\omega {\alpha }_{2}}\right]}^{-\tfrac{1}{\omega }}\end{eqnarray} \tag{ 1 }$

to the light curve, where t_j is the jet break time, α₁ is the pre-break power-law slope, α₂ the post-break slope, and ω 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 ω (a value of 3 was found to be consistent with most GRB observations by Liang et al. 2007, but some events were found to require a sharper break with ω = 10). The results of the fit parameters are presented in Table 5. The pre-break decline α₁ does not depend on the choice of ω; we found ${\alpha }_{1,{\rm{F}}606{\rm{W}}}=-0.96\pm 0.01$ and ${\alpha }_{1,{\rm{X}}}=-1.24\pm 0.02$ in both cases. The best fit to the post-break decline was ${\alpha }_{2,{\rm{F}}606{\rm{W}}}\,=-2.27\pm 0.13$ and ${\alpha }_{2,{\rm{X}}}=-2.40\pm 0.19$ assuming ω = 3, and ${\alpha }_{2,{\rm{F}}606{\rm{W}}}=-1.96\pm 0.07$ and ${\alpha }_{2,{\rm{X}}}=-2.23\pm 0.15$ when ω = 10. The optical and X-ray light curves and our best fits in both cases are shown in Figure 3.

**Figure 3.** Observed optical (extinction-corrected), X-ray, and interpolated 6.1 GHz light curve of the afterglow of GRB 160625B (points) and our power-law fits including the broken power laws described by Equation (1) (lines). The $r^{\prime}$ -band magnitudes from Troja et al. (2017) (solid circles) have been converted into flux density. X-ray flux densities from Swift/XRT (solid triangles) and Chandra/ACIS-S (open triangles) are reported at 5 keV. The post-break fit is better assuming ω = 10 (dotted–dashed green line), especially regarding the optical HST point at 140.2 days. The pre-break fit does not depend on the choice of ω.
Download figure:
Standard image High-resolution image

$r^{\prime} $ — **Figure 3.** Observed optical (extinction-corrected), X-ray, and interpolated 6.1 GHz light curve of the afterglow of GRB 160625B (points) and our power-law fits including the broken power laws described by Equation (1) (lines). The $r^{\prime}$ -band magnitudes from Troja et al. (2017) (solid circles) have been converted into flux density. X-ray flux densities from Swift/XRT (solid triangles) and Chandra/ACIS-S (open triangles) are reported at 5 keV. The post-break fit is better assuming ω = 10 (dotted–dashed green line), especially regarding the optical HST point at 140.2 days. The pre-break fit does not depend on the choice of ω.
Download figure:
Standard image High-resolution image

Table 5. Parameters of Our Best Smooth Broken Power-law and Single Power-Law Fits to the GRB 160625B Light Curves

Parameter	ω = 3	ω = 10	Single Power Law
${t}_{j,{\rm{F}}606{\rm{W}}}$	24 ± 3 days	17 ± 2 days	17 ± 4 days
${\alpha }_{1,{\rm{F}}606{\rm{W}}}$	−0.96 ± 0.01	−0.96 ± 0.01	−0.97 ± 0.01
${\alpha }_{2,{\rm{F}}606{\rm{W}}}$	−2.27 ± 0.13	−1.96 ± 0.07	−1.94 ± 0.13
Reduced χ²	5.5	4.4	1.8

t_j,X	27 ± 5 days	22 ± 4 days	22 ± 5 days
α_1,X	−1.24 ± 0.02	−1.24 ± 0.02	−1.25 ± 0.03
α_2,X	−2.40 ± 0.19	−2.23 ± 0.15	−2.20 ± 0.13
Reduced χ²	0.91	0.81	0.84

Note. The bump(s) between 8.5 and 26.5 days were ignored in the single power-law fits to the early and late decline.

Download table as: ASCII Typeset image

We also fitted the decline using a single power law before 8.5 days and another after 26.5 days, 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 χ² values of these fits are rather high; although, the small errors also contribute to this. The result is ${\alpha }_{2,{\rm{F}}606{\rm{W}}}=-1.94\pm 0.13$ , nearly exactly coinciding with the ω = 10 case but with a ∼2.5σ difference to ω = 3. Repeating this in the X-ray results in ${\alpha }_{2,{\rm{X}}}=-2.20\pm 0.13$ , which is also almost identical to the ω = 10 case. A simultaneous single power-law fit to both post-break light curves results in α₂ = −2.01 ± 0.09.

Assuming an achromatic break, we determined t_j by taking the weighted average of ${t}_{j,{\rm{F}}606{\rm{W}}}$ and t_j,X. In the ω = 10 case, the result is t_j = 19 ± 2 days. Assuming an instantaneous break (corresponding to $\omega =\infty$ ) between the single power-law fits, the resulting jet break times are consistent, ${t}_{j,{\rm{F}}606{\rm{W}}}\,=17\pm 4$ days and ${t}_{j,{\rm{X}}}=22\pm 5$ days, and the weighted average t_j = 19 ± 3 days. In the ω = 3 case, we obtained t_j = 25 ± 3 days.

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 and 209.0 days, 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 days. We also scaled the 7.4 GHz flux at 31.34 days assuming the same power law as at 22.5 days. Points earlier than 22.5 days 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 }_{6.1\mathrm{GHz}}=-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 Library,²⁴ 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 of 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 a CBM of the interstellar medium (ISM) type (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 ξ, i.e., the fraction of electrons accelerated by the shock into a nonthermal power-law distribution. Simulations indicate that 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 γ-ray energy E_iso = 3.0 × 10⁵⁴ 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}_{{\rm{K}},\mathrm{iso}}+{E}_{\mathrm{iso}})$ for the conversion of kinetic energy to γ-rays. These fits, however, fail to reproduce the measured power-law slope of ${\alpha }_{6.1\mathrm{GHz}}=-1.08\pm 0.11$ , instead predicting a break in the radio light curve around ∼100 days (associated with the passage of ν_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 we 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 ξ values, the fit at higher frequencies or earlier times is still somewhat worse; we show 22 GHz as an example.

**Figure 4.** Observed X-ray, optical (U, F606W, and H bands shown here), and interpolated 6.1 and 22 GHz light curves of the afterglow of GRB 160625B (points), and the best fits given by boxfit (lines) at the indicated participation fraction ξ. The shape of the radio light curve is not well reproduced by any of the fits. Data denoted by gray points are ignored in the fitting (see the text).
Download figure:
Standard image High-resolution image

Table 6. Best-fit Physical Parameters of the Best boxfit Fits to GRB 160625B

Parameter	ξ = 1	ξ = 0.1	ξ = 0.01
p	2.30	2.05	2.05
E_K,iso (erg)	1.8 × 10⁵⁴	1.4 × 10⁵⁴	1.3 × 10⁵⁵
_e	0.13	0.25	0.024
_B	0.030	3.0 × 10⁻⁴	5.8 × 10⁻⁵
n (cm⁻³)	1.1 × 10⁻⁵	0.18	0.96
θ_j (rad)	0.059	0.14	0.13
θ_j (deg)	3.4	7.8	7.2
θ_obs (rad)	0.012	1.1 × 10⁻³	1.1 × 10⁻³
θ_obs (deg)	0.69	0.07	0.06

E_tot (erg)	8.3 × 10⁵¹	4.1 × 10⁵²	1.3 × 10⁵³
η	0.62	0.68	0.19
χ²/dof	8.6	4.6	4.5

Download table as: ASCII Typeset image

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, β, between the optical (r) and X-ray (5 keV) bands, steepens slightly over time, from −0.79 ± 0.02 between 3 and 10 days to −0.86 ± 0.04 at 141 days. This is steeper than −0.65, expected from p ≈ 2.3 implied by the early optical and X-ray light curves (see Section 4.2.1) for ν < ν_c but shallower than −1.15, which is expected for ν > ν_c. Alexander et al. (2017) obtained an early X-ray spectral slope similar to this, ${\beta }_{X}=-{0.86}_{-0.10}^{+0.09}$ , and explained this as ν_c being located just below the X-ray band. However, according to the UKSSDC Swift Burst Analyzer²⁵ , the X-ray photon index Γ_X (and, thus, the spectral slope in X-ray) does not significantly evolve over the first 30 days but stays around ∼1.8, after which the spectrum seems to flatten to Γ_X ∼ 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 Γ_X.²⁶ If ν_c was initially just below X-ray and changed as ν_c ∝ t^−1/2, one would expect the spectrum to instead steepen over time to its ν ≫ ν_c value. We discuss this evolution further in Section 4.2.1.

**Figure 5.** Observed SED of GRB 160625B (points) at late times, interpolated as necessary to the indicated dates, and the best fits given by boxfit (lines) at the indicated participation fraction ξ, using a constant CBM density profile are shown. Data denoted by gray points are ignored in the fitting (see the text).
Download figure:
Standard image High-resolution image

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 we assumed an f_ν ∝ ν^β SED, where β = −0.6 (consistent with ν < ν_c and p ≈ 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 to be 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 β = −0.6. The GRB flux in the $g^{\prime}$ band at 1 day was estimated by subtracting the observed flux at 28 days ( $g^{\prime} =25.39\pm 0.12;$ Laskar et al. 2016) from the flux at 1.0 days ( $g^{\prime} =25.03\pm 0.15;$ Cenko et al. 2016). The host is assumed to dominate at 28 days 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'' aperture) from the flux at 1.2 days (J ≈ 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σ (Cenko et al. 2016; Laskar et al. 2016). Our F110W and F160W observations at 35.3 days made up the other epoch to be fitted simultaneously. The resulting host extinction is A_V = 2.8 ± 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_H,int = 1.52 × 10²² cm⁻² would correspond to A_V ≈ 6.9 mag (Güver & Özel 2009), suggesting a low A_V/N_H ratio by Milky Way standards but higher than that of most GRB hosts. This ratio is consistent with the A_V versus 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 β = −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 ∼0.5 days (4 × 10⁴ 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 ω = 3 is ${\alpha }_{2,{\rm{X}}}=-1.98\pm 0.10$ and for ω = 10, ${\alpha }_{2,{\rm{X}}}=-1.96\pm 0.09$ . The jet-break times, 3.2 ± 0.9 days and 3.7 ± 0.8 days, respectively, are consistent with each other as well.

**Figure 6.** Observed F160W (extinction-corrected), X-ray, and interpolated 6 and 9 GHz light curves of the afterglow of GRB 160509A (points) and our power-law fits including the broken power laws described by Equation (1) (lines). The red triangle is the upper limit of the F160W flux at 24.8 days. X-ray flux densities from Swift/XRT (solid triangles) and Chandra/ACIS-S (open triangles) are reported at 5 keV. Both choices of ω fit the late light curve equally well. The early light curve exhibits a shallower decay and another break, and thus, points before 4 × 10⁴ s (gray) are ignored.
Download figure:
Standard image High-resolution image

Table 7. Parameters of the Best Smooth Broken Power-law Fits to the GRB 160509A X-Ray Light Curve

Parameter	ω = 3	ω = 10
t_j,X	3.2 ± 0.9 days	3.7 ± 0.8 days
α_1,X	−1.06 ± 0.10	−1.20 ± 0.06
α_2,X	−1.98 ± 0.10	−1.96 ± 0.09
Reduced χ²	0.84	0.85

Download table as: ASCII Typeset image

In the radio, we obtained the fluxes at 6 and 9 GHz at the epochs earlier than 79.9 days by power-law interpolation between observed fluxes—our measurements at 36.9 days 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., ≥10 days; Laskar et al. 2016). The resulting decline slopes are ${\alpha }_{6\mathrm{GHz}}=-0.91\pm 0.11$ and ${\alpha }_{9\mathrm{GHz}}=-0.92\pm 0.13$ . Since the reverse shock may still be contributing a nonnegligible fraction of the flux at 10 days, we also performed the fit without this epoch. The results are consistent but less constraining: ${\alpha }_{6\mathrm{GHz}}=-1.07\pm 0.18$ and ${\alpha }_{9\mathrm{GHz}}=-0.92\pm 0.21$ . The slopes at other frequencies between 5 and 16 GHz, fitted from 10 to 20 days, are all consistent with these, ranging from −0.80 ± 0.10 (7.4 GHz) to −1.02 ± 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 days 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,{\rm{F}}160{\rm{W}}}=-2.09\pm 0.10$ , consistent within 1σ with the X-ray decline.

Using boxfit, we again fitted the light curve at three different values of ξ: 1, 0.1, and 0.01. As with the power-law fits, the X-ray points before 0.6 days were ignored, since boxfit cannot accommodate continuous energy injection. Radio points with a significant reverse shock contribution were also ignored (i.e., <10 days; at frequencies <5 GHz also 10.03 days; 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 ∼100 days, 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 spectra 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.

**Figure 7.** Observed X-ray and optical/infrared light curves and the interpolated 6 and 9 GHz light curves of the afterglow of GRB 160509A (points), and the best fit given by boxfit (lines) using a wind-type CBM density profile and ξ = 0.1 are shown. The observed X-ray break is not reproduced (and indeed no break is seen even much later), and therefore, a wind-type CBM is not considered further. Fits using ξ = 1 and ξ = 0.01 produce a similar light curve. Data denoted by gray points are ignored in the fitting (see the text).
Download figure:
Standard image High-resolution image

**Figure 8.** Observed X-ray and optical/infrared light curves and the interpolated 6 and 9 GHz light curves of the afterglow of GRB 160509A (points), and the best fits given by boxfit (lines) at indicated participation fraction ξ, using an ISM-type CBM density profile are shown. The radio light curve shape is again not well reproduced by the fits. Data denoted by gray points are ignored in the fitting (see the text).
Download figure:
Standard image High-resolution image

**Figure 9.** Observed SED of GRB 160509A (points) at late times, interpolated as necessary to the indicated dates, and the best fits given by boxfit (lines) at indicated participation fraction ξ, using an ISM-type CBM density profile are shown.
Download figure:
Standard image High-resolution image

Table 8. Best-fit Physical Parameters of the Best boxfit Fits to GRB 160509A

Parameter	ξ = 1	ξ = 0.1	ξ = 0.01
p	2.29	2.13	2.05
E_K,iso (erg)	8.5 × 10⁵³	3.8 × 10⁵³	3.8 × 10⁵⁵
_e	0.19	0.45	5.7 × 10⁻³
_B	0.015	1.7 × 10⁻⁵	5.8 × 10⁻⁴
n (cm⁻³)	2.1 × 10⁻⁵	18.1	6.1 × 10⁻³
θ_j (rad)	0.046	0.20	0.045
θ_j (deg)	2.6	11.5	2.6
θ_obs (rad)	0.026	0.12	0.027
θ_obs (deg)	1.5	7.0	1.5

E_tot (erg)	1.7 × 10⁵¹	2.5 × 10⁵²	3.9 × 10⁵²
η	0.50	0.69	0.02
χ²/dof	1.8	1.9	1.8

Download table as: ASCII Typeset image

4. Discussion

4.1. The Shape of the Break

In the X-ray, we find little difference in the reduced χ² values of the fits between a sharp and a smooth break for GRB160625B. In the optical, however, fixing ω = 3 results in a visible and significant residual of 4.2σ at 140.2 days, while fixing ω = 10 results in a residual of 1.5σ. The reduced χ² 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 χ² 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 days results in a consistency with the ω = 10 case. We thus conclude that while both values of ω remain plausible, a sharp break with ω = 10 is more likely. A sharp break also implies a small viewing angle θ_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 α₂ ∼ −2.3, and its error is 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 ∼50 days from the burst (∼2.5 × t_j, compared to our latest observations at ∼7 × 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 α₂ = −2.57 ± 0.04, and Strausbaugh et al. (2018) obtained ${\alpha }_{2,\mathrm{optical}}\approx 1.6$ and ${\alpha }_{2,{\rm{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 days, during the "bump" in the light curve between ∼8 days and ∼16 days. Using the same data, Strausbaugh et al. (2018) suggested a break at 12.6 days 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 α₂, leading us to suggest that 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 ν_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 ν_c break in the spectrum and that the optical spectrum eventually becomes consistent with β ∼ −1.1, i.e., the slope above ν_c, which would disfavor a CBM density fluctuation because of this insensitivity. We, however, measure β = −0.86 ± 0.04 between F606W and 5 keV at 141 days, suggesting that ν_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 χ² values of the fits with different ω are close to equal and the post-break slopes are in agreement. A higher θ_obs results in a softer break (Ryan et al. 2015); so in this case, considering that θ_obs ∼ 0.6θ_j (from boxfit), one would expect the break to be softer than for GRB 160625B where θ_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 β = −1.07 ± 0.04, is consistent with p ≈ 2.2 and with ν_c being below the X-ray band (Laskar et al. 2016). As a result, we can use $\alpha =(2\mbox{--}3p)/4$ independent of the CBM distribution (Granot & Sari 2002); in the case of ω = 3, we obtain p = 2.08 ± 0.14 and for ω = 10, p = 2.27 ± 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 ν_c, in the case of a wind-like CBM, α_wind = (1–3p)/4, while for a constant-density CBM similar to the ISM, α_ISM = 3(1−p)/4 (Granot & Sari 2002). Thus, in the optical, one obtains p = 1.63 ± 0.02 in the wind case and p = 2.29 ± 0.02 in the ISM case. Above ν_c, in both cases, α = (2–3p)/4. Comparing the optical and X-ray spectra and fluxes, Alexander et al. (2017) argue that ν_c lies below the X-ray frequencies after ∼1.2 × 10⁴ s, and thus, the early X-ray light curve gives us p = 2.29 ± 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 ∼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 ω = 10 case, the difference between the pre- and post-break power laws is Δα_F606W = −1.00 ± 0.08 in the optical and Δα_X = −0.99 ± 0.16 in the X-ray. Thus, a t^−3/4 factor can be ruled out in the optical at a >3σ level (although in the X-ray, only at a ∼1.5σ level). The difference is larger in the ω = 3 case (>4σ and >2σ, 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 (ω = 3), both the optical and X-ray post-break decline rates are consistent with the form ${f}_{\nu }\propto {t}^{-p}$ , for p ≈ 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 an f_ν ∝ 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⁻¹, steeper than the t^−3/4 expected from the edge effect (Mészáros & Rees 1999; Panaitescu & Mészáros 1999), and the resulting α₂ values are only consistent within 2σ, 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 ± 0.04 (we will denote this as σ_p). They argue this is probably caused by emission from GRB afterglows not behaving exactly as the rather simplified analytical models predict.²⁷ Taking this scatter into account, both α_{2, F606W} and α_2,X in the ω = 10 case (or simply using only the >26.5 days points and a single power law) are in fact consistent within ≈1σ_p with f_ν ∝ 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 ± 0.13 for GRB 160625B, which is consistent with our results in both bands within σ_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)

$\begin{eqnarray}&&{E}_{\mathrm{rot}}\approx 2\times {10}^{52}\,\mathrm{erg}\,\displaystyle \frac{M}{1.4\,{M}_{\odot }}{\left(\displaystyle \frac{R}{10\mathrm{km}}\right)}^{2}{\left(\displaystyle \frac{{P}_{0}}{1\mathrm{ms}}\right)}^{-2}\,,\end{eqnarray} \tag{ 2 }$

where M is the mass, R is the radius, and P₀ is the initial spin period of the newborn magnetar. Metzger et al. (2015) placed a limit of ∼1 × 10⁵³ 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 of the fits from boxfit may technically be achievable with the magnetar model, but with the (more realistic) low ξ values, the required energy approaches or exceeds even this maximum limit. The exceptionally high E_iso can be due to a relatively narrow jet and a lower explosion energy instead, but this requires a high ξ that is inconsistent with simulations by Sironi & Spitkovsky (2011) and Warren et al. (2018)—the best fit at ξ = 1 also results in an extremely low density more typical of intergalactic environments. We do point out a caveat that the parameters of the best fits show a non-monotonic dependence on ξ, 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 2019, 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 well.²⁸ 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 ξ = 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 }_{6.1\mathrm{GHz}}=-1.08\pm 0.11$ and the lack of an observed jet break. At low ξ values, the boxfit fit is somewhat better but only if one ignores the 22.5 days 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 days, which then disappears. Even so, the fit at 22 GHz is worse at all ξ values. At 48 days, the radio SED is consistent with being entirely flat, which is only plausible in the standard model around a very smooth ν_m break. While the low ξ fits do place the ν_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 10 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 (ν_m < ν_c) are −p or −1/3, depending on which side of ν_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,6.1\mathrm{GHz}}\gtrsim 10\times {t}_{j,{\rm{F}}606{\rm{W}}}$ . The possibilities given by the standard jet model that are consistent with the slope are:

1.
Post-break, ν_c < ν_m, i.e., fast-cooling: α_6.1GHz is consistent with the expected decline of α₂ = −1. However, the measured ${\alpha }_{1,{\rm{F}}606{\rm{W}}}=-0.96\pm 0.01$ does not match the pre-break decline expected at any frequency in this scenario.
2.
Pre-break, ${\nu }_{m}\lt 6.1\,\mathrm{GHz}\lt {\nu }_{c}$ : ${\alpha }_{6.1\mathrm{GHz}}$ is consistent with p = 2.4 and α = 3(1−p)/4 = −1.05 (Granot & Sari 2002). However, the spectral index between radio and optical is −0.35 ± 0.03 at 22 days and −0.49 ± 0.01 at 140 days, which is intermediate between the indices expected above and below ν_m (respectively, (1−p)/2 ≈ −0.65 and 1/3) and, thus, implies that ν_m > 6.1 GHz at 140 days, or that p ≈ 2.0.
3.
A transition to a nonrelativistic flow, ν_m < 6.1 GHz < ν_c: the expected slope is (21–15p)/10 (van der Horst 2007), resulting in p = 2.12 ± 0.08, which is consistent with our estimate within σ_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 ξ = 1 and all radio fluxes divided by 10 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 GHz > ν_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 ν_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 (θ_obs < θ_j,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 ν > ν_c by ∼50 days. This would indicate a unique cooling behavior inconsistent with the standard expectations. The observed lack of evolution of the Swift spectra until 30 days implies that the X-ray spectral slope β_X is not the result of a ν_c break right below the X-ray frequencies, as this would require the spectrum to soften over time to its slope at ν ≫ ν_c. Furthermore, the optical-to-X-ray index is observed to gradually steepen and eventually become similar to β_X. This is qualitatively consistent with the reddening of the optical spectrum noted by Strausbaugh et al. (2018). In addition, β_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 ν_c; this can result in β = 3(1−p)/4, which would imply p ≈ 2.1. This harder spectrum is expected to dominate when the _e/_B ratio is high, which would fit the low-ξ boxfit results.

4.2.2. GRB 160509A

In the case of GRB 160509A, the change in X-ray decay slope across the break, Δα_X = −0.75 ± 0.11 for a sharp break and Δα_X = −0.92 ± 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σ or 3σ level, depending on ω. However, when considering the intrinsic p scatter of σ_p = 0.25 (Gompertz et al. 2018), ${\alpha }_{2,{\rm{X}}}$ is also consistent with an f_ν ∝ 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,{\rm{F}}160{\rm{W}}}\,=2.09\pm 0.10$ is marginally consistent (1.1σ) with p ≈ 2.2, but a lack of pre-break data prevents us from determining ${\rm{\Delta }}{\alpha }_{{\rm{F}}160{\rm{W}}}$ .

The decline of the afterglow in the radio after 10 days is about f ∝ 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 ν_m (Rhoads 1999). As with GRB 160625B, we list the possibilities consistent with this decline, allowed by standard jet theory:

1.
Post-break, ν_c < ν_m, i.e., fast-cooling: α = −1 is expected and consistent with α_radio, but this scenario is incompatible with the measured IR-to-X-ray spectral index −0.74 ± 0.09 at 35 days, as the expected index is −p/2 ≈ −1.1 (a photon index consistent with this is indeed seen in the X-ray at earlier times according to the UKSSDC Swift Burst Analyzer²⁹ — <β_X > = 1.06 ± 0.04 between 1 and 10 days, indicating that ν_c is still between X-ray and optical frequencies and ν_c > ν_m at 35 days).
2.
Pre-break, ν_m < 6 GHz < ν_c, ISM-like CBM: α_radio is consistent with the expected decline (α = 3(1−p)/4 = −0.9 assuming p = 2.2; Granot & Sari 2002), but the observed spectral index of −0.40 ± 0.01 between F160W and 9 GHz at 35 days implies ν_m > 9 GHz.
3.
A transition to a nonrelativistic flow, ν_m < 6 GHz < 9 GHz < ν_c: the expected slope is (21–15p)/10, resulting in p = 2.01 ± 0.08—again, consistent with our estimate within σ_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 ξ = 1 places a smooth and, thus, off-axis, jet-break at a later time, around 35 days 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 ∼ 3.5 days, though, and at lower, more realistic values of ξ 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,{\rm{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 (∼20 days), 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 it 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 ξ. 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 8 ³⁰ —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 × 10⁵¹ erg at ξ = 1 is well below the maximum rotational energy of a millisecond magnetar (see Section 4.2.1) Once again, we deem the lower ξ values to be more realistic based on simulations (Sironi & Spitkovsky 2011; Warren et al. 2018) and the fact that ξ = 1 results in an extremely low density. Low ξ values require energies around ∼3 × 10⁵² 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 η of the prompt γ-ray emission depends on the value of ξ used in the fitting but not monotonically: with ξ = 0.01, one obtains a much lower value for η than otherwise. In both cases, the difference in χ² between the ξ = 0.1 and ξ = 0.01 fits is minimal, and in the case of GRB 160509A, so is the difference between ξ = 1 and ξ = 0.01; thus, we cannot reliably distinguish between these scenarios. In the literature, it is commonly assumed that ξ = 1, and high values of η are obtained: for example, Lloyd-Ronning & Zhang (2004) find values as high as η ∼ 1 depending on E_iso, and mostly η ≳ 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 η ≲ 0.1 (Kumar & Zhang 2015, and references therein). Our results may indicate that if very low values of ξ 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 GRB 160509A, 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 nonthermal distribution) is fixed at ξ = 0.1 or ξ = 0.01, requiring energies of ∼3 × 10⁵² or even ∼10⁵³ 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 multicomponent 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 α = −1.19 ± 0.06 until ∼60 days) 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.C. acknowledges support from the National Science Foundation via CAREER award #1455090. A.d.U.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.C. 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.

Facilities: HST(WFC3) - Hubble Space Telescope satellite, 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 & Lupton 1998; Alard 2000); boxfit (van Eerten et al. 2012).