Numerical Simulations of the Jet Dynamics and Synchrotron Radiation of Binary Neutron Star Merger Event GW170817/GRB 170817A

Our numerical setup captures the features of BNS mergers that shape the dynamics of the relativistic outflow and its radiative signature. General relativistic magnetohydrodynamics simulations of BNS mergers indicate that between 10⁻⁴ and 10⁻² solar mass (M_⊙) of neutron star materials are ejected during the coalescence, forming a quasi-spherical debris cloud. The cloud expands mildly relativistically, with typical radial velocity $\sim 0.15\mbox{--}0.25{\rm{c}}$ (e.g., Hotokezaka et al. 2013; Shibata et al. 2017). The modeling of the "kilonova" emission associated with GW170817 reveals that $\sim {10}^{-2}\,{M}_{\odot }$ of neutron-rich materials were ejected during the coalescence. We use this cloud mass in our simulations. The ejecta cloud has a slightly oblate geometry and radial stratification; most of its mass is confined in a slow-moving core, while a small amount of mass $[{10}^{-4}]{M}_{\odot }$ lies in an extended fast-moving tail,

$$\begin{eqnarray}&&{\rho }_{c}(r,\theta )={\rho }_{c}{(r/{r}_{c})}^{-2}(1/4+{\sin }^{3}\theta )\quad r\lt {r}_{c},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rho }_{t}(r)={\rho }_{t}{(r/{r}_{c})}^{-n}\quad {r}_{c}\lt r\lt 4{r}_{c},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&v(r)={v}_{c}r/{r}_{c}\quad r\lt 4{r}_{c}.\end{eqnarray} \tag{ 3 }$$

The density values ${\rho }_{c}$ and ${\rho }_{t}$ are calculated based on the total mass of the slow-moving core and the fast-moving tail. The density power-law index n is set to 8, the same value adopted in Kasliwal et al. (2017). The fast-moving tail is predicted to be the outcome of the shock that forms during the first collision between the merging neutron stars (Kyutoku et al. 2014; Nakar et al. 2018). ${v}_{c}=0.2\,{\rm{c}}$ sets the maximal velocity of the core, and ${r}_{c}=1.3\times {10}^{9}\,\mathrm{cm}$ is the core radius. This initial condition is similar to those in Kasliwal et al. (2017) and Gottlieb et al. (2018a). Both the narrow- and wide-engine models are evolved in the same merger cloud.

The central engine is initiated when the cloud has evolved for 1 s after the BNS coalescence. We make this choice because it yields GRB prompt emission compatible with the $1.7\,{\rm{s}}$ time delay between the observation of the GW chirp and the sGRB signal (Gottlieb et al. 2018a; Kasliwal et al. 2017; Nakar et al. 2018). The total engine energy is fixed at $5\times {10}^{50}\,\mathrm{erg}$ (per hemisphere), corresponding to roughly 6% of the rest mass energy of the merger cloud. Our simulation parameters are summarized in Table 1.

Table 1.  Hydroparameters for the Narrow-Engine and Wide-Engine Models

Variable	Narrow Engine	Wide Engine
${M}_{\mathrm{cloud}}(\mathrm{DS})$	$[0.01]{M}_{\odot }$	$[0.01]{M}_{\odot }$
${L}_{\mathrm{jet}}(\mathrm{SS})$	$2.6\times {10}^{50}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$	$2.6\times {10}^{50}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$
${t}_{\mathrm{jet}}$	$2\,{\rm{s}}$	$2\,{\rm{s}}$
${{\rm{\Gamma }}}_{0}$	10	1.02
η	100	20
${\theta }_{\mathrm{jet}}$	0.1	0.35
${n}_{\mathrm{ism}}$	${10}^{-4},{10}^{-5}\,{\mathrm{cm}}^{-3}$	${10}^{-5}\,{\mathrm{cm}}^{-3}$

Note. DS(SS) represents double-sided (single-sided). We use the same merger cloud, jet engine luminosity ${L}_{\mathrm{jet}}$, and engine duration ${t}_{\mathrm{jet}}$ for both the narrow-engine and the wide-engine models. The cloud mass ${M}_{\mathrm{cloud}}$ is $\sim [0.01]{M}_{\odot }$. The initial Lorentz factor ${{\rm{\Gamma }}}_{0}$, the specific enthalpy η, and the half-opening angle ${\theta }_{\mathrm{jet}}$ of the jet engine are set to different values between these two models.

Download table as:  ASCIITypeset image

We conduct 2D axially symmetric relativistic hydrodynamic simulations with the moving-mesh code JET (Duffell & MacFadyen 2013). The jet engine is injected as a source term for both models. For the narrow-engine model, we choose the jet engine as a nozzle with a circular profile (Duffell et al. 2015). For the wide-engine model, we adopt the same injection method of Kasliwal et al. (2017) and Gottlieb et al. (2018a), where cylindrical nozzles are used. Throughout, we use an ideal gas equation of state with adiabatic index 4/3 for simplicity and to allow for comparison with other simulations of this event (e.g., Nakar et al. 2018). We have performed simulations using the RC equation of state (Ryu et al. 2006) for comparison and find that our results do not qualitatively change.

