Properties of jet and surrounding material of GW/GRB~170817A

TL;DR

This paper models the afterglows of GW/GRB 170817A using a multi-component jet structure, revealing insights into jet composition, beaming, and circum-burst environment based on broad band observations over a year.

Contribution

It introduces a multi-component jet model with a structured relativistic jet and cocoon, consistent with long-term afterglow observations and previous prompt emission analysis.

Findings

Jet has a core with Lorentz factor ~100 and a mildly relativistic cocoon.

Afterglow rise linked to low circum-burst density and jet expansion.

External shocks occurred >110 days post-merger, affecting jet evolution.

Abstract

We use published data in radio, optical and X-ray bands to analyze and model afterglows of GW/GRB 170817A. Our analysis is based on a phenomenological gamma-ray burst generator model which we previously used to study the prompt gamma-ray emission of this important transient. We find a multi-component model and a few of its variants that are consistent with broad band $\sim 1$ year observations of afterglows, once the contribution of kilonova in optical/IR band is taken into account. Considering beaming and off-axis view, we interpret the components of the model as approximately presenting the profile of a relativistic structured jet with a rapidly declining Lorentz factor from our line of sight, where it had a Lorentz factor of $\mathcal {O}(100)$, to outer boundaries, where it became a mildly relativistic cocoon with a relative velocity to light of $\sim 0.4-0.97$. Properties of the jet core obtained here are consistent with conclusions from analysis of the prompt gamma-ray emission. In particular, our results show that after prompt internal shocks the remnant of the jet retains in some extent its internal collimation and coherence. Slow rise of the afterglows can be associated to low density of circum-burst material and low column density of the jet. The long distance of external shocks from the merger, which could have been in part responsible for extensive thinning of the jet through expansion and energy dissipation before occurrence of external shocks is responsible for the peak of emission being at $\gtrsim 110$~days after the merger. We discuss implications of these observations for properties of circum-burst material around binary neutron stars. This analysis confirms our previous results showing that an outflow with a Lorentz factor of $\sim 2-5$ cannot explain observed afterglows without an additional X-ray source or significant absorption of optical/IR photons.

Tables

Table 1. Table 1: Parameter set of simulated models.

Comp.	mod.	${\gamma }_0^{\prime }$	$r_0$ (cm)	$\frac{\Delta r_0}{r_0}$	${\left(\frac{r}{r_0}\right)}_{max}$	$p$	${\gamma }_{cut}$	$\kappa$	$\delta$	${ϵ}_B$	${\alpha }_B$	${ϵ}_e{Y}_e$	${\alpha }_e$	$N^{\prime }$ (cm-3)	$n_c^{\prime }$ (cm-2)
Ultra. rel. (C1)	1	130	${10}^{16}$	${10}^{-7}$	1.5	1.8	100	-0.5	0.5	$0.08$	-1	0.1	-1	$0.04$	$5\times {10}^{22}$
	2	-	-	-	15	-	100	0.3	0.1	-	0	-	0	-	-
	2	-	-	-	20	-	100	0.4	0.05	-	1	-	1	-	-
Rel. (C2)	1	5	${10}^{16}$	${10}^{-6}$	2	2.1	100	-0.5	1	$0.08$	-1	0.1	-1	$0.04$	${10}^{23}$
	2	-	-	-	40	-	100	0.4	0.1	-	0	-	0	-	-
	2	-	-	-	100	-	100	0.5	1	-	1	-	1	-	-
Mildly rel. (C3)	1	1.06	$1.5\times {10}^{16}$	${10}^{-2}$	1.5	1.8	100	-0.5	1	$0.08$	-1	0.02	-1	$0.008$	${10}^{24}$
	2	-	-	-	10	-	100	0.	0.1	-	0	-	0	-	-
	2	-	-	-	10	-	100	1	1	-	1	-	1	-	-

Table 2. Table 2: Parameters of the phenomenological relativistic shock model

Model (mod.)	Model for evolution of active region with distance from central engine.
$r_0$ (cm)	Initial distance of shock front from central engine.
$\Delta r_0$	Initial (or final, depending on the model) thickness of active region.
$p$	Slope of power-law spectrum for accelerated electrons; See eq. (3.8) of (Ziaeepour & Gardner, 2011).
$p_1,p_2$	Slopes of double power-law spectrum for accelerated electrons; See eq. (3.14) of (Ziaeepour & Gardner, 2011).
${\gamma }_{cut}$	Cut-off Lorentz factor in power-law with exponential cutoff spectrum for accelerated electrons; See eq. (3.11) of (Ziaeepour & Gardner, 2011).
$\Gamma$	Lorentz factor of jet with respect to far observer.
$\delta$	Index in the model defined in eq. (3.29) of (Ziaeepour & Gardner, 2011).
$Y_e$	Electron yield defined as the ratio of electron (or proton) number density to baryon number density.
${ϵ}_e$	Fraction of the kinetic energy of falling baryons of fast shell transferred to leptons in the slow shell (defined in the slow shell frame).
${\alpha }_e$	Power index of ${ϵ}_e$ as a function of $r$.
${ϵ}_B$	Fraction of baryons kinetic energy transferred to induced magnetic field in the active region.
${\alpha }_B$	Power index of ${ϵ}_B$ as a function of $r$.
$N^{\prime }$	Baryon number density of slow shell.
$\kappa$	Power-law index for N’ dependence on $r^{\prime }$.
$n_c^{\prime }$	Column density of fast shell at $r_0^{\prime }$.

Equations

$$\frac{dP_{e}/d\Omega}{dP_{r}/d\Omega}=\frac{1}{\Gamma^{2}\biggl{(}1+\beta\cos(\theta)\biggr{)}}$$

$$dP/(\omega d\omega )\propto r^2\Delta r/\Gamma (r){\int }_{{\gamma }_m}^{\infty }d{\gamma }_e{n}_e^{\prime }({\gamma }_e)f({\gamma }_e,r)$$

$$\displaystyle\Delta r^{\prime}=\Delta r^{\prime}_{0}\biggl{(}\frac{\gamma^{\prime}_{0}\beta^{\prime}}{\beta^{\prime}_{0}\gamma^{\prime}}\biggr{)}^{\tau}\Theta(r^{\prime}-r^{\prime}_{0})\quad\text{dynamical model, Model = 0}$$

$$\displaystyle\Delta r^{\prime}=\Delta r^{\prime}_{\infty}\bigg{[}1-\biggl{(}\frac{r^{\prime}}{r^{\prime}_{0}}\biggr{)}^{-\delta}\biggr{]}\Theta(r^{\prime}-r^{\prime}_{0})\quad\text{Steady state model, Model = 1}$$

$$\displaystyle\Delta r^{\prime}=\Delta r^{\prime}_{0}\biggl{(}\frac{r^{\prime}}{r^{\prime}_{0}}\biggr{)}^{-\delta}\Theta(r^{\prime}-r^{\prime}_{0})\quad\text{Power-law model, Model = 2}$$

$$\displaystyle\Delta r^{\prime}=\Delta r_{\infty}\bigg{[}1-\exp(-\frac{\delta(r^{\prime}-r^{\prime}_{0})}{r^{\prime}_{0}})\biggr{]}\Theta(r-r^{\prime}_{0})\quad\text{Exponential model, Model = 3}$$

$$\displaystyle\Delta r^{\prime}=\Delta r^{\prime}_{0}\exp\biggl{(}-\delta\frac{r^{\prime}}{r^{\prime}_{0}}\biggr{)}\Theta(r^{\prime}-r^{\prime}_{0})\quad\text{Exponential decay model, Model = 4}$$

Full text

††pubyear: 2019††pagerange: Properties of jet and surrounding material of GW/GRB 170817A–LABEL:lastpage

1 Introduction

There is no general consensus about physics behind the unusual afterglows of the short GRB 170817A - the electromagnetic counterpart of the first detected Gravitational Wave (GW) from merger of a BNS (LIGO, et al., 2017a, b). This burst can be singled out by the faintness of its prompt gamma-ray and early X-ray afterglow, and its later brightening leading to the detection of a X-ray counterpart only after $\sim T+10$ days (Troja, et al., 2017; Margutti, et al., 2017; Evans, et al., 2017), where $T$ is the trigger time of the Fermi-GBM (Dingus, 1995) and Integral-SPI-ACS (Winkler, et al., 2003) about $2$ sec after the chirp of the GW from the merger (Goldstein, et al., 2017; Savchenko, et al., 2017).

The initial interpretation of these observations was off-axis view of an otherwise ordinary short GRB (Lazzati, et al., 2018; Troja, et al., 2017; Murguia-Berthier, et al., 2017; Alexander, et al., 2017). Another popular explanation was emission from break out of a cocoon - a mildly relativistic ejecta with a Lorentz factor $\Gamma\sim 2-3$ (Kasliwal, et al., 2017; Gottlieb, et al., 2017; Hallinan, et al., 2017; Nakar, et al., 2018; Piro & Kollmeier, 2017). However, the decline of flux in all 3 observed energy bands, i.e. radio (Mooley, et al., 2017; Dobie, et al., 2018; Alexander, et al., 2018; Mooley, et al., 2017, 2018b), optical (Lyman, et al., 2018; Rossi, et al., 2018; Lamb, et al., 2019), and X-ray (Margutti, et al., 2018; Hajela, et al., 2018; Troja, et al., 2018a; D’Avanzo, et al., 2018; Nynka, et al., 2018; Haggard, et al., 2018) after $\gtrsim T+200$ days is much earlier than the prediction of a side viewed jet (Lazzati, et al., 2016) (see also earlier versions of (Lazzati, et al., 2018)). Predictions of jet break out (Kasliwal, et al., 2017; Gottlieb, et al., 2017) and cocoon (Mooley, et al., 2017) models for 6 GHz band are compared with data taken up to $\sim T+110$ days in (Margutti, et al., 2018). We should emphasize that predictions of the break time are model dependent and some off-axis models such as those studied by (Rossi, et al., 2004; Lamb & Kobayashi, 2017) find a peak time $\sim\mathcal{O}(200)$ days, i.e consistent with peak time of GRB 170817A afterglows. But, these works present general models and are not adjusted to special properties of GW/GRB 170817A.

