A thousand days after the merger: continued X-ray emission from GW170817

Recent observations with the Chandra X-ray telescope continue to detect X-ray emission from the transient GW170817. In a total exposure of 96.6 ks, performed between March 9 and March 16 2020 (935 d to 942 d after the merger), a total of 8 photons are measured at the source position, corresponding to a significance of about 5 sigma. Radio monitoring with the Australian Telescope Compact Array (ATCA) shows instead that the source has faded below our detection threshold (<33 uJy, 3 sigma). By assuming a constant spectral index beta=0.585, we derive an unabsorbed X-ray flux of approximately 1.4E-15 erg/cm^2/s, higher than earlier predictions, yet still consistent with a simple structured jet model. We discuss possible scenarios that could account for prolonged emission in X-rays. The current dataset appears consistent both with energy injection by a long-lived central engine and with the onset of a kilonova afterglow, arising from the interaction of the sub-relativistic merger ejecta with the surrounding medium. Long-term monitoring of this source will be essential to test these different models.


INTRODUCTION
On August 17th 2017, advanced LIGO and Virgo observed the first gravitational wave signal from a binary neutron star (NS) merger (Abbott, et al. 2017a). This event, named GW170817, was followed by a short duration gamma-ray burst, GRB170817A, and, 9 days later, by a non-thermal afterglow emission, visible across the electromagnetic spectrum (Abbott, et al. 2017b;Troja, et al. 2017;Hallinan, et al. 2017). After an initial rising phase, F ∝ t 0.8 (Troja, et al. 2018;Mooley, et al. 2018a;Lyman, et al. 2018;Margutti, et al. 2018;Ruan, et al. 2018), the afterglow peaked at ≈160 d after the merger and then started a rapid decay phase, E-mail: eleonora.troja@nasa.gov F ∝ t −2.2 (Mooley, et al. 2018b;Lamb, et al. 2019;). This behavior is markedly different from the garden-variety GRB afterglows, observed to fade within a few minutes since the burst.
The low-luminosity of the gamma-ray emission and the atypical temporal evolution of the afterglow component are widely interpreted as manifestation of a highly-relativistic structured jet seen at an angle of ≈20-30 deg from its axis (Troja, et al. 2017;Abbott, et al. 2017b;Lazzati, et al. 2018;Lyman, et al. 2018;Troja, et al. 2018;Mooley, et al. 2018b;Margutti, et al. 2018;Lamb, et al. 2019;Ryan, et al. 2019;Hajela, et al. 2019). In this model, the energy and Lorentz factor of the relativistic ejecta vary with the angle from the jet's axis (e.g. Zhang & Mészáros 2002). The initial rising slope and the peak time strongly depend on the observer's viewing angle and the jet's angular profile (Ryan, et al. 2019). However, the post-peak behavior is dominated by the emission from the jet's core and should resemble the post jet-break evolution of a standard GRB afterglow. Even in this case, the post-break evolution can exhibit a rich behaviour, and is sensitive to the nature of the spreading dynamics of the decelerating relativistic plasma and to gradients in the circumburst ambient gas mass distribution. At sufficiently late times, emission from the jet as it has decelerated to non-relativistic flow velocities will begin to dominate the total observed flux, leading to a change in slope relative to the relativistic limit (Frail, Waxman & Kulkarni 2000). If a counter jet was launched, this too will at some point become visible (van Eerten, Zhang & Mac-Fadyen 2010). However, very few GRBs are close enough to remain continuously visible for years and, for this reason, the jet's late-time evolution is rarely probed by observations at wavelengths other than radio (e.g. De Pasquale, et al. 2016;Kouveliotou, et al. 2004).
Changes in the light curve evolution can also be the product of a genuinely new feature of the outflow not previously detected. Of particular interest to the case of neutron star mergers are scenarios that relate directly to the nature of the remnant (such as prolonged energy injection from a long-lived central engine,  and to the sub-relativistic merger ejecta, producing a low-luminosity late-peaking afterglow (Nakar & Piran 2011;Hotokezaka, et al. 2018;. In the case of GW170817, evidence suggests that a substantial amount ( 0.01 M ) of fast ( 0.1 c) ejecta comes from the luminous kilonova emission AT2017gfo (Arcavi, et al. 2017;Evans et al. 2017;Drout et al. 2017;Kasen, et al. 2017;Kasliwal et al. 2017;Nicholl, et al. 2017;Pian et al. 2017;Shappee, et al. 2017;Smartt et al. 2017;Tanvir et al. 2017;Troja, et al. 2017). As these ejecta continue to expand they will drive a blastwave in the local medium, begin decelerating as more mass is swept up, and emit synchrotron radiation from the blast waveâȂŹs forward shock. This emission, which we refer to as kilonova afterglow, peaks years after the initial burst and, at the distance of GW170817, may be bright enough to be detected with current instruments.
In order to explore the late-time behavior of the relativistic jet and constrain alternative components of emission, the location of GW170817 is periodically monitored at radio and X-ray energies. In this work, we present the results of the long-term monitoring campaign with the Chandra X-ray observatory and the Australian Telescope Compact Array (ATCA), and discuss the possible origins of the observed long-lived X-ray emission. Throughout this paper, we adopt a distance of 40 Mpc and a standard ΛCDM cosmology (Planck Collaboration, et al. 2018). Unless otherwise stated, the quoted errors are at the 68% confidence level, and upper limits are at the 3 σ confidence level.

X-rays
We presented the analysis of the first year of observations in . Since then, the target GW170817 is being monitored by the Chandra X-ray Telescope with a cadence of approximately six months under Guest Ob- where H and S are the net source counts in the hard (2.0-7.0 keV) and soft (0.5-2.0 keV) energy bands, respectively. Error bars represent 1 σ uncertainties. The last three epochs (gray symbols) were binned into a single point in order to improve the signalto-noise ratio. Horizontal lines show the values expected for an absorbed power-law with photon index Γ=2.0 (dotted line), 1.5 (dot-dashed line), and 1.25 (dashed line).
server programs 20500691 (PI: Troja) and 20500299 (PI: Margutti). These three additional epochs (Table 1) track the afterglow evolution from 1.6 to 2.6 years after the merger. The temporal evolution of the X-ray counterpart is shown in Fig. 1. Each epoch was split into multiple observations. Each observation was reduced in a standard fashion using the  In order to correct for small positional errors between different observations, we used the tool reproject aspect to determine a new aspect solution based on common bright point sources. Each observation was reprocessed using the updated astrometric information. Data were filtered with the task deflare to remove background flares by applying a sigma clipping threshold of 3. Observations carried out at a similar epoch were merged into a single image using the task flux obs. The resulting total exposures are 98.8 ks (Epoch 1), 98.9 ks (Epoch 2), and 96.6 ks (Epoch 3). Aperture photometry was performed in the broad 0.5-7.0 keV energy band. Source counts were extracted from the merged images using a circular aperture containing 92% of the encircled energy fraction, whereas the background contribution was estimated from nearby source-free regions. Xray emission from the position of GW170817 is visible at all epochs. We estimated the detection significance following the Bayesian method of Kraft, Burrows & Nousek (1991), and report in Table 1 the equivalent value for a normal probability distribution.
Due to the low number of counts, the source spectral properties can not be adequately constrained. In order to check for possible spectral evolution we computed the hardness ratio (HR; Park, et al. 2006), defined as the ratio (H −S)/(H +S), where H and S are the net source counts in the hard (2.0-7.0 keV) and soft (0.5-2.0 keV) energy bands, respectively. The HR light curve (Fig. 2) shows a possible hardening of the spectrum at late times (t 1.5 yr), although with low significance.
X-ray fluxes were calculated assuming an absorbed power-law spectrum with column density fixed to the Galactic value 1.1×10 21 cm −2 (Willingale, et al. 2013) and a photon index Γ= β +1 = 1.585, where β is the spectral index derived from broadband afterglow modeling . A harder spectrum would increase our flux estimate by ≈13% (for Γ=1.25), still within the statistical uncertainties of the measurement. Our values are lower, yet consistent within the large uncertainties, than those reported in Hajela, et al. (2019). Our conversion into fluxes is based on the broadband (from radio to X-rays) spectral shape and does not change over time, whereas Hajela, et al. (2019) derives variable conversion factors based on single-epoch X-ray observations. The latter approach is subject to greater uncertainty, and does not take into account the full spectral information available from the multi-wavelength dataset.
1 An independent analysis of this data set reports a similar countrate and a 50% higher X-ray flux (Hajela, et al. 2020). We can reproduce this result only by assuming a hard spectrum with

Radio
We re-observed the position of GW170817 with ATCA (program C3240; PI:Piro) on May 3rd, 2020 (990 d since the merger) for 11 hours. The array configuration was 6A, the centre observing frequency was 2.1 GHz and the observing bandwidth was 2 GHz. The usual primary calibrator 1934-638 was not observed, instead the band-pass calibrator 0823-500 was used to bootstrap the absolute flux density scale assuming a flux density of 6.38 Jy and a spectral slope of -0.215. The source 1245-197 was used as the phase calibrator. The data set was calibrated and imaged in Miriad using standard procedures. The array configuration resulted in a E-W angular resolution of 6.5 arcsec, sufficient to separate the target from its host galaxy NGC 4493.
No detection was found at the position of GW170817 in the natural-weighted restored image. A 3 σ upper limit of 33 µJy was estimated from rms noise statistics in a region of the restored image away from bright radio sources. This measurement constrains the broadband spectral index to β <0.68.

MODEL FITTING METHODS
Throughout this paper we continue our practice from Troja, et al. (2018 ;; Ryan, et al. (2019) of performing Bayesian fits using the model and afterglowpy software 2 described in Ryan, et al. (2019). This approach combines a decelerating spreading shell model (van Eerten, Zhang & MacFadyen 2010) that includes a range of options for lateral and radial energy structure with the emcee (version 2.2.1) Python package for Markov-Chain Monte Carlo analysis (Foreman-Mackey et al. 2013). For the jet model with a Gaussian distribution of lateral energy, the parameters are: fraction of post-shock internal energy in magnetic field εB, fraction of post-shock internal energy in the accelerated electron population εe, power-law slope of the electron population −p, homogeneous circumburst medium number density n0, on-axis isotropic equivalent energy E0, jet orientation θv, jet core width θc, and jet total width θw. We also perform fits that include an additional constant X-ray component, specified by a flux density FX . This accounts for additional sources of emission, such as a long-lived engine or a separate source at close proximity on the sky. We use the same prior on jet orientation as reported in earlier work (Troja, et al. 2018), drawn from Abbott, et al. (2017c) with Γ=0.57, drastically different from the spectral properties of the GW afterglow. 2 https://github.com/geoffryan/afterglowpy a Hubble constant as determined by Planck Collaboration, et al. (2018). The additional component FX is given a flat prior and bounded by 0 < FX < 2 × 10 −4 µJy.
In order to explore the non-thermal emission from the sub-relativistic ejecta, we consider a quasi-spherical "kilonova afterglow" model. While the bulk of the kilonova material coasts at a sub-relativistic velocity, it is expected a less massive tail of material outflows with substantially higher velocities (Bauswein, Goriely & Janka 2013;Hotokezaka, et al. 2013). The material is postulated to have an energy distribution which is a power-law in the four-velocity: E>u(u) = Etot(u/umin) −k . We use the same MCMC routines as with the structured jet analysis and the isotropic outflow model from Troja, et al. (2018), reparameterized for a kilonova-like outflow. This model is specified by a power-law k stratification of ejecta velocities, a total ejecta mass Mej = 2k/(k + 2)u −2 min Etotc −2 , a maximum ejecta fourvelocity umax, a minimum velocity βmin, as well as the environmental and synchrotron parameters n0, p, εe, and εB. It is not a given that εe, εB and p are identical for jet and kilonova component.
The structured jet fits used a parallel tempered ensemble MCMC sampler with 20 geometrically spaced temperatures between 1 and 10 6 . Each temperature rung was occupied by 100 walkers, and the chain was run for 20,000 iterations. The kilonova afterglow fits were run using a standard ensemble sampler with 300 walkers for 64,000 iterations. Further details of the method can be found in the references listed above. Our models were compared to the X-ray, radio and optical afterglow light curves using the same data set described in (Troja, et al.  To compare different models we utilize the Widely Applicable Information Criterion (WAIC; Watanabe 2010). The WAIC is an estimate of the "expected log predictive density" (elpd ): a score measuring the likelihood new data will be well described by the current model (Gelman, Hwang, & Vehtari 2013). The elpd measures the predictive power of a fit, it rewards a tight match to the data while penalizing over fitting and extraneous parameters. The WAIC is proven to be asymptotically equal to the elpd for a wide range of models and is straightforward to compute from MCMC posterior samples, whereas the elpd itself can only be computed if the true model is known. We use the pWAIC2 estimator for the effective number of parameters (Gelman, Hwang, & Vehtari 2013).
Following Vehtari, Gelman, &, Gabry (2017) we compute the WAIC score for each model at every data point. The total WAIC score WAIC elpd for a model is the sum of the scores for each data point. Each model score and score difference ∆WAIC elpd have a standard error computed from the variance over the contributions from each data point. This standard error is likely optimistic but within a factor of 2 of the true value (Bengio & Grandvalet 2004). In a twoway comparison, a model is favoured if its ∆WAIC is several times larger than its standard error.

RESULTS
Two and a half years after the merger, Chandra continues to detect X-ray emission at the location of GW170817. A comparably long-lived X-ray emission is rare in GRBs, and was reported only for long duration bursts, such as GRB 130427A (De Pasquale, et al. 2016) and GRB 980425 (Kouveliotou, et al. 2004). For a spectral index β=0.585, the extrapolation of the observed X-ray emission corresponds to F 606W ≈29.7±0.3 AB mag in the optical and ≈5±2 µJy at 3 GHz. For comparison, at the GW location HST /WFC3 can reach a 5 σ point-source sensitivity of F 606W ≈28 AB mag in four orbits (Lamb, et al. 2019), whereas a 6 hr long VLA observation can reach a 5 σ sensitivity of ≈10-15 µJy in S-band 3 . X-ray observations therefore remain the most powerful probe into the faintest stage of the GW counterpart.
In the latest epochs, the measured X-ray flux is higher than model predictions based on the earlier dataset , suggesting a shallower temporal decay. Contamination from an unrelated X-ray source seems unlikely. The probability of a background AGN of comparable flux is about 10 −4 arcsec −2 (Georgakakis, et al. 2008). The density of luminous X-ray sources within the galaxy is also relatively small, as can be directly seen from Fig. 1. The population of X-ray binaries in elliptical galaxies is in part associated to globular clusters, however deep HST observations find no globular cluster at the transient position (Troja, et al. 2017;Lamb, et al. 2019). The density of field X-ray binaries depends on the specific star formation rate (sSFR). Present systematic studies cover the range of log (sSFR)> −12.1 (Lehmer, et al. 2019), while NGC 4993 has a much lower value, log(sSFR< −13) (Im, et al. 2017). Assuming that the relationship established at higher vales of sSFR holds, 10 X-ray binaries with LX 3×10 38 erg s −1 are expected in NGC 4993. Taking into account the distribution of Xray sources as a function of their radial offset (Mineo, et al. 2014), we derive a chance alignment of ≈ 10 −3 arcsec −2 at the position of GW170817. Any significant departure from the jet model is likely inherent to the source, and could be caused by several factors, which we discuss below. Figure 3 compares the X-ray dataset to the range of jet model light curve predictions. The fit results are summarized in Table 2. Our previous best fit (Ryan, et al. 2019), based on the full first year data set, is shown by the dashed curve. The discrepancy between the new data and these earlier predictions is approximately 2σ, with the previous fit notably under-predicting the new observations. A refit of the full updated data set is also shown in Fig. 3, the solid bands denoting the distribution of X-ray flux estimated by the model. within the uncertainties with those from the first year data, although both the viewing angle θv and circumburst density n0 center on higher values than before. Both these increases can be understood on simple grounds. The early rise of the jet fixes the ratio θv/θc but leaves their absolute values relatively unconstrained (Ryan, et al. 2019;Nakar & Piran 2020). As the jet is slowly approaching the Sedov regime, the brighter than expected late X-ray emission requires a wide jet to contribute more flux. Indeed, Table 2 shows our fit value for the opening angle θc increased from 0.07 rad to 0.09 rad when the new observations were included. Since the early afterglow fixes θv/θc, the required viewing angle increases as well. The circumburst density is increased to keep the jet break at 160 d, compensating for the increased viewing angle which would otherwise push the jet break to a later time (Ryan, et al. 2019).

Jet
In Table 2, we compare the results of our modelling to additional observing constraints, which were not input into the fit. Ghirlanda, et al. (2019) constrained the size of the radio centroid at T0 + 207 d to δ < 2.5 mas at 90% confidence. All our models are safely within this limit. A more stringent constrains comes from the apparent velocity βapp of the center-of-brightness on the sky. A value of βapp = 4.1 ± 0.5 has been obtained from Very Large Baseline Interferometry (VLBI) by Mooley, et al. (2018c), measured between 75 d and 230 d after the burst. The model fit to the first year of data estimates βapp = 3.5 +1.2 −0.8 , consistent with the ob-served value. However, the updated fit significantly under predicts the observed centroid movement, estimating only βapp = 2.2 +0.5 −0.4 . This is largely due to the increased viewing angle, to which the superluminal apparent velocity is a sensitive function.
In our Gaussian jet model the observed motion of the radio afterglow centroid, which requires smaller viewing angles, appears therefore in slight tension with the late X-ray flux, which instead favors larger viewing angles. The tension could be alleviated if the afterglow light curve were able to flatten faster than our current modelling allows. Such an effect could originate from the dynamics of the GRB jet, changes to the emitted synchrotron spectrum, or possibly an additional emission component.
Because the spreading of GRB jets occurs during an intermediate dynamical regime between ultra-narrow highly relativistic flow and broad non-relativistic flow, the evolution of the jet during the spreading stage is more sensitive to the details of outflow geometry than either asymptotic limit of behaviour would suggest. This affects both multidimensional hydrodynamical simulations of jets and semianalytical models. Our model is based on a semi-analytical model for jet spreading (Ryan, et al. 2019;van Eerten, Zhang & MacFadyen 2010), and shares this sensitivity. For that reason, we also test the extreme assumption of no spreading at all. Such a jet is non-physical, but serves to bracket the range of jet model light curve predictions.
We ran a fit to the full dataset with a non-spreading Gaussian jet. The best fit (maximum posterior) light curve is shown in Figure 3 (solid line) and the summary of fit results are presented in Table 2. The non-spreading jet has a slower decay after the jet break and is more easily able to accommodate the late data points while requiring an earlier and broader peak. Changing the model assumption about jet spreading mostly affects our inferred values for the angles and circumburst density (see Table 2). These end up smaller, consistent with the previous estimates derived from the dataset at 360 d but outside the uncertainties from the fit to the full dataset. The apparent velocity increases to βapp = 4.3 +1.4 −0.9 due to the smaller viewing angle, and is consistent with the observed value. Although this model does not describe a realistic jet configuration, this fit serves to demonstrate that the interpretation of afterglow data at these late times is highly sensitive to the dynamics of jet spreading.
Both for jets with and without lateral spreading, the full transition to the non-relativistic regime takes tNR ≈ 10 4 days to complete and will not impact the light curve at the current time scale of observations for a reasonable range of model parameter values. The same holds for the appearance of the counter jet, which our models project to temporarily lead to a near-flat light curve between 3000-5000 days after the burst (at around 10 −16 erg cm −2 s −1 at X-ray frequencies and around 0.2 µJy at 3 GHz).
Rather than the divergence between model and data being due to limitations of the model, the jet dynamics might also genuinely change under changing external conditions, specifically a change in circumburst density. Analytical modeling for a homogeneous environment show that the flux below the cooling break scales proportional to circumburst number density n according to n 1/2 and n 0.4 (p = 2.2) in the relativistic and non-relativistic limit respectively (see Table 2. Fit results for the jet models. Col. 1 reports the parameters name and units. Col. 2: a Gaussian structured, spreading jet fit to the first 360 days of observations. Col. 3: identical jet model fit to all 940 days of data. Col. 4: a Gaussian jet with spreading artificially stopped. This model is not physical, but serves to bracket the diversity of possible behaviours of spreading jets. Col. 5: a spreading Gaussian jet with an additional constant X-ray flux. Notes -Marginalized posterior values for each fit parameter, the median and 68% confidence interval, from the MCMC runs are given in columns 2 -5, rows 1-8. Rows 9-11 give the marginalized posterior values for the total energy Etot, apparent velocity βapp measured between the VLBI observations (Mooley, et al. 2018c), and rms width of the centroid during the EVN observations Ghirlanda, et al. (2019) respectively, also with median and 68% confidence interval. The last three rows give the reducedχ 2 value of the maximum-posterior estimate (and degrees of freedom for each fit), the WAIC estimate of the expected log predictive density (elpd), and difference between the WAIC values and the spreading Gaussian Jet fit with standard error. A higher elpd indicates a model better able to predict the data.
e.g. Leventis, et al. 2012). In other words, it would merely take a factor four increase in density at distances beyond about a parsec from the merger site (the approximate distance traveled by the jet when observed at its light curve peak around 160 days) in order for the light curve baseline to drift towards a factor two increase, consistent with the latest observations. A change in the light curve slope can also occur if the synchrotron cooling break frequency enters or exits the Xray band. However, our structured jet modeling shows that both radio and X-ray light curves remain in the same spectral regime between injection break νm and cooling break νc throughout our observations, and that νc shifts upwards again after a closest approach to the X-ray band during the light curve peak around 160 days. This is exactly the same evolutionary pattern for νc as predicted across jet breaks from ultra-high resolution numerical hydrodynamics simulations (starting from top-hat initial conditions, see  (Figure 2), although, given the large error bars, it is not possible to draw too strong a conclusion about this similarity. We therefore find it unlikely that the cooling frequency νc affects the latest X-ray observations. Finally, it could be the case that the synchrotron parameters themselves evolve over time. For example, the value of p evolving closer to 2, as expected for non-relativistic shock speeds (Blandford & Ostriker 1978;Bell 1978), would indeed lead to a harder spectrum (from 0.585 for p = 2.17 to 0.5 for p = 2) and shallower temporal slope (since α in t −α equals p for a fast spreading jet, 3p/4 for a non-spreading jet and (15p − 21)/10 in the non-relativistic limit, see Zhang, & Meszaros 2004;Panaitescu & Kumar 2004;Frail, Waxman & Kulkarni 2000, respectively). However, this would also affect the overall flux normalization, which contains an εe(p − 2) term. Although the impacts of these effects on the light curve will be mitigated by the spread in emission arrival times from the blast wave, it would still require εe to co-evolve such that a substantial shift in baseline flux level is to be avoided. An updated broadband measurement of the slope of p could directly answer the question whether p is indeed evolving, but for now we conclude that the tentative flattening of the light curve has not been established as a generic prediction of such a scenario. . Posterior distribution on the ejecta velocity index k, assuming the kilonova afterglow contributes to the observed X-ray flux at 2.5 yrs (orange). Radio upper limits were also included in the fit. In purple, the posterior distribution on k if the X-ray flux remains above the current level up to 5 years after the merger. Table 2 presents the results of fitting an additive constant Xray flux to the spreading, Gaussian structured jet afterglow. In such a scenario the viewing angle θv and circumburst density n0 are somewhat reduced compared to their values from the jet alone, and consistent with the values derived from the 1 year dataset. The additional flux density at 5 keV is constrained to FX = (2.8 ± 1.2) × 10 −5 µJy, corresponding to (8 ± 3) × 10 −15 erg cm −2 s −1 , about half the observed flux at T0 + 939 d. The smaller viewing angle causes a larger apparent velocity, βapp = 2.7 +1.0 −0.6 , consistent with the observations of Mooley, et al. (2018c). The improvement in WAIC score between the jet plus constant and the standard jet is marginal (1.2 ± 1.4), and does not warrant the addition of another parameter in the model.

Kilonova Afterglow
We use the latest X-ray and radio observations to constrain the range of valid kilonova afterglow models 4 . In lieu of running a combined fit with both structured jet and kilonova afterglows, we use the structured jet plus constant fit (Sect. 4.1.1) as a measure of the possible contribution of the kilonova afterglow to the current epoch. We run a simple MCMC fit with the kilonova afterglow model to the X-ray flux FX ≈ 8 × 10 −15 erg cm −2 s −1 , as well as the latest radio upper limits. There are no other constraints apart from priors and the requirement the light curve be currently rising.
We focus our study on the emission arising from the fastest ejecta, often referred as the "blue" kilonova component, as it is expected to peak ealier and initially be brighter (Alexander, et al. 2017;. Our prior on Mej is a normal distribution with mean 2.25 × 10 −2 M and width 0.75 × 10 −2 M , as derived from the modeling of AT2017gfo (e.g. Arcavi, et al. 2017;Evans et al. 2017;Nicholl, et al. 2017;Kasen, et al. 2017;Pian et al. 2017;Tanvir et al. 2017;Troja, et al. 2017). Our prior on the minimum outflow velocity βmin is a normal distribution with mean 0.3 and width 0.05, as lower values would lead to delayed and dimmer peaks below our detection limits . The velocity distribution index k was given a uniform prior between 1 and 10. The circumburst density n0 was given a log-uniform prior between 10 −3 and 10 −1 cm −3 in agreement with the constraints from the jet model. The electron spectral index p was given a uniform prior between 2 and 3, while εe and εB were given log-uniform priors between 10 −5 and 1. We note these parameters are under no obligation to take identical values in both the structured jet and kilonova afterglow.
We find the current data set admits a broad range of kilonova models (Figure 4) and is insufficient to provide strong constraints on any of the parameters, including the velocity distribution index k. Figure 5 shows the posterior probability distribution on k. Essentially any value is consistent with current observations. Preliminary constraints (disfavoring k <6) were derived by Hajela, et al. (2019), our exploration of the parameters space finds instead a broader range of possible solutions. This result is consistent with the analysis presented in Hajela, et al. (2019), in particular their Fig. 5 showing a wide range of allowed values, but does not support the conclusion k ≥6. Higher values of k result in fainter initial emission and a steep rise to the peak flux. Lower values of k are instead brighter at earlier times with a slow rise to the final flux. These are easily brought in agreement with the current observations with a slight reduction of εe and εB (Figure 4). Continued monitoring of this target would therefore be critical to determine the rising slope of the kilonova afterglow component, and constrain the ejecta velocity profile.
Unfortunately, due to the large number of parameters and uncertainty in the physical properties of the kilonova blast wave, it is difficult to make robust conclusions about its afterglow emission at this time. Ultimately, the large uncertainty in the synchrotron parameters εe and εB dominate the analysis, and will only be overcome with successful observations. As shown in Figure 6, the same observing settings thus far adopted to monitor GW170817 probe the top 30% of the estimated flux distribution and could detect the kilonova afterglow under favorable conditions.

Energy injection from a pulsar
Another possibility of flattening the light curve is to invoke energy injection of a long-lived NS. This possibility was suggested by  to account for the X-ray variability around 160 days and to interpret some of the features in the kilonova AT2017gfo associated with GW170817 (Yu, Liu & Dai 2018;Wollaeger, et al. 2019). A long-lived NS central engine is allowed by the EM and GW observational data, as long as the surface dipole field strength is not very strong (Ai, et al. 2018) and the NS equation of state is stiff enough (Ai, ). For such a NS, the spindown time scale can be of the order of years, so that significant energy injection is still possible at the time of our observations. Indeed,  predicted the flattening of the lightcurve based on their model parameters to interpret the X-ray variability.
We consider a general energy injection law from the central engine, L(t) ∝ t −q , where q < 1 is needed to give a noticeable change of blastwave dynamics (Zhang & Mészáros 2001). We consider two possibilities. The first is that the spindown luminosity is injected into the blastwave as a Poynting flux. For GW170817/GRB 170817A, the current epoch is already in the post-jet-break phase since the light curve is already in the rapid decay regime. Let us assume that the blastwave is still in the relativistic regime and that sideways expansion is not important, one can derive an analytical model for decay slopes. For a constant density medium (which is relevant for NS-NS mergers), one has (Zhang 2018) 5 Γ ∝ t −(2+q)/8 , r ∝ t (2−q)/4 , νm ∝ t −(2+q)/2 , νc ∝ t (q−2)/2 . The peak flux density can be estimated as Fν,max ∝ r 3 B Γ[θ 2 j /(1/Γ) 2 ] ∝ r 3 Γ 4 ∝ t (2−5q)/4 . 5 The relevant parameters are for the jet core and an on-axis observer. For a structured jet with a large viewing angle like the case of GRB 170817A, these scalings are relevant after the jet core enters the line of sight, i.e. during the rapid decay phase. For the kilonova afterglow, we report the two upper bounds (95% and 68% confidence intervals) on the estimated flux distribution (cf. Figure 4). On the right we report the 5 σ sensitivity of typical observing settings for ATCA, VLA, and Chandra (CXO), as well as for the next generation ngVLA (Corsi, et al. 2019) and the Athena X-ray observatory.
For νm < ν < νc which is relevant for X-rays at such a late epoch, the flux density evolution should satisfy This expression is consistent with the pre-jet-break energy injection theory (Zhang, et al. 2006) if the edge effect correction factor [θ 2 j /(1/Γ) 2 ] is removed. For q = 0 relevant to pulsar injection in the pre-spindown phase, this gives Fν ∝ t (2−p)/2 , which is nearly flat (for our best fit p = 2.17, this gives t −0.085 ). This is consistent with the numerical result presented in Figure 6. For this first scenario, energy injection should be achromatic. The same flattening feature should appear in the radio band as well.
The second scenario invokes an internal dissipation of the pulsar wind, which has been manifested by the so-called "internal plateaus" as observed in both long (Troja, et al. 2007;Lü & Zhang 2014) and short (Rowlinson, et al. 2010;Lü, et al. 2015) GRBs. The temporal profile should directly follow ∝ t −q , which is also flat for q = 0. The light curve should be chromatic, as seen in GRB afterglows (Troja, et al. 2007;Rowlinson, et al. 2010;Lü & Zhang 2014;Lü, et al. 2015), and the radio band may not show a simultaneous flattening as the X-ray band. Since all the other flattening mechanisms (discussed earlier in Sect.3.1 and 3.2) also predict achromatic behaviors, a detection of chromatic behavior between X-ray and radio will provide a definite clue about a long-lived central engine.
If the flattening is indeed caused by energy injection of a long-lived pulsar, the spindown time scale should be at least this long, i.e. (Dai & Lu 1998;Zhang & Mészáros 2001) T sd ∼ (2 × 10 7 s) B −2 p,13 P 2 0,−3 > 1, 000 d, where Bp = 10 13 G Bp,13 is the surface polar magnetic field strength, and P0 = 1 ms P0,−3 is the initial spin period of the pulsar. This condition is readily satisfied if Bp is below a few times of 10 12 G, which is consistent with the constraints from other observations from this event (Ai, et al. 2018;Ai, Gao & Zhang 2020;. Within the energy injection model, lightcurve flattening appears when the injected energy exceeds the original energy in the blastwave, and the ceases when the total available spin energy is injected. According to our structured jet modeling, the total kinetic energy in the jet is ∼ 10 50 − 10 52 erg with medium value 5 × 10 50 erg. This is smaller than the typical available spin energy of a new-born millisecond pulsar from an NS-NS merger (typically a few 10 52 erg, but could be smaller due to possible a secular gravitational wave loss, Fan, Wu & Wei 2013;Gao, Zhang & Lü 2016). As a result, such an energy injection is expected if the merger product is indeed a long-lived neutron star. The injection energy may be up to a factor of a few to a few hundreds of the existing energy in the jet, so that the injection episode may last for years according to this model.

CONCLUSIONS
Whereas optical and radio emission from GW170817 have now faded below detection threshold, its X-ray counterpart continues to be visible at 2.5 years after the NS merger. Earlier predictions of the structured jet model systematically underestimate the latest Chandra detections. A Gaussian structured jet can still reproduce the afterglow temporal evolution by increasing the viewing angle to ≈30 • , although this updated model underpredicts the centroid motion, as constrained by high-resolution radio imaging. Alternatively, the slow X-ray decline could indicate a genuine new feature of the afterglow, originating from the dynamics of the GRB jet, changes to the emitted synchrotron spectrum, or possibly an additional emission component. The latter contribution is constrained by our modeling to FX ≈8×10 −15 erg cm −2 s −1 , corresponding to an X-ray luminosity LX ≈1.5×10 38 erg s −1 (0.3-10 keV). Continued energy injection by a long-lived central engine would cause a persistent flattening of the X-ray lightcurve. Depending on the origin of this emission (internal or external), the same flattening could be observed in the radio band. The observed behavior could mark instead the onset of a non-thermal "kilonova afterglow", produced by the interaction of the sub-relativistic merger ejecta with the surrounding medium. We find that the current dataset is not sufficient to meaningfully constrain any of the parameters, including the velocity distribution index k. Our results do not support earlier predictions of k ≥ 6 and find a wide range of allowed values. Future multi-band observations of this component would be essential to determine the velocity profile of the sub-relativistic ejecta, thus complementing earlier kilonova studies, based on the thermal optical/nIR emission.