Our simulations take place on a spherical grid with ${N}_{\theta }=160$ zones, evenly distributed in polar angle over the half-sphere. The central engine is modeled by injecting relativistic flow near the cloud center. The inner boundary is located at $4\times {10}^{7}\,\mathrm{cm}$. The radial grid is logarithmically spaced so that the cell aspect ratio is close to 1. Each simulation cell face moves radially with the flow. We adopt an adaptive mesh refinement (AMR) scheme that dynamically refines regions with high Lorentz factor.

In modeling a successful jet, we inject hot relativistic material within a narrow opening angle (see Table 1). The jet drills through the dense core of the merger cloud and breaks out highly overpressurized. It drives sideways expansion in the fast-moving lower-density tail of the merger cloud. Eventually, the outflow escapes the cloud altogether, at a radius $\sim 2.4\times {10}^{11}\,\mathrm{cm}$. GRB prompt emission photons are released from the vicinity of this breakout radius (Gottlieb et al. 2018a; Kasliwal et al. 2017; Nakar et al. 2018). Along the propagation direction, the relativistic GRB ejecta shocks the slower-moving merger debris ahead of it. The internal collision compresses the outflow into a very thin ultrarelativistic core. Meanwhile, the rapid lateral expansion of the sideways shock accelerates a mildly relativistic cocoon of neutron star materials, extending to a large lateral angle, as shown in Figure 1.

Zoom In Zoom Out Reset image size

After having emerged from the cloud, the jet has developed an angular dependent structure. The angular distribution of the total energy (kinetic plus thermal) is shown qualitatively in Figure 1 and quantitatively in Figure 2. The jet contains an ultrarelativistic core (${\rm{\Gamma }}\sim 100$) and a mildly relativistic sheath (${\rm{\Gamma }}\sim 10$). In Figure 2, we differentiate between the relativistic shell (${\rm{\Gamma }}\gt 1.1$) and the entire domain, labeled as shell and domain, respectively.

Zoom In Zoom Out Reset image size

We find that the angular energy distribution (${dE}/d{\rm{\Omega }}$) for both components is well described by the quasi-Gaussian profile,

$$\begin{eqnarray}&&{dE}/d{\rm{\Omega }}={\epsilon }_{0}{e}^{-{(\theta /{\theta }_{c})}^{\alpha }}.\end{eqnarray} \tag{ 4 }$$

This angular energy distribution is different from the top-hat model typically used in the fitting of GRB afterglow light curves (e.g., van Eerten et al. 2012); it better resembles the model described in Zhang et al. (2004). The total energy and the energy-averaged Lorentz factor,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{avg}}\equiv \displaystyle \frac{\int {\rm{\Gamma }}{EdV}}{\int {EdV}},\end{eqnarray} \tag{ 5 }$$

of the relativistic shell maintain their initial angular structure for a long period of time, $\sim {10}^{8}\,{\rm{s}}$. In Equation (5), E is the local energy density (measured in the lab frame) and Γ is the Lorentz factor of the fluid element. The maximum isotropic equivalent energy of the structured jet is ${E}_{\mathrm{iso},\mathrm{peak}}=4\pi {\epsilon }_{0}\approx {10}^{53}\,\mathrm{erg}$. Within an opening angle ${\theta }_{c}=0.15$, the average isotropic equivalent energy of the relativistic core is ${E}_{\mathrm{iso},\mathrm{avg}}=6\times {10}^{52}\,\mathrm{erg}$, larger than the average isotropic equivalent energy ${E}_{{\rm{k}},\mathrm{iso}}\approx (1\mbox{--}3)\times {10}^{51}\,\mathrm{erg}$ inferred for typical short GRBs, but still within the observed range (Fong et al. 2015).

The angular structure develops as a result of overpressurized relativistic ejecta escaping the merger cloud into the relatively dilute ambient medium. This results in significant lateral expansion (as depicted in Figure 1), in addition to radial acceleration. The jet propagating into the ambient medium consists of a shock-heated, baryon-clean core, surrounded by a shock-heated sheath of NS merger ejecta materials.

After the jet is launched by the central engine, it accelerates by converting its internal energy into kinetic energy. Over the course of $10\,{\rm{s}}$ (as measured in the lab frame, see the bottom panel of Figure 3), the jet attains its terminal Lorentz factor $\eta =100$ (which is also the specific internal enthalpy of the engine material at the jet base). During its propagation through the ejecta cloud, the jet performs work on it. In order to determine how the energy is partitioned during this phase of the evolution, we have computed the thermal energy E_t and the kinetic energy E_k for each of three components: the jet material, the shocked merger cloud (sometimes referred to as "cocoon" material), and the ISM. E_t and E_k are given by

$$\begin{eqnarray}{E}_{t,i} & = & \displaystyle \int [(e+p){{\rm{\Gamma }}}^{2}-p]{s}_{i}\,{\rm{d}}{V},\\ {E}_{k,i} & = & \displaystyle \int [{\rm{\Gamma }}({\rm{\Gamma }}-1)\rho ]{s}_{i}\,{\rm{d}}{V},\end{eqnarray} \tag{ 6 }$$