Gradually it became clear that the presence of a highly relativistic component in the outflow at late times is inevitable. For instance, using relativistic hydrodynamic simulations of (Duffell & MacFadyen, 2013; Duffell, et al., 2015) for determining characteristics of the BNS merger outflow and asymptotic formulation of external shocks and synchrotron emission by (Sari et al., 1998) for fitting the afterglows, (Margutti, et al., 2018; Xie, et al., 2018) found that an outflow with a narrow relativistic core having a Lorentz factor of $\sim 100$ and a sheath/side lob with a Lorentz factor of $\Gamma\sim 3-10$, where our line of sight passes through with an angle of $\sim 20^{\circ}$, can explain observations up to $\lesssim 140$ days. Following the detection of superluminal motion of the radio afterglow due to an oblique viewing angle, (Mooley, et al., 2018a; Ghirlanda, et al., 2019) concluded $\Gamma\sim 4$ for the outflow at $\sim 230$ days and estimated an initial Lorentz factor of $\sim 10$ for the jet at the time of prompt gamma-ray. They also estimated an off-axis angle of $\theta_{v}\sim 20^{\circ}$ for the line of sight. In a further work (Mooley, et al., 2018b) ruled out the cocoon/jet break out model suggested by some of these authors in their previous works and fit the spectrum with a phenomenological non-linear 2-component broken power-law expression (Beuermann, et al., 1999) in contrast to a simple power-law used e.g. in (Mooley, et al., 2017; Margutti, et al., 2018; Troja, et al., 2018b). Finally, using all the observations in radio, optical/IR and X-ray during the first year after the merger event, (Lamb, et al., 2019) found that the data can be fit by a 2-component jet model consisting of an ultra-relativistic component with $\Gamma\gtrsim 100$ and a relativistic component with $\Gamma\sim 5$.

In (Ziaeepour, 2018b) we used a phenomenological shock and synchrotron emission model to show that an outflow with a Lorentz factor of $\sim 2-3$ underestimates X-ray flux. The same phenomenological formulation was used in (Ziaeepour, 2018a) to model the prompt gamma-ray emission of GW/GRB 170817A. The range of Lorentz factor studied in (Ziaeepour, 2018b) is the same as those employed in the early literature on this burst - specially those associating the unusual afterglows to a mildly relativistic cocoon111In the literature “mildly relativistic” indicates a large range of Lorentz factors from $\Gamma-1\sim\mathcal{O}(0.1)$ to $\Gamma\sim\mathcal{O}(1)$.. The reason for this choice was the assumption that at $t>T+10$ days, i.e. well after internal shocks, the weak ultra-relativistic jet responsible for the faint prompt gamma-ray burst had dissipated its energy and its Lorentz factor had to be much smaller than $\Gamma\sim 10-100$ concluded in (Ziaeepour, 2018a). In line with the same argument, the observed deficiency of X-ray with respect to optical and radio concluded from simulations with aforementioned low Lorentz factor implied either an additional source of X-ray - for instance a contribution from the decay of radioactive elements produced by the kilonova - or a significant absorption of optical photons.

In the present work we use the same phenomenological model as the one employed in (Ziaeepour, 2018a, b), but we drop the assumption of a dissipated jet. We show that a multi-component model, including both ultra-relativistic and mildly relativistic components, and a kilonova, can explain all the data. The components of the model approximately present angular profiles of density and Lorentz factor of the polar ejecta from the BNS merger and its evolution. We use this model and properties of circum-burst material to investigate reasons behind the late brightening of the afterglows.

In Sec. 2 we describe the model and compare it with afterglow models of GW/GRB 170817A in the literature. Interpretation of the 3-component model is discussed in Sec. 3. Parameters of the phenomenological model are summarized in Appendix A. Degeneracies of these parameters are discussed in details in (Ziaeepour & Gardner, 2011; Ziaeepour, 2018a) and are not repeated here. Nonetheless, to investigate whether conclusions made in Sec. 3 can be significantly impacted by them, in Appendix B we present light curves of several variants of our best model and discuss their properties. Our results are summarized in Sec. 4.

2 Afterglow model

The phenomenological model of relativistic shocks and synchrotron-self-Compton emission of (Ziaeepour, 2009; Ziaeepour & Gardner, 2011) used in the present work is reviewed in (Ziaeepour, 2018a, b) and we do not repeat it here. Nonetheless, for the sake of self-sufficiency definition of parameters of the model are given in Table 2. Some details about how kinematics of the shock and dynamics of the emissions are modelled and related are given in Sec. 3.3.2

It is important to remind that there is a significant difference between our approach, in which a synthetic burst is generated for a set of input parameters characterizing the jet and its surrounding, and modelling of afterglows according to asymptotic power-law behaviour of light curves and spectra based on the original calculations of general aspects of synchrotron emission from external shocks by (Sari et al., 1998). Notably, in our approach simulated bursts explicitly depend on distance and column density of the jet and thereby give an assessment of these quantities. Although analysis of GW/GRB170817A afterglows e.g. by (Margutti, et al., 2017, 2018; Alexander, et al., 2018; Hotokezaka, et al., 2018) are based on the jet characteristics obtained from Magneto-Hydro-Dynamics (MHD) or relativistic hydrodynamics simulations, synchrotron emission is calculated according to the asymptotic formulation of (Sari et al., 1998), which estimates power-low behaviour of afterglows in a given energy band by comparing emission’s frequency with characteristic and cooling frequencies of accelerated electrons. However, backreaction of shocks and energy dissipation are not explicitly taken into account. The model used here is also phenomenological and an approximation. In particular, some of important processes and quantities, which cannot be easily formulated from first principles, are presented by parametrization and their initial values are chosen by hand. Nonetheless, the model takes into account in a systematic manner backreaction and evolution of physical properties important for the synchrotron/self-Compton emission. Moreover, the model is applicable to both internal (Ziaeepour & Gardner, 2011; Ziaeepour, 2018a) and external shocks, and thereby allows to compare and to verify consistency of parameters obtained from the two types of emission in the same framework.

Due to large number of parameters in the model and CPU time necessary for each simulation it is not possible to perform a systematic search for the best fit to data222To give an idea about calculation time, on a 3.06 GHz Intel Duo CPU T9900 processor a simulation with 3 time intervals

• regimes - takes about 20-30 minutes without calculating inverse Compton and about 2 times longer with inverse Compton. If only 5 parameters are changed on a lattice with 5 nodes for each, total calculation time without inverse Compton would be about $10^{5}$ minutes or about 130 days. In practice more than 5 parameters should be adjusted to find best models and the time necessary for a systematic search would be much longer than above estimation.. Nonetheless, despite their apparently arbitrariness, physically acceptable values of parameters are not completely random. The distance of external shocks from central source is determined by wind nebula surrounding progenitor neutron stars and its termination shock. In pulsars wind nebula extends to $\sim 10^{15}-10^{17}$ cm (Chatterjee & Cordes, 2004; Slane, 2017; Posselt, et al., 2018). But its dependence on the properties of neutron stars and their evolution is not well understood. Our simulations show that an initial distance of $\sim 10^{16}$cm, which is in the logarithmic middle of the range given here leads to acceptable fit to the GW/GRB 170817A data.

The density of circum-burst material on which the jet/outflow is shocked has a lower limit corresponding to the ISM density of $\lesssim 0.04$ cm*-3* in the host, concluded from the absence of significant neutral hydrogen in NGC 4993 (Hallinan, et al., 2017). The spectrum of accelerated electrons in the shock is also fairly constrained to $\sim 2$ by Particle In Cell (PIC) simulations (Spitkovsky, 2008). We do not consider any external magnetic field in the simulations presented here. The remaining parameters define the geometry of ejecta and surrounding material, and are adjusted by trial and error to fit the data as good as possible.

We find that a 3-component model consisting of a diluted ultra-relativistic jet with $\Gamma\sim\mathcal{O}(100)$, a relativistic outflow with $\Gamma\sim\mathcal{O}(10)$ and a mildly relativistic outflow/cocoon with $\Gamma-1\sim\mathcal{O}(0.1)$ can explain observations in radio and X-ray bands and satisfies upper limit constraint imposed on optical/IR data, see below for more details. From now on we call these components C1, C2, and C3, respectively. Parameters of the model for these components are listed in Table 1. We remind that as the parameter space was not systematically searched, values of parameters for the best model given here should be treated as order of magnitude estimations. For this reason we do not provide any uncertainty for them. Moreover, to see how variation of parameters can affect light curves and spectrum, and whether there is large degeneracies in the parameter space, which may invalidate interpretation of the model and its comparison with data, in Appendix B we present light curves of several variant models for each of the above components and compare them with the model presented in Table 1. Notably, the model presented in Fig. 7-a for component C3 has a Lorentz factor of $4$, that is similar to estimation from apparent superluminal motion of radio counterpart (Mooley, et al., 2018a), and fits radio data up to $\sim T+200$ days roughly as good as C3 in Table 1. However, at later times the latter provides a better fit. See Sec. 3.2 for further discussion of this subject.

2.1 Data used for comparison with models

For comparing models with observations we use published data from various sources:

Radio data is in $5-6$ GHz radio band taken from (D’Avanzo, et al., 2018) except for the last two points which are taken from (Dobie, et al., 2018) and (Mooley, et al., 2018b). The radio data in (D’Avanzo, et al., 2018) is reproduced from observations of (Hallinan, et al., 2017; Mooley, et al., 2017).

Optical/IR observations in $r$ and $i$ bands or energetically close bands i.e. $R$, HST $F606W$ and HST $F814W$. For $<T+10$ days the data is taken from (Soares-Santos, et al., 2017). The $R$ magnitude at $\sim T+10$ is taken from Pian, et al. (2017), $r$ band magnitude at $\sim T+110$ is from (Lyman, et al., 2018), HST $F606W$ at $\sim T+134$ from (D’Avanzo, et al., 2018) (originally from (Margutti, et al., 2018)), HST $F814W$ magnitude at $\sim T+160$ from (Rossi, et al., 2018) and HST $F606W$ magnitude at later times from (Lamb, et al., 2019). Magnitude $m$ is changed to flux density in ph/sec/cm2 using $F(\text{ph/sec/cm}^{2})=10^{-0.4(m+48.6)+23}(\Delta\lambda/\lambda)\times 1.509\times 10^{3}$.

