The outflow structure of GW170817 from late time broadband observations

We present our broadband study of GW170817 from radio to hard X-rays, including NuSTAR and Chandra observations up to 165 days after the merger, and a multi-messenger analysis including LIGO constraints. The data are compared with predictions from a wide range of models, providing the first detailed comparison between non-trivial cocoon and jet models. Homogeneous and power-law shaped jets, as well as simple cocoon models are ruled out by the data, while both a Gaussian shaped jet and a cocoon with energy injection can describe the current dataset for a reasonable range of physical parameters, consistent with the typical values derived from short GRB afterglows. We propose that these models can be unambiguously discriminated by future observations measuring the post-peak behaviour, with slope -1.0 for the cocoon and -2.5 for the jet model.


I N T RO D U C T I O N
The discovery of GW170817 and its electromagnetic counterparts (GRB170817A and AT2017gfo; Abbott et al. 2017a;Coulter et al. 2017;Goldstein et al. 2017) ushered in a new era of multimessenger astrophysics, in which both gravitational waves and photons provide complementary views of the same source. While observations at optical and infrared wavelengths unveiled the onset and evolution of a radioactive-powered transient, known as kilonova, observations at X-rays and, later, radio wavelengths probed a different component of emission, likely originated by a relativistic outflow launched by the merger remnant. Troja et al. (2017) explained the observed X-ray and radio data as the onset of a standard short GRB (sGRB) afterglow viewed at an angle (off-axis). However, as already noted in Troja et al. (2017) and Kasliwal et al. (2017), a standard top-hat jet model could explain the afterglow data set collected at early times, but failed to account for the observed gamma-ray emission. Based on this evidence, Troja et al. (2017) suggested that a structured jet model (e.g. Zhang et al. 2004;Kathirgamaraju, Barniol Duran & Giannios 2018) provided a coherent description of the entire broad-band data set. Within this framework, the peculiar properties E-mail: eleonora@umd.edu; eleonora.troja@nasa.gov of GRB170817A/AT2017gfo could be explained, at least in part, by its viewing angle (see also Lamb & Kobayashi 2017;Lazzati et al. 2017). An alternative set of models invoked the ejection of a mildly relativistic wide-angle outflow, either a jet-less fireball (Salafia, Ghisellini & Ghirlanda 2018) or a cocoon (Nagakura et al. 2014;Hallinan et al. 2017). In the latter scenario, the jet might be chocked by the merger ejecta (Mooley et al. 2017), and the observed gamma-rays and broad-band afterglow emission are produced by the expanding cocoon. The cocoon may be energized throughout its expansion by continuous energy injection. In this paper detailed models of structured jet and cocoon, from its simplest to more elaborate version, are compared with the latest radio to X-ray data. Predictions on the late-time evolution are derived, and an unambiguous measurement capable of disentangling the outflow geometry, jet versus cocoon, is presented.

X-rays
The Chandra X-ray Observatory and the Nuclear Spectroscopic Telescope ARray (NuSTAR) re-observed the field of GW170817 soon after the target came out of sunblock. Chandra data were 9 GHz a Units are 10 −15 erg cm −2 s −1 for X-ray fluxes, and μJy for radio fluxes. b NuSTAR 3σ upper limit.

Figure 1.
Exposure corrected X-ray images of the field of GW 170817 taken at 9 d (left; Troja et al. 2017) and 108 d (right) after the merger. The X-ray counterpart, marked by the crossed lines, significantly brightened between the two epochs.
reduced and analysed in a standard fashion using CIAO v. 4.9 and CALDB 4.7.6. The NuSTAR data were reduced using standard settings of the pipeline within the latest version of NuSTAR Data Analysis Software. Spectral fits were performed with XSPEC by minimizing the Cash statistics. A log of observations and their results is reported in Table 1. We found that the X-ray flux at 108 d is ≈5 times brighter than earlier measurements taken in August 2017 (Fig. 1), thus confirming our previous findings of a slowly rising X-ray afterglow (Troja et al. 2017). By describing the temporal behaviour with a simple power-law function, we derive an index α ≈ 0.8, consistent with the constraints from radio observations (Mooley et al. 2017). During each observations, no significant temporal variability is detected on time-scales of ≈1 d or shorter.
In order to probe any possible spectral evolution we computed the hardness ratio (Park et al. 2006), defined as HR = H−S/H+S, where H are the net source counts in the hard band (2.0-7.0 keV) and S are the net source counts in the soft band (0.5-2.0 keV). This revealed a possible softening of the spectrum, with HR ≈ −0.1 for the earlier observations ( 15 d) and HR ≈ −0.5 for the latest observations (>100 d). However, the large statistical uncertainties prevent any firm conclusion.