where ρ, p, and e are the comoving mass density, pressure, and internal energy density, respectively, and Γ is the Lorentz factor of the fluid. The subscript i labels the individual component, and the scalar field s_i represents the fraction of each component filling the local volume dV; within each cell, ${\sum }_{i=1}^{3}{s}_{i}=1$. This decomposition is accomplished by assigning three passive scalars to individual computational cells. The jet material is injected with $s=(1,0,0)$, the merger cloud material initially has $s=(0,1,0)$, and the ISM has $s=(0,0,1)$. As the simulation evolves, individual cells generally acquire some of each component due to mixing at the grid scale. To obtain ${E}_{k,i}$ and ${E}_{t,i}$ for each component i, we integrate Equations (6) over the volume.

Zoom In Zoom Out Reset image size

The top panel of Figure 3 displays the time evolution of the kinetic and thermal energies in these three components. At the very beginning of jet propagation, the kinetic energy of the jet increases as it accelerates by expending its thermal energy supply. We also observe that simultaneously, the thermal energy content of the cloud material increases. This is the result of ${PdV}$ work, as well as shock heating done by the jet on the cloud material as it drills through. Around $100\,{\rm{s}}$, the jet reaches the outskirts of the merger cloud (see Figure 1), its kinetic energy saturates, and it stops performing work on the cloud material. The jet continues to cool adiabatically (see the dotted lines showing thermal energy in Figure 3) as it propagates into the circumburst environment.

The BNS merger event responsible for GW170817 occurred in the outskirts of an elliptical galaxy (Blanchard et al. 2017; Levan et al. 2017). Low ISM densities are not unusual in such environments, and therefore we have adopted values in the range $n={10}^{-4}\mbox{--}{10}^{-5}\,{\mathrm{cm}}^{-3}$ for the circumburst number density, which is assumed to be a constant for our discussion. In the comoving frame of the relativistic shell, the upstream ISM particles stream inward with Lorentz factor Γ. When an ISM particle crosses the shock front, the direction of its velocity becomes random after multiple collisionless interactions. In the lab frame, the average energy of each downstream ISM proton is ${{\rm{\Gamma }}}^{2}{m}_{p}{c}^{2}$. Detailed studies of jet dynamics and radiation have been covered in GRB reviews (e.g., Piran 1999; Mészáros 2006; Nakar 2007; Berger 2014; Kumar & Zhang 2015). During the coasting phase of the relativistic jet, its bulk Lorentz factor does not change substantially. However, it performs work on the ISM while at the same time accumulating mass. The total energy of the swept-up ISM is given by

$$\begin{eqnarray}&&{E}_{\mathrm{iso},\mathrm{ism}}\approx \displaystyle \frac{4\pi }{3}{{nR}}^{3}{{\rm{\Gamma }}}^{2}\propto {t}^{3}.\end{eqnarray} \tag{ 7 }$$

In Figure 3, the energy of the ISM is shown to be increasing throughout the coasting phase ∝t³, in agreement with Equation (7). ${E}_{\mathrm{iso},\mathrm{ism}}$ becomes comparable to the energy of the jet at lab time ${t}_{\mathrm{lab}}\sim 2.5\times {10}^{8}\,{\rm{s}}$. This is ∼3 times longer than the predicted deceleration time ${R}_{d}/{\rm{c}}$, according to the estimate of Kumar & Zhang (2015),

$$\begin{eqnarray}{R}_{d} & = & {(3{E}_{\mathrm{iso}}/4\pi {{nm}}_{p}{c}^{2}{{\rm{\Gamma }}}^{2})}^{1/3}\\ & \approx & {10}^{17}{E}_{\mathrm{iso},53}^{1/3}{n}^{-1/3}{{\rm{\Gamma }}}_{2}^{-2/3}\,\mathrm{cm},\end{eqnarray} \tag{ 8 }$$

which yields a deceleration time of $7\times {10}^{7}\,{\rm{s}}$ for our parameters.

The jet's transition from the coasting phase to the deceleration phase is accompanied by the formation of a strong forward shock, which then propagates into the ISM. A weaker reverse shock, which propagates into the jet ejecta, also forms. The thermal energy of every component increases at the lab time ${t}_{\mathrm{lab}}\gt {10}^{7}{\rm{s}}$. Throughout the deceleration phase, the bulk Lorentz factor decays as ${\rm{\Gamma }}\propto {R}^{-3/2}\propto {t}^{-3/2}$ (Blandford & McKee 1976; Kobayashi et al. 1999). In the bottom panel of Figure 3, the on-axis energy-averaged Lorentz factor ${{\rm{\Gamma }}}_{\mathrm{avg}}$ is shown to decay roughly as $\propto {t}^{-3/2}$, in agreement with the analytical estimate.