X-ray data is in $0.3-8$ keV for Chandra observations and $0.3-10$ keV for data taken by the Swift-XRT and by the XMM-Newton. For the epoch before $\sim T+9$ days only 2 upper limit for X-ray flux is available: an upper limit for time interval $\lesssim T+10$ days from the Neil Gehrels Swift-XRT (Evans, et al., 2017) and an upper limit at $\sim T+2.2$ from Chandra observatory (Troja, et al., 2017). The X-ray counterpart was detected and followed up later and we use data from Chandra observations at $\sim T+9$ days (Margutti, et al., 2017), $\sim T+16$ days (Troja, et al., 2017), $\sim T+110$ (Margutti, et al., 2018), XMM-Newton observation at $\sim T+134$ days (D’Avanzo, et al., 2018), Chandra observations at $\sim T+260$ days (Nynka, et al., 2018; Troja, et al., 2018a) and at $\sim T+359$ days (Haggard, et al., 2018).

Despite heterogeneity of the data we did not attempt to homogenized it for not to add further uncertainties. This fact and mismatch between observations and simulated bands should be taken into account when models are compared with data.

After the detection of the prompt Gamma-ray spikes by the Fermi-GBM (Goldstein, et al., 2017) and the Integral-IBIS (Savchenko, et al., 2017) GW/GRB 170817 is not detected in high energy electromagnetic bands. Only an upper limit on any extended emission in $15-50$ keV band, based on averaging of background counts per 16 days during from $T\gtrsim 2$ sec to $\sim T+1$ year, is obtained from Swift Survey data (Krimm, 2018).

2.2 Comparison with observations

Fig. 1-a shows light curves of each component of the model and Fig. 1-b the total light curves in each band. We notice that X-ray and radio light curves have much better fit to data than simulated optical light curve. However, we know that optical data, specially at early times, is dominated by kilonova emission (Smartt, et al., 2017; Cowperthwaite, et al., 2017; Covino, et al., 2017; Nicholl, et al., 2017) and cannot be modeled with shock/synchrotron emission. Thus, presentation of early optical data in Fig. 1 and other light curve plots in this work is for the sake of completeness. Although a priori kilonova contribution can be modeled and removed, the residual would be model dependent. For this reason, we show optical data as it is observed, but use it as an upper limit for the GRB contribution which should be respected by any model of GRB 170817A afterglows. As the model presented in Fig. 1 (and some of its variants discussed in Appendix B) fit well radio and X-ray data, we presume that their prediction for optical emission should be reliable. Under this assumption, these models show that after $t_{kn}\sim T+200$ to $\sim T+300$ days - depending on the model - kilonova contribution in optical/IR emission was not anymore significant and the afterglow was dominated by synchrotron emission from external shocks of the relativistic outflow333Note that in addition to model dependence one has to take into account the heterogeneity of the optical/IR data and mismatch with the simulated energy band. For these reasons the time of dominance of synchrotron with respect to kilonova contribution given in the text should be considered as nominal rather than exact.. This estimation for the time of kilonova fainting $t_{kn}$ is consistent but somehow larger than that by (Waxman, et al., 2019), who find $t_{kn}\gtrsim 100$ days. However, due to degeneracies in the parameter space of their kilonova emission model (Waxman, et al., 2019) do not fit it to data. Moreover, synchrotron light curves used for estimating GRB contribution are one possibility between many, see the discussion of degeneracies in Appendix B. Other works on the evolution of the kilonova emission, e.g. (Kasen, et al., 2017; Smartt, et al., 2017; Metzger, et al., 2018; Kasen & Barnes, 2019) are limited to initial tens of days after merger and cannot be compared with late optical data. Thus, we conclude that at present a satisfactory model for late optical emission of kilonova is not available.

Another assessment of the performance of the model can be made by comparing simulated spectrum with the photometric spectrum reconstructed from observations shown in Fig. 2-a. It shows a good a consistency between simulated spectrum and the data. As expected, amplitude of optical emission at $\sim T+\mathcal{O}(10)$ days is higher than the model. Moreover, due to the dominant contribution of kilonova emission the pseudo-spectrum of energy flux shown in Fig. 2-b at this epoch is significantly different from those of later times. In this plot if we neglect optical data and use only radio and X-ray (dotted lines in Fig. 2-b), spectra of all observation epochs have similar behaviour. Of course, a 2-point pseudo-spectrum is a very crude presentation of the broad-band spectrum. Nonetheless, it is not affected by kilonova emission. Fig. 2-c) shows evolution of spectral slope, which is similar to afterglows of other GRBs, namely softer during earliest observations around $T+\mathcal{O}(10)$ days, gradually becomes harder until the peak of emissions around $T+110$ days, and finally softens at later times. The last data point in Fig. 2-c is obtained from radio and X-ray data with largest uncertainties (see Fig. 1-b) and apparent increase of slope and hardening of spectrum during last observation is very uncertain.

2.3 Comparison with other analyses

As we described in the Introduction, some other authors have analyzed late afterglows of GW/GRB 170817A by using multi-component models including an ultra-relativistic jet. For instance, (Margutti, et al., 2018; Alexander, et al., 2018) consider a structured jet with a top-hat ultra-relativistic component with $\Gamma\sim 100$ in the inner $\theta\lesssim 9^{\circ}$, where $\theta$ is angle with respect to symmetry axis of outflow, and a relativistic flow which its Lorentz factor decreases in the angular interval $10^{\circ}\lesssim\theta\lesssim 60^{\circ}$ and has a mean $\Gamma\sim 10$. The authors consider a line of sight angle $\theta_{v}\sim 20^{\circ}$ from the jet axis. This model is very similar to the model described in the latest version of (Lazzati, et al., 2018). But Lorentz factor and energy profile of the jet in the two works are different. (Lamb, et al., 2019) consider two profiles for the jet, one similar to 2-component model of (Margutti, et al., 2018; Alexander, et al., 2018). The other model assumes a jet with Gaussian energy and Lorentz factor profiles. Both models find a small central/core angle of $\lesssim 5^{\circ}$. In their 2-component model high and low Lorentz factors are $\sim 100$ and $\sim 5$, respectively. However, in the Gaussian model on-axis Lorentz factor can be as large as $\gtrsim 900$. (Mooley, et al., 2018b) model light curves of afterglows with a phenomenological function, which effectively has two coupled components. A priori they can be related to their Lorentz factor and density profiles. However, the complicated and non-linear form of this model make any physical interpretation difficult.

A question arises here: Is the presence of an ultra-relativistic component in our models and those cited in the previous paragraph inevitable ? The existence of such a component means that our line of sight passed through the ultra-relativistic core of the jet. However, several authors find that the data can be fit with synchrotron emission of a jet with a Lorentz factor $\lesssim 10$ in our direction. For instance, (Xie, et al., 2018) find several variants of such models and show that they fit the data up to $\lesssim 150$ days. However, none of their models fit earliest X-ray data and aside from radio bands their data set include only a couple of points in X-ray and optical/IR. The model of (Margutti, et al., 2018; Alexander, et al., 2018) reviewed in the previous paragraph has an ultra-relativistic core with $\Gamma\sim 100$ up to $\theta\sim 9^{\circ}$, but the viewing angle of $\lesssim 20^{\circ}$ found by these analyses means that the line of sight passes through a region with a Lorentz factor $\ll 100$ - they do not indicate the exact value of Lorentz factor on the line of sight. Similar to (Xie, et al., 2018) their model under estimate early X-ray data. To these examples one should add other models, namely cocoon (Mooley, et al., 2018a, b; Nakar, et al., 2018) and cocoon break out (Kasliwal, et al., 2017; Gottlieb, et al., 2017), which have much smaller Lorentz factors. None of these models is consistent with late X-ray data and over-estimate its flux. Their consistency with optical data is also uncertain. Interestingly, with publication of later data, a general trend in the literature toward models with larger Lorentz factors is discernible - specially in works using observations up to $\sim 200$ days and beyond. For example, in earlier versions of (Lazzati, et al., 2018) a Lorentz factor of $\lesssim 10$ was associated to the jet. But in the final version of this work, which includes more data, the value of Lorentz factor is increased to $\mathcal{O}(100)$.  (Lamb, et al., 2019) who use observations from earliest time to roughly one year after merger find a Lorentz factor $\gtrsim\mathcal{O}(100)$. In conclusion, it is unlikely that the necessity of an ultra-relativistic Lorentz factor along our line of sight found here be due to degeneracies or quirk of our models.

In addition to differences in Lorentz factor on the line of sight, predictions for time of kilonova decay $t_{kn}$ from modelling of GRB 170817A in the literature are not the same and are somehow shorter than $t_{kn}\sim 200-300$ days obtained from the model of Table 1 and its variants in Appendix B. Specifically, GRB afterglow models of (Lamb, et al., 2019; Wu & MacFadyen, 2019; Troja, et al., 2018c) seem to fit optical data at $\gtrsim T+100$ days. A priori this means that a $t_{kn}\lesssim 100$ days can be concluded from these works. However, they do not show early optical/IR data points and do not explicitly discuss the contribution of kilonova in the optical flux. In fact, it seems that they adjust the models to fit optical data after $\sim 100$ days. In this case, the time of kilonova decay in their models is a prior rather than a posterior conclusion derived from data. By contrast, as described earlier, in selection of plausible models here optical data was only considered as an upper limit and no attempt was made to fit simulated light curves to it.

3 Interpretation

3.1 Origin of ultra-relativistic component

In what concerns the Lorentz factor of the ultra-relativistic component of our model and similar models in the literature, they are consistent with analysis of (Ziaeepour, 2018a). Indeed, modelling of the prompt gamma-ray emission of GW/GRB 170817A shows that the best estimation for the Lorentz factor of the jet at the end of main gamma-ray spike is $\gtrsim 100$.