Radio
The target GW170817 was observed with the Australia Telescope Compact Array in three further epochs after 2017 September. Observations were carried out at the centre frequencies of 5.5 and 9 GHz with a bandwidth of 2 GHz. For these runs the bandpass calibrator was 0823-500, the flux density absolute scale was determined using 1934-638 and the source 1245-197 was used as phase calibrator. Standard MIRIAD procedures were used for loading, inspecting, flagging, calibrating, and imaging the data. The target was clearly detected at all epochs, our measurements are reported in Table 1.
The broad-band spectrum, from radio to optical (Lyman et al. 2018) to hard X-rays ( Fig. 2), can be fit with a simple power-law model with spectral index β = 0.575 ± 0.010 and no intrinsic absorption. A fit with a realistic afterglow spectrum (Granot & Sari 2002) constrains the cooling break to ν c 1 keV (90 per cent confidence level).

Jet and cocoon
The standard model for sGRB afterglows describes these in terms of synchrotron emission from a decelerating and decollimating relativistic jet. More recent studies have argued for the additional presence of a slower moving cocoon (e.g. Nagakura et al. 2014;Lazzati et al. 2017;Mooley et al. 2017) also in the case of sGRBs. Numerical studies of jet breakouts have revealed a range of possibilities for jet velocities and initial angular structure. In the case of jetted outflow seen from a substantial angle θ ν (i.e. larger than the opening angle θ c of a top-hat flow, or in the wings of a jet with energy dropping as a function of angle), relativistic beaming effects will delay the observed rise time of the jet emission.
The dynamics of the jet component can be treated semianalytically at various levels of detail (e.g. Rossi & Rossi 2006), depending on additional assumptions for the angular structure of the jet. Early X-ray and radio observations already rule out (Troja et al. 2017) the universal jet structure (i.e. a power-law drop in energy at larger angles), and we will not discuss this option further here. For the Gaussian structured jet, we assume energy drops according to E(θ) = E 0 exp[−θ 2 /2θ 2 c ], up to a truncating angle θ w . We approximate the radial jet structure by a thin homogeneous shell behind the shock front. Top hat and Gaussian jet spreading are approximated following the semianalytical model from Van Eerten, Zhang & MacFadyen (2010), with the Gaussian jet implemented as a series of concentric top-hat jets. This spreading approximation was tuned to simulation output (Van Eerten, Zhang & MacFadyen 2010) that starts from top-hat initial conditions but develops a more complex angular structure overtime. Since the off-axis jet emission will fully come into view after deceleration, deceleration radius, and initial Lorentz factor value no longer impact the emission and are not included in the model.
The cocoon is similarly treated using a decelerating shell model, now assuming sphericity and including pre-deceleration stage. For a direct comparison to Mooley et al. (2017), we have added a mass profile that accounts for velocity stratification in the ejecta and thereby provides for ongoing energy injection. The total amount of energy in the slower ejecta above a particular four-velocity is assumed to be a power law E >u (u) = E inj u −k for u ∈ [u min , u max ] (note that for relativistic flow u → γ ). The energy from a slower shell is added to the forward shock once this reaches the same velocity. The total cocoon energy (and therefore light curve turnover time) is thus dictated by u min . We assume an initial cocoon mass of M ej . Both jet and cocoon emerge into a homogeneous environment with number density n.
We estimate the emission of the ejecta using a synchrotron model (Sari, Piran & Narayan 1998). Electrons are assumed to be accelerated to a power-law distribution in energy of slope −p, containing a fraction e of the post-shock internal energy. A further fraction B is estimated to reside in shock-generated magnetic field energy. We integrate over emission angles to compute the observed flux for an observer at luminosity distance d L and redshift z.
The key distinguishing features of the various models are their rise slope, peak time, and decay slope. Above synchrotron injection break ν m , a cocoon will show a steeply rising flux F ν ∝ t ∼3 . A top-hat jet will show a steeper rise while a Gaussian energy profile will show a more gradual rise than a top-hat jet (see Extended fig. 3 of Troja et al. 2017). Gaussian and top-hat jets will have peak times following (Troja et al. 2017): while a cocoon outflow whose energy is dominated by the slow ejecta will peak at a time according to t peak ≈ 81 kE 50 n −3 u 8 min 1/3 d. (cocoon) Here E 0,50 is the on-axis equivalent isotropic energy in units to 10 50 erg, E 50 the total cocoon energy following energy injection in the same units, and n −3 circumburst density in units of 10 −3 cm −3 . Cocoon and jet models differ in their expected post-peak downturn slopes. For the cocoon, ∼t −1.0 (a decelerating spherical fireball) is expected, while for jet models ∼t −2.5 is expected due to a combination of jet spreading dynamics and the entire jet having come into view (Fig. 3). A jet plus cocoon might yield an intermediate slope value, while still confirming the collimated nature of the sGRB.