Figure 4 depicts the evolution of radial profiles of three physical variables at a representative off-axis polar angle $\theta =0.4$. After the acceleration phase, a highly relativistic thin shell is formed. The radial profile of the density follows a power-law decay with respect to the dynamical radius. At the shock front, the thermal energy density of the fluid is significant. Both the inner and outer boundaries of the simulation domain track the radial movement of BNS merger ejecta over the entire duration of the simulations.

Zoom In Zoom Out Reset image size

Synchrotron emission using the model of Sari et al. (1998) can be directly calculated from multidimensional hydrodynamical simulation data (e.g., van Eerten et al. 2010; De Colle et al. 2012). The main parameters determining the synchrotron radiation from the forward shock are the fraction ${\epsilon }_{B}$ of post-shock energy residing in magnetic fields, and the fraction ${\epsilon }_{e}$ in nonthermal electrons. We further adopt the convention that ${\xi }_{e}$ is the fraction of the electrons sharing the electron internal energy ${\epsilon }_{e}e$, and that the energy distribution of the relativistic nonthermal electrons is given by ${dN}/d\gamma \propto {\gamma }^{-p}$. We assume that ${\xi }_{e}=1$, and the electron spectral index p is taken as a free parameter. We perform simulations of the successful structured jet propagating in low-density environments with two different values for the ISM density, $n={10}^{-4},{10}^{-5}\,{\mathrm{cm}}^{-2}$. By varying the value of the observer viewing angle ${\theta }_{\mathrm{obs}}$ and the microphysical parameters (${\epsilon }_{e},{\epsilon }_{B},p$), we obtain two sets of off-axis light curves that match the broadband afterglow observations of GRB170817A.

The results of these fits are shown in Figure 5 (see also Margutti et al. 2018). Here we present light curves calculated from the structured jet simulation, contrasted with semianalytical light curves computed using BOXFIT (van Eerten et al. 2012) and a simpler top-hat jet profile. The top-hat profile has the same isotropic equivalent energy, ${E}_{\mathrm{iso},\mathrm{avg}}\,=6\times {10}^{52}\,\mathrm{erg}$, as the self-consistently simulated jet, and we adopt an opening angle of ${\theta }_{c}=0.15$ taken from modeling the simulated jet according to Equation (4). Given the same radiation parameters, the light curves calculated from each model peak at roughly the same time and exhibit similar peak fluxes. However, the early part of the afterglow light curve differs significantly between these two models. In particular, the off-axis light curve from the structured jet brightens earlier than the top-hat jet. The slope of the late decaying light curve from these two models is similar.

Zoom In Zoom Out Reset image size

The late appearance of the X-ray and radio emission completely rules out any on-axis ultrarelativistic jet models. Indeed, if a relativistic top-hat jet had been pointed away from us, the afterglow emission would have been first detected at a later time, when the emission from the decelerated jet entered our line of sight. The rising light curve from the structured jet is robustly shallower than that of off-axis top-hat models (Mooley et al. 2018), and is thus detectable at earlier times. The off-axis light curves from the structured jet naturally explain the GRB170817A afterglow emission.

A wide-engine scenario has been invoked in previous studies (Kasliwal et al. 2017; Gottlieb et al. 2018b; Nakar et al. 2018). During its propagation, the jet is slowed by its interaction with the merger ejecta, and the interaction eventually drives a quasi-spherical mildly relativistic outflow (see Figure 6(a)).

In our 2D jet simulations, the relativistic shell was well resolved during its propagation inside and outside the merger ejecta, via an AMR scheme. We find that by resolving the relativistic shell, the ejecta material that accumulates on top of the jet head is able to get pushed aside in the narrow-engine model. No strong "plug" instability effect has been observed (Lazzati et al. 2010; Mizuta & Ioka 2013; Gottlieb et al. 2018b). For the wide-engine model, the jet engine, with a large opening angle and small initial momentum, does not have enough power to penetrate the heavy ejecta material. It diverges halfway through and gets deflected to a large polar angle (see Figure 6(b)).

An angular structure is also formed in the wide-engine scenario, as shown qualitatively in Figure 6 and quantitatively in Figure 7. The energy angular distribution could again be fitted by a quasi-Gaussian model with an opening angle ${\theta }_{c}=0.33$, larger than the opening angle in the narrow-engine scenario. Furthermore, the wide jet is found to have a lower peak isotropic equivalent energy ${E}_{\mathrm{iso},\mathrm{peak}}=4\pi {\epsilon }_{0}$ $\,\approx \,9\,\times {10}^{51}\,\mathrm{erg}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Whereas in the narrow-engine model, roughly 20% of the jet energy is deposited in the merger cloud, we find that number to be $\sim 80 \%$ in the wide-engine scenario. This is revealed in the different kinetic energy these two components end up having after acceleration (see the top panel in Figure 8). The bottom panel in Figure 8 shows the time evolution of the energy-averaged Lorentz factor (Equation (5)) as a function of the polar angle. ${{\rm{\Gamma }}}_{\mathrm{avg}}$ ranges from 2 to 10 between polar angles of 0.0 and 0.4 (also see Figure 7).

Zoom In Zoom Out Reset image size