A priori it is not expected that a relativistic jet preserves its Lorentz factor during propagation from the site of internal shocks, which in the case of GW/GRB 170817A was at $r_{i}\equiv r_{0}\sim\mathcal{O}(1)\times 10^{10}$ cm from merger (Ziaeepour, 2018a), up to its collision with the ISM or circum-burst material at $r_{e}\equiv r_{0}\sim\mathcal{O}(1)\times 10^{16}$ cm, in which afterglows are generated. However, weaker late internal shocks, cooling of shocked material, and shocks on circum-merger material between $r_{i}$ and $r_{e}$ can dissipate kinetic energy of the jet and decelerate it. Indeed, weak gamma-ray and X-ray spikes and continuous emission are observed in both long (Willingale, et al., 2006) and short bursts (Montanari, et al., 2005; Norris, et al., 2010) and are interpreted as emission from internal shock of side lobes - high-latitude regions - of the jet (Takami, et al., 2007) or weaker internal shocks (Willingale, et al., 2006; Ziaeepour & Gardner, 2011).

Shocks on the material close to the BNS can be simulated by considering short distance external shocks similar to models presented in Fig. 5-d, -e, & -f, which have $r_{e}=r_{0}=10^{13}$ cm or $r_{e}=r_{0}=10^{15}$ cm. Because the value of parameter $\kappa$, which parametrizes variation of ISM/circum-burst material density is the same as C1 model (and its variants), it can be considered as extension of the same distribution. Thus these simulations presents interaction of the ultra-relativistic jet with material near the merger. An important properties of these simulation is that their emissions peak at early times. Thus, along with weak late internal shocks, they can a priori produce a plateau regime around $<T+10^{5}$ sec, which is observed in some short GRBs such as GRB 070724A, GRB170127B, and GRB 111020A. The latter GRB is a an interesting case because it occurred at redshift $0.02$ (Tunnicliffe, et al., 2014) and was only a few times brighter than GRB 170817A. Unfortunately, for the latter burst no data is available in this time interval and these emissions cannot be investigated444If such data existed, amplitude of models in Fig. 5-d & e could be adjusted to the data..

Simulations presented in this work and those reviewed in Sec. 2.3 are the evidence that despite energy loss, a fraction of ultra-relativistic jet survives without significant dissipation555One of shortcomings of the phenomenological model of (Ziaeepour, 2009; Ziaeepour & Gardner, 2011) is that it is deterministic and cannot give a probability or distribution for dissipated and non-dissipated particles in the jet. . Comparing the column density of ultra-relativistic component C1 in Table 1 (and its variant in Fig. 5-b) with the initial column density of the jet before internal shocks, which according to analysis of (Ziaeepour, 2018a) was $10^{25}$ $cm^{-2}$, shows that it is $\sim 200$ folds less than what it was pre-prompt emission. Assuming that dissipation of the second fainter and softer prompt spike corresponds to relativistic component C2 listed in Table 1 and its variants discussed in Sec. B.3, we observe that its column density had become $\sim 5-50$ folds lower, i.e. reduced from $5\times 10^{23}$ $cm^{-2}$ to $10^{23}$ $cm^{-2}$ in C2 in Table 1 or to $10^{22}$ $cm^{-2}$ for models in which Lorentz factor $\gamma^{\prime}_{0}=\Gamma=30$, i.e. the same as the second prompt gamma-ray peak. Evidently, there is no proof for any relation between C2 model and the second spike in the prompt gamma-ray because time resolution of late observations does not allow to discriminate between afterglows of different density shells. Nonetheless, it is reasonable to presume that late time slower component include contribution from slower and more easily dissipatable part of the prompt outflow. Alternatively, rather than considering these models as competitive, they may be considered as finer decomposition of jet at late times, where the model with $\Gamma=30$ presents a highly reduced remnant of second peak and C2 in Table 1 presents dissipated and laterally expanded remnant of both prompt shells. In this case, column densities and/or thickness of active regions of these components must be slightly smaller than what is given in Table 1 and Fig. 6 to make the total emission consistent with data, see Sec. B.3 for more details. In any case, these variant models and their parameter space are consistent with our warning in Sec. 2 that values of parameters must be considered as order of magnitude estimations. Another important conclusion from variant models presented in Appendix B is that it does not seem possible to explain observed data with a 1-component or even 2-component model666The necessity of a multi-component model is the clear evidence that the phenomenological formulation used here is too simplistic. Notably, it does not take into account anisotropic density and Lorentz factor of particles inside the jet. Asymptotic formulation of (Sari et al., 1998) has also this problem..

If the jet were an adiabatically expanding cone, its column density had to decline by a factor of $\sim(r_{i}/r_{e})^{-2}\sim 10^{12}$ where $r_{i}$ and $r_{e}$ are distance from center when internal and external shocks occurred, respectively. The much smaller dilation factors we find for the 3-component model means that the material inside the jet had preserved in some extent its internal coherence and collimation - most probably through imprinted electric and magnetic fields in the plasma - and its expansion was not completely adiabatic and free. Alternatively, accretion of material between $r_{i}$ and $r_{e}$ could a priori compensate reduction of column density due to an adiabatic expansion. Assuming a power-law variation of jet’s column density by accretion per unit fly length, i.e. $dn^{\prime}_{a}(r)/dr=(\alpha n^{\prime}_{a}(r_{0})/r_{0})(r/r_{0})^{\alpha-1}$ and neglecting loss of kinetic energy to radiation, from conservation of momentum we find: $n^{\prime}_{c}(r)=n^{\prime}_{0}(r/r_{0})^{-2}+n^{\prime}_{a}(r_{0})[(r/r_{0})^{\alpha}-1]$ and $\beta(r)=n^{\prime}_{0}\beta(r_{0})/[n^{\prime}_{0}+n^{\prime}_{a}(r_{0})(r/r0)^{2}((r/r_{0})^{\alpha}-1)]$. For $r/r_{0}=10^{6}$ and $n^{\prime}_{c}(r)/n^{\prime}_{0}\sim 200$, $\beta(r)\sim 2\times 10^{-10}\beta(r_{0})$. This means that the jet would be practically stopped. Therefore, small reduction of column density could not be due to accretion of circum-burst material. On the other hand, much smaller column density of C1 with respect to the prompt jet demonstrates that energy dissipation during prompt shocks due to the accretion of material close to the merger and/or internal interactions and loss of internal coherence had indeed occurred and led to lateral expansion, and thereby dilation of the jet’s core and reduction of its Lorentz factor, specially in its outer boundaries. Therefore, our assumption about dissipation of the ultra-relativistic jet in (Ziaeepour, 2018b) was justified. Furthermore, thanks to long and extended follow up of this transient, for the first time we are clearly detecting the tiny remnant of the ultra-relativistic core of the jet in a short GRB and can distinguish its signature from slower part of the polar outflow of the BNS merger.

3.2 Jet profile and viewing angle

Since the observation of unusually faint prompt gamma-ray of GRB 170817A, understanding the underlying physics has been the subject of significant debate. Specifically, the viewing angle of the ejecta had a central role in the proposed explanations.

Observation of gravitational waves led to an estimation for orbital inclination angle of $18^{\circ}\lesssim\theta_{in}\lesssim 27^{\circ}$ (Mandel, 2018). Moreover, superluminal motion of radio afterglow with an apparent speed of $\beta_{app}=4.1\pm 0.5$ is observed by (Mooley, et al., 2018a; Ghirlanda, et al., 2019) and in combination with information from light curves and standard shock and synchrotron model and jet simulation (Mooley, et al., 2018a) estimated an off-axis angle of $\sim 20^{\circ}\pm 5^{\circ}$ for the line of sight. The above value of $\beta_{app}$ constrains Lorentz factor of the source to $\Gamma\gtrsim 4$. In this case, C3 component in Table 1 is not consistent with the observed superluminal motion. Nonetheless, a variant of C3 with $\Gamma=4$ and smaller ISM/circum-burst density and active region thickness presented in Fig. 7-a satisfies superluminal constraint. It fits radio data as good as C3 in Table 1 up to $\sim T+230$ days, which coincides with the second observation epoch of (Mooley, et al., 2018a). However, at later times the model with smaller Lorentz factor is a better fit to the data.

On the other hand, some other issues must be considered. Formulation of the phenomenological model developed in (Ziaeepour, 2009; Ziaeepour & Gardner, 2011) assumes $\Gamma\gg 1$. Therefore, uncertainty of C3 simulation with $\Gamma$ close to 1 and its parameters should be larger than other components. Moreover, estimation of the error on $\beta_{app}$ measurement might have been too optimistic. Fig. 2 (and Fig. 1 of extended data) of (Mooley, et al., 2018a) shows both error ellipses and synthesized beam shape, and explicitly says that the source was not resolved. Giving the fact that the size of the image displacement is comparable to beam size, the size of 1-sigma uncertainty ellipse seems too small. In addition, lateral expansion of the cocoon between two observations used for determining $\beta_{app}$, which were separated by $\sim 150$ days, might have changed the position of image centroid (Lind & Blandford, 1985; Grano & Piran, 2012). It is expected that due to scattering and energy dissipation in the core of the jet, the cocoon becomes slower, denser and expands laterally (Nishikawa, et al., 2016). Therefore, its effective off-axis angle with respect to line of sight might have changed between two observations, leading to larger apparent displacements. Thus, although $\beta_{app}$ was certainly nonzero, it could be much smaller than published value, and thereby consistent with a smaller Lorentz factor. An alternative possibility is to consider both C3 in Fig. 7-a and C3 in Table 1. In this case the former presents the state of out flow at angles further than C2 and the latter its state at even larger angles. Considering the uncertainty of data and parameters, such a 4-component model is consistent with radio data. In any case, most of our conclusions remain valid for both $\Gamma\sim 1.1$ and $\Gamma\sim 4$ for the component C3, and both 3 and 4 component models. For this reason we only consider 3-component model and discuss differences in interpretations for low and high Lorentz factor C3 whenever necessary.

As a final evidence for plausibility of a mildly relativistic component with small Lorentz factor, we remind direct observation of a cocoon in the long GRB 171205A and its associated supernova SN2017uk (Izzo, et al., 2017), which measures $\beta\sim 0.3$ corresponding to $\Gamma\sim 1.05$, i.e. similar to C3 component in Table 1.