Bayesian model fit including LIGO/VIRGO constraints
We perform a Bayesian Markov-Chain Monte Carlo (MCMC) model fit to data from synthetic detections generated from our semianalytical cocoon and jet models. At fixed distance, the off-axis structured jet models considered here are fully determined by the set of eight parameters jet = {θ v , E 0 , θ c , θ w , n, p, e , B }. The isotropic cocoon with a power-law velocity distribution, on the other hand, requires nine parameters cocoon = {u max , u min , E inj , k, M ej , n, p, e , B }. We generate samples of the posterior for both models using the affine-invariant ensemble MCMC sampler implemented in the EMCEE package (Goodman & Weare 2010;Foreman-Mackey et al. 2013). For both the cocoon and jet models we initialize the MCMC walkers in a small ball near the maximum of the posterior, calculated through trial runs. We run each model with 300 walkers for 128 000 steps, dropping the first 36 000 steps as an initial burn-in phase, generating ∼3 × 10 7 posterior samples.
We assign independent priors for each parameter, uniform for θ c , θ w , k, and p, and log-uniform for E 0 , u max , u min , E inj , M ej , n, e , and B . The viewing angle θ v is given a prior p(θ v ) ∝ sin θ v . This is proportional to the solid angle subtended at this viewing angle and is the expected measured distribution of randomly oriented sources in the sky if observational biases resulting from jet dynamics and beaming are not accounted for. Parameters are given wide bounds so as to ensure the parameter space is fully explored. The era of multimessenger astronomy allows us to directly link the observational constraints from the different channels. The upper bound on the viewing angle is therefore chosen to include the 95 per cent confidence interval from the LIGO analysis of GW170817A assuming either value (CMB or SNe) of the Hubble constant (Abbott et al. 2017b).
The LIGO/VIRGO analysis of GW170817A includes the inclination i, the angle between the total angular momentum vector of the binary neutron star system and the line of sight (Abbott et al. 2017a, b). Assuming this angle to be identical to the viewing angle θ v , we can incorporate their posterior distributions p GW (θ v ) into our own analysis of the Gaussian structured jet. We do this by applying MNRASL 478, L18-L23 (2018) Downloaded from https://academic.oup.com/mnrasl/article-abstract/478/1/L18/4969981 by California Institute of Technology user on 04 June 2018 a weighting factor p GW (θ v )/p(θ v ), where p(θ v ) is our prior, to the MCMC samples. This is valid so long as the MCMC adequately samples the region of parameter space favoured by p GW , which we confirmed for our analysis.
LIGO/VIRGO report three distributions p GW , one using only the gravitational wave data and two incorporating the known redshift of host galaxy NGC 4993, utilizing the value of the Hubble constant reported by either the Planck collaboration or the SHoES collaboration (Planck Collaboration XIII 2016; Riess et al. 2016;Abbott et al. 2017b). We incorporate the latter two distributions into our analysis. The reported distribution functions ( fig. 3 of Abbott et al. 2017a, b) were digitized and found to be very well described by Gaussian distributions in cos θ v . We take the distributions to be where for the Planck H 0 μ 0 = 0.985 and σ = 0.070 and for the SHoES H 0 μ 0 = 0.909 and σ = 0.068. The overall normalization of p GW need not be specified in this approach.

Constraints to jet and cocoon models
The slow rise ∝ t 0.8 observed in the radio and X-ray data is inconsistent with the steep rise predicted by the basic cocoon model without energy injection (Hallinan et al. 2017). The top-hat jet model is also characterized by a steep rise and a narrow plateau around the peak, as the result of the jet energy being limited to a narrow core. This is a well defined property of this model, rather independent of the other afterglow parameters, and the late-time observations of GW170817 allow us to robustly reject it. Some degree of structure in the GRB outflow is therefore implied by the late-time observations. The universal jet model is too broad and overpredicts the early time X-ray flux (Troja et al. 2017). As discussed in the next section, our models of Gaussian jet and cocoon with energy injection provide equally good fits to the current data set, and are basically indistinguishable until the peak time. However, as discussed in Section 3, the post-peak slopes are expected to differ, with F ν ∝ t ∼−1.0 for the cocoon and F ν ∝ t ∼−2.5 , and intermediate values indicating a combination of directed outflow plus cocoon. Future observations will be critical to distinguish the nature of the outflow (collimated versus isotropic).