In the literature (e.g., Mooley et al. 2018), the stratified quasi-spherical explosion model utilizes an outflow profile: $E{(\gt {\rm{\Gamma }}\beta )\propto ({\rm{\Gamma }}\beta )}^{-\alpha }$. The energy power-law index value $\alpha =5$ has been found to match early observations ($\lt 100\,\mathrm{days}$). In the left column of Figure 9, we show the cumulative distribution of energy $E(\gt {\rm{\Gamma }}\beta )$ as a function of four velocity. For the wide-engine model, we find that $E(\gt {\rm{\Gamma }}\beta )$ is not well characterized by a single power law. Rather, α increases from roughly 0.3 on the low-velocity end toward 1 at ${\rm{\Gamma }}\beta \sim 10$ (qualitatively similar to results of Hotokezaka et al. 2018). This is in contrast with the narrow-engine model we explored in Section 3, where a significant fraction of the energy was seen to reside at high Lorentz factor.

The right column of Figure 9 shows the four-velocity distribution histogram of the total energy and the thermal energy at ${t}_{\mathrm{lab}}={10}^{8}\,{\rm{s}}$. A large amount of shock-generated thermal energy resides in the ultrarelativistic shell (${\rm{\Gamma }}\gt 10$) for the narrow-engine model. The shock-generated thermal energy of the cloud and the jet engine material has higher Lorentz factor compared with the thermal energy of the shock-heated ISM.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 10, we show broadband afterglow light curves computed from the wide-engine model, compared with observations at the radio, optical, and X-ray frequencies. For this model, we were only able to obtain a successful fit with a very low external density of $n={10}^{-5}\,{\mathrm{cm}}^{-3}$.

Off-axis light curves from the narrow-engine model and on-axis light curves from the wide-engine model are able to match the rising light curve observed in the first ∼100 days of GRB170817A. The rising light curve component in all of these cases is produced by stratification. In the narrow-engine model, angular stratification is of importance. The high latitude (here defined as near the axis) relativistic material decelerates and adds flux to the rising light curve (see Section 5) without producing a sudden brightening, even when the jet core decelerates and comes into view for off-axis observers. If an angularly structured jet is responsible for GRB170817A, the jet core is probably already being observed. In contrast, for the wide-engine model, which produces a quasi-isotropic explosion, the outflow is radially stratified. The slower materials catch up with the decelerating blast wave, driving the rising radiation. In both cases, the light curve comes from the mildly relativistic material (Nakar & Piran 2018; see discussion in Section 5). Both models predict that the afterglow light curve will decay ∼200 days after the merger and share roughly the same decay pattern.

In Figure 12, we display the X-ray synchrotron emission light curve calculated from the narrow-engine simulation with ISM density $n={10}^{-4}\,{\mathrm{cm}}^{-3}$. The light curve covers seven orders of magnitude in observer time, starting from 1 minute and extending to $\sim 30\,\mathrm{years}$ after the BNS merger. The microphysical parameters ${\epsilon }_{e}$ (relativistic electron energy fraction), ${\epsilon }_{B}$ (magnetic energy fraction), and p (slope of the electron distribution) are set to standard values: ${\epsilon }_{e}=0.1,{\epsilon }_{B}\,=0.01,\mathrm{and}\,p=2.2$. In order to check the accuracy of our numerical light curves, we compare the peak time and peak flux to estimates from existing analytical models. The first model (Estimate_A) is based on an adiabatic double-sided top-hat jet with total kinetic energy E_k, an initial opening angle ${\theta }_{j}$, and a simple hydrodynamical evolution model (Nakar et al. 2002; Granot et al. 2017). The peak time of the off-axis afterglow light curve occurs when the bulk Lorentz factor of the top-hat jet drops to ${\rm{\Gamma }}=1/{\theta }_{\mathrm{obs}}$. The peak time and peak flux are given as

$$\begin{eqnarray}&&{t}_{\mathrm{peak}}({\theta }_{\mathrm{obs}})=0.7(1+z){\left(\displaystyle \frac{{E}_{{\rm{k}},51}}{{n}_{0}}\right)}^{1/3}{\left(\displaystyle \frac{{\theta }_{\mathrm{obs}}}{0.1}\right)}^{2}\ \mathrm{days},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}{F}_{{\nu }_{m}\lt \nu \lt {\nu }_{c}}^{\mathrm{peak}}(t) & = & 0.6\displaystyle \frac{{g}_{1}(p)}{{g}_{1}(2.2)}{(1+z)}^{(3-p)/2}{D}_{L28}^{-2}\\ & & \times \,{\epsilon }_{e,-1}^{p-1}{\epsilon }_{B,-2}^{\tfrac{p+1}{4}}{n}_{0}^{\tfrac{p+1}{4}}{E}_{{\rm{k}},50.7}{\nu }_{14.7}^{(1-p)/2}{\theta }_{\mathrm{obs},-1}^{-2p}\ \mathrm{mJy}.\end{eqnarray} \tag{ 14 }$$

In another model (Estimate_B), the projected surface area and solid angle of emission are taken into consideration (Lamb & Kobayashi 2017). The peak time and peak flux are given by