Giving the fact that synchrotron emission is highly directional and emission from a relativistic source is beamed, off-axis view of the jet has significant consequences for observations. For instance, assuming a uniform Lorentz factor across the emission surface, a far observer receives oblique radiation only from a cone with half angle $\theta\sim\arcsin(1/\Gamma)$. Therefore, if the jet of GW/GRB 170817 was structured, radiation received by an off-axis observer would be from jet’s slower wings rather than faster core, and thereby dominated by photons with lower energies because they had been produced by weaker and less boosted shocks in the jet’s high latitude (boundary) region. Defining maximum visibility angle $\theta_{max}\equiv\arcsin(1/\Gamma)$, its value for components C1, C2, and C3 of the model listed in Table 1 are $\theta_{max_{1}}\sim 0.5^{\circ}$, $\theta_{max_{2}}\sim 11.5^{\circ}$ and $\theta_{max_{3}}\sim 65^{\circ}$ or $14.5^{\circ}$ for C3 in 1 and its alternative in Fig. 7-a), respectively.

Fig. 1 shows that each of the components of the model dominates the total emission in one of the three energy bands with observational data, and higher the Lorentz factor of the component, higher is the energy of its dominant emission among the three components 777For the time being, the simulation code used in this work uses an analytical expression for determining synchrotron/self-Compton flux. In the simulation code terms depending on higher order of $\theta$, angle between emitting element and line of sight, are neglected. This is a good approximation when $\Gamma\gg 1$. Under this approximation $\theta$ dependence is only through a $(cos(\theta)+\beta)\leq(cos(\theta)+1)$ factor, which must be integrated between $\theta_{0}\geq-\theta_{max}$ and $\theta_{1}\leq\theta_{max}$. However, angular size of emitting surface may be smaller than $2\theta_{max}$. In this case integration over maximum visible angle over-estimate the flux. But the difference would be at most a factor of few and comparable to other uncertainties of the model. Indeed, for this and other simplifications and approximations applied to the model that we should consider parameters as order of magnitude estimations. Taking into account the observed inclination of the BNS orbit and superluminal motion of radio counterpart, we interpret the 3-component model as a structured jet, in which each component approximately presents characteristics of the jet and its shocks on the surrounding material from our line of sight up to outer boundary of the outflow. The value of $\theta_{max}$ for the components are consistent with this interpretation888Although formulation of relativistic shocks and synchrotron emission in (Ziaeepour, 2009; Ziaeepour & Gardner, 2011) is more systematic than standard approach, it remains very much simplified and phenomenological. For instance, it considers a uniform column density and Lorentz factor for the jet. However, GRMHD simulations show that these assumptions are not true and column density is a function of both Lorentz factor and angle. It is why three components are necessary to explain the data. Therefore, this division and characteristics of components must be considered as an effective and simplified description.. Accordingly, their Lorentz factor presents azimuthal variation of velocity inside the polar outflow of the merger up to a $\cos\theta$ factor, that is $\Gamma_{simul,i}=\Gamma_{i}(\theta_{i})\cos\theta_{i}$, where $\theta_{i}$ is angle between centroid of component $i$ of the model and the line of sight. The value of $\theta_{max}$ for the components constrains $\theta_{i}\lesssim\theta_{max}$ and thereby the effect of projection on the estimation of parameters. Specifically, for C1 it is negligible and for C2 is about $10\%$. If we consider the model with $\Gamma=4$ for C3, the effect of off-axis is again about $10\%$. But for C3 model in Table 1 it can be as large as $42\%$.

Fig. 3 shows a schematic presentation of the structured jet, components of the model, and their positions with respect to our line of sight and symmetry axis of the outflow, assumed to have the highest Lorentz factor. Simulations of particle acceleration by transfer of Poynting to kinetic energy in the polar outflow of mergers show that the direction of maximum acceleration somehow deviates from the rotation and magnetic field axes (Komissarov, et al., 2009). Thus, using the estimation of orbit inclination, the angle between our line of sight and direction of maximum Lorentz factor $\theta_{v}$ can be constrained as $\mathcal{O}(1)^{\circ}\lesssim\theta_{v}\lesssim 27^{\circ}$, independent of jet model. This interval is also consistent with estimations of viewing angle by (Margutti, et al., 2018; Alexander, et al., 2018; Mooley, et al., 2018a)999In Fig. 3 jet axis and the line of sight are assumed to be in the same projected side with respect to the rotation axis. Thus $\theta_{v}<\theta_{in}$. More generally, considering a cone around the line of sight with rotation axis on its surface, any line passing inside of its cross-section will have a smaller angle with line of sight than the rotation axis. Therefore, for $\theta_{v}>\theta_{in}$ the jet axis must be outside the cone. In this case the line of sight would be even further from jet axis than what is depicted in Fig. 3, and closer to non-relativistic or mildly-relativistic cocoon. But, as discussed in the Introduction, such setup is not able to explain observations.

The knowledge of orbital inclination is not enough to constrain Lorentz factor profile up to its invisible core of the jet. Nonetheless, the sharp reduction of Lorentz factor from C1 which must be approximately along the line sight to C3 at $\lesssim 65^{\circ}$ or $\lesssim 14^{\circ}$ - depending on which C3 model is used - means that the visible part of the jet has an exponential $\Gamma\sim\Gamma_{max}\exp(-(\theta+\theta_{v})/\theta_{\tau})$ or Gaussian $\Gamma\sim\Gamma_{max}\exp(-(\theta+\theta_{v})^{2}/2\sigma^{2}_{\theta})$ profile, where $\theta$ is off-axis with respect to the line of sight, and $\theta_{v}$ is the angle between jet’s symmetry axis and line of sight. In absence of any knowledge about centroid of components, we must find the range of parameters which satisfy constraints on $\theta_{i}$, and on parameters of the profile, namely $\theta_{v}>0$, $\sigma^{2}_{\theta}>0$, and $\theta_{\tau}>0$, respectively for Gaussian and exponential profiles.

Fig. 4 shows the consistency region for a Gaussian profile for both the model of Table 1 and its C3 variant shown in Fig. 7-a. Assuming $\theta_{1}\sim 0.25^{\circ}$, for the model in Table 1 we find: $\theta_{2}\sim 9^{\circ}-11^{\circ}$, $\theta_{3}\sim 12^{\circ}-15^{\circ}$ for $\Gamma_{max}<1000$. Interestingly, this figure shows that the allowed values for viewing angle $\theta_{v}$ are strongly correlated with $\Gamma_{max}$ and are restricted to $5^{\circ}\lesssim\theta_{v}\lesssim 7^{\circ}$ for $\Gamma_{max}\sim 250$ and $14^{\circ}\lesssim\theta_{v}\lesssim 18^{\circ}$ for $\Gamma_{max}\sim 1000$. For second model the range of allowed angles are even more restricted: $\theta_{v}\sim 7.5^{\circ}-8^{\circ}$, for $\Gamma_{max}\sim 250$ and $\theta_{v}\sim 8^{\circ}-15^{\circ}$, for $\Gamma_{max}\sim 1000$, $\theta_{2}\sim 10.5^{\circ}$, $\theta_{3}\sim 11.5^{\circ}$, for $\Gamma_{max}<1000$. Considering the estimation of orbit inclination, the range of values for viewing angle obtained here confirms simulations of particle accelerations (Komissarov, et al., 2009) that show a deviation between magnetic field - assumed to be the same as rotation axis - and maximum of acceleration. In the case of GW/GRB 170817A the deviation was $\sim[7^{\circ}-15^{\circ}]$.

It is straightforward to see that in the case of exponential distribution $\theta_{v}$ and $\Gamma_{max}$ are fully degenerate. This degeneracy induces a relation between Lorentz factor and centroid of three components of the model, i.e. $(\theta_{2}-\theta_{1})/(\theta_{3}-\theta_{2}))=\ln(\Gamma_{1}/\Gamma_{2})/\ln(\Gamma_{2}/\Gamma_{3})$. Using this relation and imposing $\Gamma_{max}\sim 500$ and $\theta_{1}=[0.25^{\circ}-0.5^{\circ}]$, we find $\theta_{v}=[0.06^{\circ}-0.3^{\circ}]$, $\theta_{2}=[0.75^{\circ}-1.5^{\circ}]$, $\theta_{3}=[20.7^{\circ}-41.4^{\circ}]$, $\theta_{\tau}=[0.06^{\circ}-0.3^{\circ}]$ for the model in Table 1. For the alternative model $\theta_{2}$ is the same, $\theta_{v}\sim 0.06^{\circ}$, $\theta_{3}=[0.87^{\circ}-1.74^{\circ}]$, $\theta_{\tau}=[0.08^{\circ}-0.4^{\circ}]$. Small off-axis angles found for this profile are inconsistent with observed superluminal motion. Even for $\Gamma_{max}\sim 1000$ off-axis angles remain too small to be consistent with a superluminal motion. Therefore, we rule out an exponential profile.

The 3-component model provides also information about density profile of the jet. Indeed in agreement with simulations of jet acceleration to ultra-relativistic velocities (Komissarov, et al., 2009), the density in the core of the accelerated jet, where Lorentz factor is maximum, is much lower than its surrounding. Moreover, as described earlier, at the time of external shocks the jet had been already subjected to partial dissipation. Thus, a large fraction of its content was scattered, slowed down, and acquired significant transverse momentum. These processes should make side lobes - the cocoon - slower and denser, and increase gradient of variation in the outflow. These expectations are consistent with the model presented here. In addition, fraction of kinetic energy of the jet transferred to electric and magnetic fields during shocks in outer boundary of the outflow are expected to be smaller. C3 component which is only mildly relativistic satisfies these expectations. However, Table 1 shows that the distance of this component from the center at the onset of external shocks is $\sim 50\%$ larger than the other components. This may have multiple origins: uncertainties of simulations when $\Gamma\gtrsim 1$ and/or delay in the formation of shocks due to small velocity of C3 and low density of ISM/circum-burst material.

3.3 Delayed brightening