Bayesian analysis results
The MCMC fit results for the Gaussian jet and cocoon models are summarized in Table 2. Both models are consistent with current data and have similar quality of fit with χ 2 per degree of freedom near unity. Corner plots showing the 1D marginalized posterior distributions for each parameter and the 2D marginalized posterior distributions for each pair of parameters are included in the on-line materials (Figs 5 and 6). The Gaussian jet model prefers a narrow core (θ c ∼ 5 deg) with wings truncated at several times the width of the core. The viewing angle is significant (θ v ∼ 30 deg) but degenerate with θ c , E 0 , and n 0 . The large uncertainties on individual parameters are partially a result of these degeneracies within the model. The viewing angle correlates strongly with θ c and n 0 and anticorrelates with E 0 . Fig. 4 (left-hand panel) shows the range of possible X-ray and radio light curves of the Gaussian structured jet, pulled from the top 68 per cent of the posterior. There is still significant freedom within the model, which will be better constrained once the emission peaks.
Incorporating constraints on θ v from LIGO/VIRGO tightens the constraints on θ c , E 0 , and n 0 due to the correlations between these parameters. Using either the Planck or SHoES H 0 decreases the likely θ v and θ c significantly, with comparatively small adjustments to E 0 and n 0 .
The cocoon model strongly favours an outflow of initial maximum four-velocity u ∈ [3.7, 18.6] whose energy is dominated by a distribution of slower ejecta. The mass of the fast ejecta M ej and the low-velocity cut-off u min are very poorly constrained. The total energy of the outflow is strongly dependent on u min , which can only be determined by observing the time at which the emission peaks. Fig. 4 (right-hand panel) shows the range of possible X-ray and radio light curves of the cocoon model.
The posterior distributions of e and B under both models are consistent with theoretical expectations, with the cocoon model providing weaker constraints than the jet. Both models very tightly constrain p = 2.156 ± 0.015, a consequence of simultaneous radio and X-ray observations, and place the cooling frequency above the X-ray band. The cocoon model tends to prefer a larger total energy and smaller circumburst density than the Gaussian jet, although the posterior distributions of both quantities are broad.

Implications for the prompt γ -ray emission
In the case of a Gaussian jet, it is natural to consider the question whether the sGRB and its afterglow would have been classified as typical, had the event been observed on-axis. The afterglow values for E 0 , θ c , and θ w are indeed consistent with this notion. As discussed previously by Troja et al. (2017), we expect the observed γ -ray isotropic equivalent energy release E γ ,obs ∼ 5 × 10 46 erg to scale up to a typical value of E γ ,OA ∼ 2 × 10 51 erg, once the orientation of the jet is accounted for. This implies a ratio θ v /θ c = 2 ln[E γ,OA /E γ,obs ] ≈ 4.6. From our afterglow analysis we infer a value of 5.6 ± 0.9 (95 per cent uncertainty) for this ratio, accounting for the correlation between the two angles shown by our fit results. Our inferred range of θ v / θ c lies marginally above the typical value, but remains consistent with expected range of sGRB energetics.
While the structured jet model implies that GRB 170817A would have been observed to be a typical sGRB when seen on-axis, the cocoon model implies that it belongs to a new class of underluminous gamma-ray transients. In the former case, the origin of standard sGRBs as due to neutron star mergers can be considered confirmed, with the added benefit of being able to combine multimessenger information about jet orientation into a single comprehensive model fit to the data (Table 2). In the latter case, the future detection rate of multimessenger events is potentially higher and dominated by these failed sGRBs. Continued monitoring at X-ray, optical, and radio wavelengths is important to discriminate between these two different scenarios.

C O N C L U S I O N S
Our modelling of the latest broad-band data confirms that a jetted outflow seen off-axis is consistent with the data (Troja et al. 2017). The late-time data favour a Gaussian shaped jet profile, while homogeneous and power-law jets are ruled out. A simple spherical cocoon model also fails to reproduce the observed behaviour and, to be successful, a cocoon with energy injection from earlier shells catching up with the shock is required. Both models can describe the data for a reasonable range of physical parameters, within the observed   range of sGRB afterglows. A Gaussian jet and a re-energized cocoon are presently indistinguishable but we predict a different behaviour in their post-break evolution once the broad-band signal begins to decay, with F ν ∝ t ∼−1.0 for the cocoon and F ν ∝ t ∼−2.5 , and intermediate values indicating a combination of directed outflow plus cocoon. While the cocoon model invoke a new class of previously unobserved phenomena, the Gaussian jet provides a self-consistent model for both the afterglow and the prompt emission and explains the observed properties of GW170817 with a rather normal sGRB seen off-axis.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at MNRASL online.