$$\begin{eqnarray}{t}_{\mathrm{peak}}({\theta }_{\mathrm{obs}}) & = & 195\left[\displaystyle \frac{(5+p){(7-p)}^{1/3}}{{(p-1)}^{4/3}}\right]{({\theta }_{\mathrm{obs}}-{\theta }_{j})}^{8/3},\\ & & \times \,{n}_{-1}^{-1/3}{E}_{{\rm{k}},52}^{1/3}\ \mathrm{days},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{f}_{{\nu }_{m}\lt \nu \lt {\nu }_{c}}^{\mathrm{peak}}=C(p)f({\theta }_{\mathrm{obs}},{\theta }_{j}){({\theta }_{\mathrm{obs}}-{\theta }_{j})}^{2(1-p)}{\nu }^{(1-p)/2}\\ &&\times \,{E}_{k}\ {n}^{(1+p)/4}{\epsilon }_{B}^{(1+p)/4}{\epsilon }_{e}^{p-1}{D}_{L}^{-2}\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{Hz}}^{-1},\end{eqnarray} \tag{ 16 }$$

where the expressions for ${g}_{1}(p),C(p),f({\theta }_{\mathrm{obs}},{\theta }_{j})$ are given in Granot et al. (2017) and Lamb & Kobayashi (2017), and ${E}_{k}={E}_{\mathrm{iso}}{\theta }_{j}^{2}/2$ is the jet kinetic energy (double-sided). Here, we model the ultrarelativistic core of the structured jet simply as a uniform top-hat jet with kinetic isotropic equivalent energy ${E}_{\mathrm{iso},52}\sim 6$, and jet half-opening angle ${\theta }_{j}\sim 0.15$. The ISM density is set to the value adopted in the simulation, $n={10}^{-4}\,{\mathrm{cm}}^{-3}$. As shown in Figure 12, the analytical estimate of the peak time and peak flux at different viewing angles from Estimate_B is in agreement with the calculated light curve.

Zoom In Zoom Out Reset image size

The on-axis X-ray light curve shown in Figure 12 (top left panel) displays three temporal power-law segments: (1) an early steep decay phase ${F}_{\nu }\propto {t}^{-\alpha }$, with temporal index ${\alpha }_{1}\sim 2.3$, (2) a shallow decay (plateau) phase with index ${\alpha }_{2}\sim 0$, and (3) a later decay phase with index ${\alpha }_{3}\sim 2$. These light curve components share similarities with the on-axis X-ray light curves for GRBs observed by Swift (Zhang et al. 2006; Kumar & Zhang 2015).

The off-axis light curves shown in Figure 12 (top right and bottom two panels) exhibit an early rapidly fading phase followed by a later rebrightening. Both on- and off-axis light curves have three common stages: an early declining afterglow, an intermediate transition phase (the rising part), and a late afterglow (the late declining part). The early declining emission mainly comes from the shock-heated cloud (for the on-axis light curve, it is the jet instead) and decays on a timescale of minutes to days depending on the viewing angle. The flux contribution from the external shock in the ISM steadily increases. Both the intermediate transition phase and the late afterglow light curve come from the shock-heated ISM (see Figure 12).

Figure 13 shows the temporal and spatial decomposition of the computed light curve. First, we separate the entire simulation duration into three lab time periods: ${t}_{\mathrm{lab}}\gt {10}^{4}$ s, ${t}_{\mathrm{lab}}\gt {10}^{6}$ s, and ${t}_{\mathrm{lab}}\gt {10}^{7}$ s. Most of the early declining flux observed during approximately the first day of observer time is emitted during the lab frame time period ${10}^{4}\,{\rm{s}}\lt {t}_{\mathrm{lab}}\lt {10}^{6}\,{\rm{s}}$ for both on- and off-axis light curves. As seen in Figure 3, before lab time ${t}_{\mathrm{lab}}\sim {10}^{7}\,{\rm{s}}$, almost all of the thermal energy in the domain is in the jet engine and the merger cloud material. The early declining emission is due to the cooling of the post-shock jet engine and merger cloud material. At later times, the relativistic shell experiences strong forward and reverse shocks as it sweeps up and shocks the ISM. The internal energy from the shock-heated ISM then begins to play an important role in the synchrotron emission. The flux during the intermediate transition and late afterglow is mainly emitted during the lab time period ${t}_{\mathrm{lab}}\gt {10}^{7}\,{\rm{s}}$, consistent with the dynamical evolution of the thermal energy. The turning point between the early declining and intermediate transition phases depends on the Lorentz factor of the emitting shell. For on-axis observers, photons radiated at lab frame time ${t}_{\mathrm{lab}}$ and lab frame position ${\boldsymbol{r}}$ will reach an observer at observer time ${t}_{\mathrm{obs}}$ (see, e.g., Piran 1999; Mészáros 2006):

$$\begin{eqnarray}&&{t}_{\mathrm{obs}}=(1+z)({t}_{\mathrm{lab}}-{\boldsymbol{r}}\cdot {\boldsymbol{n}}/c)\sim (1+z){t}_{\mathrm{lab}}/2{{\rm{\Gamma }}}^{2}.\end{eqnarray} \tag{ 17 }$$