A significant difference between the afterglow of GW/GRB 170817A and other short bursts is its seemingly intrinsic faintness at early times, i.e. at $\gtrsim T+10^{5}$ sec, when X-ray and optical follow up began, and its brightening at late times. However, none of these conclusions are certain. Of course, as we showed in (Ziaeepour, 2018a), prompt emission of this burst was intrinsically the faintest among all GRBs with known redshift. Therefore, its early afterglow had to be equally faint. Indeed, X-ray flux upper limits of $\sim 10^{-13}$ erg sec*-1* cm*-2* from observations of the Swift-XRT at $\gtrsim T+1.6$ days in $0.3-10$ keV (Evans, et al., 2017) and $\sim 10^{-14}$ erg sec*-1* cm*-2* in $0.3-8$ keV from Chandra observations at $\sim T+2.2$ days (Troja, et al., 2017) are intrinsically the most stringent among short bursts with redshift (Fong, et al., 2018). However, early X-ray afterglow of both short and long GRBs have variety of slopes (Kann, 2012; Fong, et al., 2015), see also simulations in Appendix B. Therefore, in absence of early observations, GW/GRB 170817A cannot be classified as an early dark burst.

Late brightening of the afterglow of this burst is not unique either. For instance, a slight brightening may have been observed at $\sim T+10$ days for GRB 130603B and its accompanying kilonova (Fong, et al., 2013). However, the excess found by (Fong, et al., 2013) was with respect to the model fitted to the data rather than in the observed flux. In any case, no other GRB has ever followed up for as long as GW/GRB 170817A. Therefore, we do not know how uncommon is the brightening of afterglows up to $\gtrsim T+200$ days. Afterglows of GRBs up to few thousands of seconds are most probably a superposition of weak internal shocks in what remains from the relativistic jet after the main prompt shock, and emission from external shocks (Willingale, et al., 2006). Because a highly relativistic jet follows closely the prompt photons, the delay between their arrival to external material in general would be short - of order of few tens of seconds or less for ultra-relativistic jets and typical distances to circum-burst material. Assuming a narrow prompt spike, the delay between the prompt and onset of afterglow emission is $\Delta t\sim r_{e}/2c\Gamma^{2}$, where $c$ is speed of light. For the model of Table 1 the delay is $\sim 4$ sec. For alternative models discussed in Appendix B with shorter distances to central object this delay is even shorter. Therefore, the early afterglows, if they were observed, were superimposed on a decaying tail emission from internal shocks and time evolution of total emission depended on relative strength of decaying prompt tail and ascending afterglow. The absence of early brightening of afterglows in most GRBs means that they are dominated by tail emission of internal shocks. There are however exceptions. For instance, in long GRB 110213A (D’Elia, et al., 2011) rapid decay of the prompt X-ray made ascending and peak of X-ray afterglow visible. But such cases are rare. More frequent case, specially in short GRBs, is the presence of a plateau. Best examples are GRB 070724A (Ziaeepour, et al., 2007) and GRB 070809 (Marshall, et al., 2007), which have an early X-ray light curve with plateau similar to long GRBs.

Our simulations show that in the case of GW/GRB170817A the slow rise of the afterglow is due to low density of ISM/circum-burst material and low column density of the jet. The weakness of the jet was in part due to intrinsically low density and low Lorentz factor of polar ejecta - at least in our direction - and in part the result of large distance of surrounding material from center, that is $\sim 1000$ AU (see Table 1) rather than e.g. $\sim 100$ AU for the termination shock in solar system or $\sim 200$ AU for recently detected NIR emitting material around the isolated neutron star RXJ0806.4-4123 (Posselt, et al., 2018). The initial suggestions about off-axis view and brightening of afterglows when the jet become dissipated and its content scattered to our line of sight, is not consistent with relatively early break of light curves in all three energy bands after $\sim T+110$ days. In fact, high Lorentz factor of X-ray emitting component found here and in the literature means that it is emitted very close to our line of sight. Therefore, the effect of distance and densities may be stronger than off-axis. It is however cautious to consider that this conclusion may be somehow biased and correlated to spatial resolution of our simulations, because as Table 2 shows, the initial/final width of synchrotron emitting region is defined as a fraction of initial distance of the shock from central source. Nonetheless, estimation of shock distance for this GRB in the literature and for other GRBs with various methods, such as cooling of thermal emission (Olivares, et al., 2012) show that the range of distances obtained from simulation of GRBs (Ziaeepour & Gardner, 2011) with the code used here are realistic.

3.3.1 Comparison with standard afterglow model

There is no essential difference between model of (Ziaeepour, 2009; Ziaeepour & Gardner, 2011) and Blandford-McKee approach used by (Sari et al., 1998). In both cases the shock is assumed to be radiative, i.e. falling material to the shock front (in its rest-frame) loses its kinetic energy as synchrotron radiation and makes the jet heavier. The only difference is that (Ziaeepour, 2009) formulation takes into account the back-reaction of this process in a systematic way and also allows evolution of induced fields, which cannot be formulated from first principles, in a parametric manner.

We remind that asymptotic behaviour of external shock emission formulated in (Sari et al., 1998), which is used in works reviews in Sec. 2.3, considers a uniform spherical ejecta. Therefore, conclusions about viewing angle of the jet is based on the value of Lorentz factor and beaming of emissions from a relativistic source (Rybicki & Lightman, 2004):

[TABLE]

where $P_{e}$ and $P_{r}$ are emitted and received power in the observer frame, and $\theta$ is the angle between the line of light and the boost direction.

Comparing simulated spectrum shown Fig. 2-a, which is up to uncertainties fully consistent with observations shown in the same plots, with Fig. 1-a & b of (Sari et al., 1998)101010We remind that spectra shown in Fig. 2 includes an additional energy factor with respect to those in Fig. 1-a & b of (Sari et al., 1998)., we find that the total spectrum is consistent with a fast cooling model with $\nu_{c}<\nu_{m}\sim\nu_{x}$. However, as power of electrons Lorentz factor distribution $p\sim 2$ in all the components models of the model, the spectrum can be also classified as slow cooling with $\nu_{m}<\nu_{c}\sim\nu_{x}$. Spectra in Fig. 2-a show that for each component the dominant emission band is very close to the peak of its spectrum - thus consistent with being dominant - and satisfies relations similar to those of the total spectrum given above. In particular, similarity of components spectra and their sum demonstrates why many of analysis in the literature reviewed in Sec. 2.3, specially those modelling GRB 170817A with a single relativistic component, find that mission in radio, optical/IR, and X-ray bands belong to a same power-law spectrum with a negative slope.

Synchrotron cooling frequency $\nu_{c}$ depends on the induced magnetic field, which in turn depends on the density of ISM; larger the density of accumulated (fallen) material, larger magnetic field and higher synchrotron characteristic frequency. It can be shown, using equations (4) and (6) of (Sari et al., 1998), that $\nu_{c}\propto n^{-1}t^{-3/2}$, where $n$ is density of accelerated charges. Therefore, for a given $\nu_{c}$ lower density means slower emission variation and longer duration of emission because of smaller induced magnetic field. However, in (Sari et al., 1998) formulation the effect of jet column density is not explicit. We remind that the formation of a shock needs a significant density difference in the opposite sides of a boundary. Therefore, column density and density of ISM/circum-burst material must be considered together. Slow evolution of emissions needs relatively low ISM/circumburst material and/or low density jet. The length of the jet/ejecta is also important. For the same column density, a longer jet is more diluted and its collision with circum-burst material generates less turbulence and less energy is transferred to induced fields. Moreover, if densities of colliding shells are low, the faster shell must sweep a longer distance inside the slow shell to accumulate material and form a shock front. In model C1 in Table 1 ISM/circum-burst material density is in the middle of the range estimated for short bursts (Fong, et al., 2015). However, column density of jet is low. In the variant model in Fig. 5-b, which is also a good fit to data, ISM/circum-burst material density is lower, but a longer active region is assumed. Thus, for the ultra-relativistic component C1, the low density of colliding material may be the reason for long distance of external shocks from central source. By contrast, other components, which had larger column densities, have been radially more extended.

The total spectrum of the model in Fig. 2-a does not show how $\sim\nu_{c}$ and $\nu_{m}$ change with time. Nonetheless, evolution of spectrum discussed in Sec. 2 show that, as expected, after the peak of emission $\sim\nu_{c}$ and $\nu_{m}$ of each component become smaller than what was their values before the peak.

In summary, we conclude that the early X-ray afterglow of GW/GRB 170817A at $\lesssim T+10^{5}$ sec had to be faint and dominated by weak internal shocks. It declined quickly. A consequence of this conclusion, which unfortunately due to the lack of data cannot be verified, is that the claimed excess of UV emission at $\sim T+1.6$ days (Smartt, et al., 2017; Covino, et al., 2017; Cowperthwaite, et al., 2017) was indeed from kilonova rather than the afterglow of GRB 170817A. The origin of excess may be heating of a strong wind by neutrinos and strong magnetic field of the short lived HMNS (Metzger, et al., 2018)

3.3.2 Slope of afterglow rise

We notice that after initial fast rise the slope of light curves shown in Figs. 1, 5, 6, 7 in the region far from the peak is $\sim 1\pm 0.3$ and does not vary strongly111111We thank the anonymous referee for bringing this point to our attention.. To understand this behaviour we first remind that the asymptotic formulation of standard afterglow model by (Sari et al., 1998) predicts a constant rising slope of $\sim 0.5$ for energy bands close but on the low energy wing with respect to the spectral peak. As we discussed in Sec. 3.3.1, in the multi-component models of Table 1 and Appendix B the dominant emission band is close to the spectral peak. Therefore, having a roughly similar rising slope in these models is consistent with the asymptotic afterglow model. However, the latter underestimate the slope.

One reason for difference between predictions of the two models can be the fact that (Sari et al., 1998) formulation depends only on the density of accelerated particles in the jet, which generate the shocked induced electromagnetic fields. The density of circum-burst material $N^{\prime}$ on which the jet is shocked does not explicitly appear in the formulation. To see the importance of $N^{\prime}$ we notice that in the C1 alternative model of Fig. 5-a, in which only $N^{\prime}$ is less than that of C1 in Table 1, the rise slope of the X-ray light curve is $\sim 0.7$ and the peak flux is smaller. As explained in Sec. 3.3.1, time/radius dependence of quantities which determine strength of the shock, such as ISM/circum-burst density and fraction of kinetic energy transferred to electrons and to induced magnetic field - parametrized by $\kappa$, $\alpha_{e}$ and $\alpha_{B}$, respectively - influence the slope of emission rise and the observed differences with the standard model. However, they are not the only parameters which influence the rise of emission. For instance, in the main C1 model the slope after initial fast rise is $\sim 1$. But it increases gradually up to the peak of the X-ray light curve.