Here, ${\boldsymbol{n}}$ is a unit vector pointing in the direction toward the observer. Along the jet propagation direction, the bulk Lorentz factor of the on-axis relativistic shell is ∼100. A photon emitted from the shell at lab time ${t}_{\mathrm{lab}}\sim {10}^{7}\,{\rm{s}}$ will thus be received by on-axis observers at observer time ${t}_{\mathrm{obs}}\,\sim 8\,\mathrm{minutes}$. This roughly determines the turning point between the initial steep decay and the shallow decay phase of the on-axis light curve shown in Figures 12 and 13.

Zoom In Zoom Out Reset image size

Previous studies suggest that the initial steep decay phase of the on-axis light curve is linked to the tail of the GRB prompt emission and has internal shock origin (e.g., Barthelmy et al. 2005; Duffell & MacFadyen 2015). Our result supports this interpretation. We refer the reader to Section 6 for further discussion about the early declining light curve. The shallow decay phase has been previously interpreted in the context of a refreshed shock model (Rees & Mészáros 1998; Sari & Mészáros 2000; Zhang et al. 2006). Based on our simulation results, we find that the duration of the shallow decay phase depends on the initial bulk Lorentz factor and the isotropic equivalent energy of the relativistic jet. It also depends on the ambient density. The typical timescale of the plateau phase observed for classical GRBs is ${10}^{3}\sim {10}^{4}\ s$ (Kumar & Zhang 2015). The BNS case considered here, an energetic jet propagating in a very-low-density environment, results in a duration longer than this timescale.

For the off-axis light curve, the early rapidly fading and later rebrightening behavior distinguishes it from the on-axis light curve. In Figure 13, we divide the simulation domain into angular regions and calculate the flux contribution from each one of them. Off-axis observers will first detect radiation from the part of the outflow that is moving toward them, i.e., in the direction of the observer's line of sight. As time goes on, the decelerated relativistic shell at higher latitudes (i.e., closer to the polar axis) contributes to the rebrightening light curve at lower latitudes, driving the flux level smoothly to greater values (e.g., Lazzati et al. 2018).

The early rebrightening portion of the light curve comes from the off-axis mildly relativistic material moving toward the observer along the line of sight and is essentially "on-axis" emission with respect to the observer. At ${\theta }_{\mathrm{obs}}=0.4$, the slope of the rebrightening light curve from the $0.2\lt \theta \lt 0.6$ region is moderate, at ∼1.3 (top right panel), while the slope of the rebrightening light curve from the $0\lt \theta \lt 0.2$ region is significantly larger, at ∼3. The difference in the slope value results from whether the light curve is observed "on-axis" or "off-axis" with respect to the line of sight. Observers located outside of the beaming cone of the relativistic shell ${\theta }_{\mathrm{obs}}\gt 1/{\rm{\Gamma }}$ will see an "off-axis" light curve. The observed "off-axis" light curve should rise faster than ${t}_{\mathrm{obs}}^{3}$ (Nakar & Piran 2018). For the GW170817 BNS merger event, the fact that the observed multiband light curve is much shallower, scaling as ${F}_{\nu }\propto {t}_{\mathrm{obs}}^{0.78}$, implies that "on-axis" emission was always observed for this event (Nakar & Piran 2018).

For the structured jet model, that "on-axis" emission comes from the mildly relativistic sheath at an off-axis angle ${\theta }_{\mathrm{obs}}\sim 20^\circ$. The energy-averaged Lorentz factor at this angle is around ${\rm{\Gamma }}\sim 3$ (Figure 2, lower right panel), in agreement with the analytical constraint, ${\rm{\Gamma }}\sim 1.5\mbox{--}7$, from Nakar & Piran (2018). When the central ultrarelativistic core decelerates and become "on-axis," the light curve stops increasing and smoothly turns over. The peak flux is determined by the central relativistic core. This is consistent with the peak time and peak flux estimates discussed in Section 5.1.

A similar rebrightening feature occurs in the observations of short GRBs (see e.g., Campana et al. 2006; Gao et al. 2015), long GRBs (e.g., Margutti et al. 2010), and X-ray flashes (e.g., Huang et al. 2004). The analysis of the off-axis light curve made here may provide an alternative interpretation for these rebrightening events.

The early declining light curves are computed with the synchrotron radiation model applied to the optically thin shock-heated merger ejecta. The early emission (hereafter a merger flash) decreases in time due to adiabatic cooling of the previously accelerated electrons. The flash is overtaken in all wave bands by rising synchrotron radiation from the external shock after roughly a day.

For the GW170817 BNS merger event, any early declining phase has been missed. The optical flux of the early synchrotron radiation is faint compared to the observed kilonova optical data (e.g., R-band). Early X-ray emission at several hours is below the instrument-detection limit of Chandra. These are shown in Figure 14, which displays the detection limits of various instruments, along with the observational data and two sets of fitting light curves.

Zoom In Zoom Out Reset image size