A quantitative explanation of the rise slope in the framework of the model of (Ziaeepour, 2009; Ziaeepour & Gardner, 2011) needs an analytical solution of the model, which due to its complexity is not available. Nonetheless, we can use some of analytical results obtained in (Ziaeepour, 2009; Ziaeepour & Gardner, 2011) to understand which physical properties are involved in the rise of light curves and its dependence on the parameters of the model.

The phenomenological model of (Ziaeepour, 2009; Ziaeepour & Gardner, 2011) divides the process of relativistic shocks and synchrotron/self-Compton to two parts: 1) kinematic of the ejecta; 2) synchrotron/self-Compton emission. It calculates the evolution of kinematic by taking into account the total energy dissipated as synchrotron/self-Compton and obtain perturbative solutions for Lorentz factor $\Gamma(r)$ and column density $n_{c}(r)=n(r)\Delta r$, where $r$ is an average distance of a narrow relativistic ejecta from center and $\Delta r$ is the width of shocked region - the active region. Under this approximation time $t$ and $r$ are not independent and $r(t)-r(t_{0})=c\int_{t_{0}}^{t}\beta(t^{\prime})dt^{\prime}$. The results of this step is used to calculate synchrotron/self-Compton emission, from well known textbook formulation of these processes, with the difference that it takes into account variation of electron and baryon densities in the jet/fast shell, their varying Lorentz factor, and possible radial density inhomogeneities of the slow shell/ISM. Time/radius dependence of these quantities connect kinematical and dynamical parts of the model. Additionally, in the case of internal shocks a total boost of colliding shells must be take into account. However, this phenomenological formulation is unable to determine how the width of active region $\Delta r$ changes with time/radius. Nonetheless considering formation and decay of a shock front, its initial width of active region should be zero. It should grow gradually to a maximum and decline to zero at the end of shocks. A series of phenomenological models presenting such evolution is defined in Appendix A.

Under these approximations the synchrotron energy flux can be written as:

[TABLE]

where $n^{\prime}_{e}$ is the density of electrons and $f$ is a function of $r$ and electrons Lorentz factor $\gamma_{e}$, see e.g. equation (62) in (Ziaeepour, 2009) for full expression. The $r^{2}$ factor changes densities to total emission from volume of the active region. It reflects the fact that for the same density, column density, and jet opening angle, the total amount of emitting material is larger if the average radius of the shell is larger. Therefore, in addition to quantities related to micro-physics of the ejecta and environment, geometrical setup, namely the initial distance of external shocks $r_{0}$, extension of the ejecta and evolution of the thickness of active region $\Delta r$ are crucial for the rise and peak amplitude of emissions, see equation 30 and Appendix A in (Ziaeepour, 2009) for technical details. In particular, for the same densities and other parameters, longer distance to center leads to longer rise time and higher and later peak emission.

The term $n^{\prime}_{e}\Delta r$ in (2) can be considered as column density of electrons with Lorentz factor $\gamma_{e}$. Normalization factor of the distribution of electrons Lorentz factor is proportional to density of active regions $n^{\prime}=n^{\prime}_{c}/\Delta r$ and inversely proportional to minimum Lorentz factor of electrons $\gamma_{m}\propto\epsilon_{e}\propto(r/r_{0})^{\alpha_{e}}$. However, due to non-linear term $f(\gamma_{e},r)$ factors of $\Delta r$ in (2) do not cancel each others. Evolution of $\Gamma(r)$ is also dependent on the phenomenological exponent of $\delta$, $\alpha_{e}$, $\alpha_{B}$ and $\kappa$, see equation (19) in (Ziaeepour & Gardner, 2011). As the values of these indices are the same in all the models discussed here - except the model shown in Fig. 5-e in which $\kappa=0$ - it is not a surprise that light curves have roughly the same rise slope. Nonetheless, due to the nonlinearity of evolution equations the slope in not exactly the same in all the models. In addition, as mentioned earlier, it changes with time. This fact reflects the influence of complicated dynamics, which cannot be estimated by a simple asymptotic power-law. For instance, the model in Fig. 5-e (dotted line) has a faster rise than C1 in Table 1 because $\kappa=0$, meaning that the density of circum-burst material in this model does not decline at large distances from the merger. By contrast, in the model of Fig. 5-a, which has the same $\kappa$ as C1, the rise of flux is slower because the density of circum-burst material is much smaller than in C1. This shows how initial conditions and constant quantities in the model influence the dynamics of the emission.

3.4 Material surrounding BNS

In the phenomenological model of (Ziaeepour, 2009; Ziaeepour & Gardner, 2011) the density of circum-burst material and its variation with distance are defined by parameters $n^{\prime}$ and $\kappa$, respectively. Table 1 shows that $\kappa\neq 0$ and $n^{\prime}$ is not the same for all components of the model. This is probably an evidence that circum-burst material was not only the ISM, which a priori should be independent of the merger and approximately uniform. Thus, additional material should have been present.

The origin of circum-burst material and its properties can be traced back to the evolution of progenitor neutron stars. In young neutron stars and pulsars the distance to wind Termination Shock (TS) is $R_{TS}=\sqrt{\dot{E}/(4\pi\eta P)}$ where $\dot{E}$ is the change in the rotational kinetic energy and $P$ is average pressure in the wind nebula surrounding the neutron star (Slane, 2017). After initial expansion of nebula and establishment of an approximately steady state condition pressure inside the nebula is balanced with the ISM pressure and extension of pulsar nebula is stabilized to $R_{TS}\sim\mathcal{O}(0.1)$ pc (Slane, 2017). However, in old neutron stars reduction of glitching activities and dissipation of magnetic field gradually decreases $\dot{E}$ and may reduce $R_{TS}$. Recent observation of thermal material at a relatively short distance of $\sim 200$ AU$~{}\sim 3\times 10^{15}$ cm around the isolated neutron star RXJ0806.4-4123 with an age of $\sim 10$ Myr (Posselt, et al., 2018) is an example of such cases. On the other hand, during BNS formation and merger if neutron stars were initially at a distance $\gtrsim R_{TS}$ and if there were still a remnant of their wind nebula around them, its content as well as any other material would be disrupted. Moreover, during early stages of inspiral glitching activities and mass ejection might have been resumed due to perturbation of neutron stars crust and might have replenished the nebula. These processes can explain the putative additional surrounding material and its anisotropic distribution according to the model presented here.

To estimate column density of this material we can use distribution used in the model i.e. $N^{\prime}(r)=N^{\prime}(r_{0})(r/r_{0})^{-}\kappa$. For model C1 this estimation gives a column density of $\sim 4\times 10^{14}$ cm*-2*, i.e. much smaller than column density of C1 which has the smallest column density among components of the model. This is smaller than swept material in the first $\sim 3\times 10^{5}$ sec after the onset of the external shocks and is completely negligible. See also Fig. 5-d & e for simulation of shocks on the ISM/circum-burst material with the same density as C1 but at low lower distances from center . On the other hand, if we consider much denser circum-burst material at shorter distances, much higher X-ray flux generated by the shocks violates upper limits at $\sim T+2$ days, see Fig. 5-f for an example with $N^{\prime}=4$ cm*-3*, i.e. 2 orders of magnitude larger than C1 in Table 1, in which the early X-ray flux is by about 2 orders of magnitude larger than upper limits. In such a model the column density of material between sites of internal and external shocks would be $\sim 4\times 10^{16}$ cm*-2*. Using these simulations and estimations, we conclude that the column density of material inside $r_{e}\sim 10^{16}$ was $<\mathcal{O}(1)\times 10^{15}$ cm*-2* or equivalently its average density was $<0.4$ cm*-3*.

In conclusion, diversity of short GRB afterglows reflect age, history, and environmental differences of progenitor neutron stars. In the case of GW/GRB 170817A we notice that $N^{\prime}$ is much smaller than typical values predicted for young neutron stars/pulsars (Slane, 2017). This means that the progenitor stars were most probably old and had lost most of material they had ejected during their youth. This conclusion is consistent with analysis of the prompt gamma-ray of this transient (Ziaeepour, 2018a). The relatively large distance of circum-burst material may have several reasons. But with available information about progenitors it is not possible to pin down the dominant cause.

4 Outline

In this work we used the same phenomenological formulation which had been used to model the prompt gamma-ray emission of GW/GRB 170817A to analyze its afterglows. We found a 3-component model presenting a structured relativistic jet which its collision with circum-burst material at a distance of $\sim 1000$ AU generated observed afterglows. It reproduces radio and X-ray light curves and photometric spectrum in these bands. Its optical emission is consistent with the dominance of the kilonova emission in this band up to $\gtrsim 200-300$ days after the merger event, where this range covers degenerate models studied here.

Additionally, this analysis helps understand physical conditions around the progenitor BNS before and after their merger. In particular, it shows that a small fraction of the prompt ultra-relativistic jet had survived internal shocks, and despite dissipation and lateral expansion it had preserved in some extent its internal coherence up to a distance of $\sim 10^{16}$ cm from central source, where it collided with surrounding material. Another result of the model studied here is that despite oblique view of the jet, our line of sight passed through an ultra-relativistic region. This conclusion is consistent with works in the literature which fit observations up to $200$ days post merger and beyond. It seems that an ultra-relativistic component is indispensable for a proper fit of the whole data. Other components of the model, which their emissions are dominantly in low energies, are interpreted as approximately presenting side lobes/cocoon of a structured jet at the time of external shocks.

The line of sight angle is model dependent and correlated to the jet profile and its maximum - on-axis - Lorentz factor. In the models studied here it can be as small as $\sim 5^{\circ}$ for a maximum Lorentz factor of $\sim 250$ or as large as $\sim 18^{\circ}$ if on-axis Lorentz factor approached $1000$. Considering observed orbital inclination of the BNS and superluminal motion of the radio counterpart, and marginalizing on the models and their parameters, the line of sight angle can be estimated as $5^{\circ}-18^{\circ}$.

We identified low density of circum-burst material from merger, its low density and low density of the jet as factors leading to slow rise of afterglows. Relatively long distance of the circum-burst materials from center might have been in part the reason for extended expansion of the jet up to longer distances than in typical short bursts and its dilation, which in addition to slow rise of afterglows, pushed the peak of emission to $\gtrsim 110$ days, rather than a couple of days seen in short GRBs with a X-ray plateau. Unfortunately, in absence of long duration follow up of short GRBs we are not able to access whether late shocks and brightening is an exception in this burst or a common behaviour of many short GRBs, i.e. those for which an early X-ray plateau is not observed.

Acknowledgment:

The author thanks Hans Krimm for providing the unpublished upper limit of the Swift-BAT on the late gamma-ray emission from GW/GRB 170817, and Amy Lien and Scott Barthelmy from Swift Science Team for their help to acquire this data.

Appendix A Definition of parameters and evolution of active region

Note that in the formulation model in (Ziaeepour, 2009; Ziaeepour & Gardner, 2011) $\Gamma$ is the Lorentz factor slow shell with respect to a far observer and $\gamma^{\prime}_{0}$ indicates the relative Lorentz factor fast and slow shells. In external shocks on the ISM or circumburst material, their velocity is negligible and $\Gamma=1$. For this reason, here we have replaced $\gamma^{\prime}_{0}$ with $Gamma$ which is usually used in the literature and indicates the Lorentz factor of outflow at the beginning of external shocks.

As explained in Sec. 3.3.2 the evolution of $\Delta r^{\prime}(r^{\prime})$ cannot be determined from first principles. For this reason we consider the following phenomenological models:

[TABLE]

The column called Model in Table 2 refers to these evolution models for $\Delta r^{\prime}(r^{\prime})$. The initial width $\Delta r^{\prime}(r^{\prime}_{0})$ in Model = 1 & 3 is zero. Therefore, they are suitable for description of initial formation of an active region in internal or external shocks. Other models are suitable for describing more moderate growth or decline of the active region. In Table 1 the column $mod.$ indicates which evolution rule is used in a simulation regime - as defined in the foot notes of this table - using model number given in 3-7.

Appendix B Alternative models for the structured jet of GW/GRB 170817A

As mentioned in Sec. 2 the large number of parameters in the phenomenological formulation used here to model afterglows of GW/GRB 170817A does not allow to perform a systematic search in the parameters space to optimize selected models. Moreover, based on physical arguments there must be correlation between parameters. Unfortunately in absence of a first-principle formalism they cannot be easily removed and induce degeneracies in the model.

To investigate how degeneracies affect the model of afterglows presented here and whether they can alter our conclusions, Figs. 5, 6, 7 show some variants of the model presented in Table 1. In most of these simulations only distance of shock from merger, Lorentz factor, column density of the jet, thickness of shocked (active) region, and density of ISM/circumburst material, are changed and other parameters are kept the same as ones shown in Table 1. In some models index of energy distribution of electrons $p$ and fraction of energy transferred to electric field $\epsilon_{e}$ are also changed.

Although some of these models fit the data as good as the model of Table 1, they have some issues which we will discuss case by case. For instance, a variant of C1 with larger $\Delta r_{0}/r_{0}$ and smaller ISM/circum-burst density presented in Fig. 5-b is degenerate with C1 in Table 1, but future observations of X-ray afterglow can distinguish between them. A variant of C2 with larger $\Delta r_{0}/r_{0}$ and smaller $\epsilon_{b}$ shown in Fig. 6-d has a slightly better fit to the data. But it is not exactly a replacement of C2 and can be considered as a finer division of the structured jet to components.

In the following subsections we summarize properties of variant models separately for each component.

B.1 Variants of C1

Fig. 5 shows variants of component C1. In model a) density of ISM/circum-burst material is smaller and consequently its X-ray emission is not sufficient. Model b) has larger $\Delta r_{0}/r_{0}$ but smaller ISM/circum-burst density, which falls in the lower half of the range of ISM/circum-burst material density estimated for other short GRBs (Fong, et al., 2015). This model is an example of a model degenerate with C1 in Table 1. However, it has flatter X-ray light curve at later times, i.e. after $T+4\times 10^{7}$ sec $\sim T+470$ days. Thus, future observations may distinguish between this model and C1 in Table 1. None of models in Fig. 5-a or -b can explain the radio afterglow. Fig. 5-c includes two models with smaller Lorentz factors, larger $\Delta r_{0}/r_{0}$, and larger ISM/circumburst and/or jet column densities. Both models have inconsistent X-ray and optical light curves and insufficient radio emission. We also notice that although $r_{0}$ in these models is the same as C1 in Table 1, the time of peaks in the same bands are not the same in the 3 models.

Fig. 5-d presents two models, one with a shorter initial radius $r_{0}$ for external shocks and the other with a longer distance. The most important difference of these models with C1 in Table 1 is the time of peak emission. Although the amplitude of X-ray light curves of these models are low, it can be adjusted by increasing $\Delta r_{0}/r_{0}$ and/or column density. Therefore, the essential problem of these models is the position of light curves peaks. Model e) has $r_{0}$ less than C1 but larger than the first model in Fig. 5-d and a jet column density one order of magnitude larger than C1. Similar to models in Fig. 5-d the time of peak emission and its amplitude are inconsistent with data. Finally, model f) has similar distance as the model in Fig. 5-e but higher jet column density and ISM/circum-burst material density. It is too bright in X-ray and optical and peaks are too early.

B.2 Variants of C2

As mentioned in Sec. 2, optical and IR emissions are dominantly from kilonova and only during the latest observations at $>T+200$ days contribution of the GRB might have become significant. For this reason optical data should be used as an upper limit for GRB contribution. It is also more difficult to select the best model for C2 and quantify characteristics of the component/section of the jet which emitted dominantly in optical bands.

Fig. 6 shows light curves of components and their sum for 4 variant models to C2 in Table 1. Model a) has a smaller $\Delta r_{0}/r_{0}$ and predicts a lower contribution from GRB in the optical/IR emission than C2. Model b) has a larger $\Delta r_{0}/r_{0}$ than C2 but smaller ISM/circum-burst density and fraction of kinetic energy transferred to accelerated electrons. This model fits the last 3 optical/IR observations and has a slightly better fit - $\chi^{2}\approx 1.57$ for C2 against $\chi^{2}\approx 0.05$ for this model (only the last 3 data points are considered). However, the model in Fig. 6-b may be also interpreted as slightly over-estimating optical emission at the epoch of latest 3 observations.

Models c) and d) are simulated with $\Gamma=30$, i.e. similar to the second prompt gamma-ray peak at the end of internal shocks (Ziaeepour, 2018a). The purpose for this choice is to see whether C2 can be a remnant of shells or sector of jet creating the second prompt peak. Interestingly, they fit X-ray light curve sightly better than model C1 in Table 1 because they have a significant X-ray mission at late times. But c) slightly over-estimates optical emission. This inconsistency can be resolved by e.g. slightly reducing the width of active region or column density. However, such adjustments decrease its X-ray contribution too and neutralizes advantage of this model with respect to the model of Table 1. Model d) has both X-ray and optical well consistent with the data. However, the density of ISM/circum-burst material in this model is 5 folds less than in C1 component. Giving the fact that radiation from models c) and d) should come from an azimuthal angle $\lesssim 2^{\circ}$ further from C1, such a large variation of ISM/circum-burst density seems unrealistic. By contrast, if in place of C1 we use its variant in Fig. 5-b, which has much smaller ISM/circum-burst density closer to that of model d) here, then together with d) for C2 they make an overall consistent model. In any case, column density of c) and d) models is much smaller than other components. In a picture in which the 3-component model crudely presents profile of the jet, a component with a tiny column density can be included to more prominent ones and neglected. As mentioned in the footnote 5, a better model of the jet should include continuous variation of column density and Lorentz factor with azimuthal angle.

B.3 Variants of C3

Fig. 7 shows variant models for component C3. Models 7-a, -b, -c have a Lorentz factor of 4, i.e. similar to estimation of (Mooley, et al., 2018a) from observation of superluminal movement of radio counterpart. Model a) has a shorter initial distance $r_{0}$ from center, thinner active region, and lower ISM/circum-burst material density, and denser jet. Radio light curve of this model is consistent with data up to $\sim T+2\times 10^{7}$ sec $\sim 230$  days. This date corresponds to the second epoch of observations reported in (Mooley, et al., 2018a), which was used to measure apparent superluminal displacement of the source. However, C3 model in Table 1 fits later observations better. Nonetheless, considering uncertainties of the data, this model a priori remains an acceptable alternative. On the other hand, a column density of $\sim 10^{25}$ cm*-2* in a) seems too large to be realistic. Indeed this is equal to the column density of the ultra-relativistic component at much closer distances to the merger and before prompt internal shocks. Energy dissipation in the outflow during its propagation from site of internal shocks up to location of external shocks should have reduced column density unless it had kept its coherence and collimation, and thereby its kinetic energy and Lorentz factor. But reduction of $\Gamma$ from 130 to just 4 is inconsistent with these assumptions. Another way to reconcile model a) with the overall picture is to assume that it is dominated by less accelerated, initially unshocked material ejected by the merger at azimuthal angles between ultra-relativistic component and slow kilonova disk/torus. A component with these properties may be present in some GRMHD simulations, for instance those by (Dionysopoulou, et al., 2015; Kiuchi, et al., 2015). In any case, replacing C3 in Table 1 with a) does not significantly modify conclusions of Sec. 3.

Models b) and c) have larger ISM/circum-burst density and over-produce both optical and radio emissions. Models d), e) and f) have a Lorentz factor slightly larger than C3 in Table 1 and $r_{0}$ shorter by $50\%$. None of them produces enough radio emission even with a column density $n^{\prime}_{c}$ as large as $10^{25}$ cm*-2*. These models show the interplay between Lorentz factor, densities, and thickness of active synchrotron emitting region.