In X-ray, the late-XRT observations use a detection limit of $2\times {10}^{-14}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ (for a $10\,\mathrm{ks}$ exposure). For early-XRT, the detection limit is assumed to scale with the square root of the exposure time. For Chandra, we adopt a constant detection limit of ${10}^{-15}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$. In Figure 14, the X-ray detection limits have been converted to the flux limits in units of mJy assuming the default X-ray photon energy is $1\,\mathrm{keV}$. In the optical, R-band imaging detection limit for the Hubble Space Telescope's (HST) is set to 27. In the radio, the detection limit of the Very Large Array (VLA) is set to $10\,\mu \mathrm{Jy}$, assuming a $10\,{\rm{h}}$ reaction time.

The associated early X-ray light curve would be detectable by Swift XRT up to 30 minutes following the GRB prompt emission. The hard X-ray light curve ($15\,\mathrm{keV}$) becomes not or barely detectable by the Swift Burst Alert Telescope (BAT) and Fermi-GBM after one minute.⁴ However, under favorable conditions, the detection of the early declining afterglow in radio, optical, and X-ray at large off-axis angles may be possible for nearby BNS mergers.

As seen in Figures 5 and 10, both the narrow- and wide-engine models are capable of producing the late (${t}_{\mathrm{obs}}\gtrsim 1$ day) afterglow emission of GRB170817A. However, the nondetection of hard X-ray ($15\,\mathrm{keV}$) emission following GRB170817A by Fermi-GBM on the minute timescale may disfavor the wide-engine model because, even if seen at $\sim 20^\circ$ off axis, the quasi-isotropic explosion would have been detected by GBM at $\sim 15\,\mathrm{keV}$ for ∼10 minutes, as shown in the right plot of Figure 15. In the wide-jet scenario, which results in a quasi-isotropic explosion, fitting the observed late afterglow light curve with a lower-density ambient medium requires larger values of both ${\epsilon }_{B}$ and ${\epsilon }_{e}$. Also, unlike the off-axis structured jet, the quasi-isotropic outflow does not have flux contribution from the decelerating relativistic shell at high latitude. It needs a larger energy reservoir to drive the same level of afterglow flux. These two reasons would place the X-ray merger flash within the detection threshold of GBM, in disagreement with the $\lesssim 2\,{\rm{s}}$ duration of the detected short GRB signal. Nondetection of the X-ray merger flash at one minute by GBM also favors a higher-density ($n\sim {10}^{-4}\,{\mathrm{cm}}^{-3}$) over a lower-density ($n\sim {10}^{-5}\,{\mathrm{cm}}^{-3}$) ISM environment.

Zoom In Zoom Out Reset image size

GW170817 occurred in a part of the sky not accessible to Swift due to Earth occultation. However, had this event been accessible to the Swift satellite and XRT had slewed to its location within minutes, we show in the left plot of Figure 15 that XRT could have detected a declining merger flash at $\sim 1\,\mathrm{keV}$ lasting for ∼10 minutes.

Future BNS merger detections are expected to occur more frequently at larger distances, $\gtrsim 100\,\mathrm{Mpc}$. Had GW170817 occurred at that distance rather than at $40\,\mathrm{Mpc}$, the late X-ray rebrightening signal might not have been detected by Chandra. Therefore, it is important to understand the rapidly fading merger flash or what other types of electromagnetic transients might be detectable from BNS mergers at larger distances.

The emission model used to create the light curves in Figures 14 and 15 assumes the presence of synchrotron radiating nonthermal electrons. For it to be applicable, we require a mechanism to produce and sustain the nonthermal electron population, ${\epsilon }_{e}\lesssim 0.1$. We must also invoke the presence of magnetic energy at the level ${\epsilon }_{B}\sim {10}^{-2}$.

The presence of nonthermal electrons in the outflow can be supplied by shocks or magnetic reconnection. For example, the internal (subphotospheric) shock at $\sim 2\times {10}^{11}\,\mathrm{cm}$ might enable first-order Fermi acceleration that would supply nonthermal electrons. Another possibility is that particles are accelerated by reconnection of residual magnetic field in the plasma outflowing from the merger sight.

The presence of magnetic energy at the level ${\epsilon }_{B}\sim {10}^{-2}$ assumed in our synchrotron radiation modeling is justified by the presence of shocks, or by magnetic dynamo activity around the central engine. Indeed, subequipartition-level magnetic fields are expected to be produced downstream of the internal shock via the Weibel instability (although this depends on the uncertain kinetic physics of radiation-mediated shocks). Magnetic energy might also exist in the neutron star merger ejecta from either the premerger neutron star magnetic field or dynamo amplification during the merger itself (Zrake & MacFadyen 2013). Although magnetic energy density decreases as the merger ejecta expands, ${\epsilon }_{B}$ does not evolve significantly. This is because the energy density of the tangential magnetic field (B_ϕ and B_θ) in the coasting shocked cloud decreases like ${r}^{-2}$, while the gas internal energy decreases like ${r}^{-2\gamma }$ where the adiabatic index is $\gamma \geqslant 4/3$. Therefore, under expansion alone, ${\epsilon }_{B}$ either stays the same or marginally increases with radius.