Synchrotron Radiation from the Fast Tail of Dynamical Ejecta of Neutron Star Mergers

The discovery of the gravitational waves from a neutron star merger, GW170817, has opened a new era of multi-messenger astronomy (Abbott et al. 2017a, 2017b). The gravitational-wave signal was followed by multifrequency electromagnetic emission including a γ-ray pulse, uv/optical/IR macronova/kilonova, and long lasting nonthermal emission ranging from radio to X-rays. The macronova/kilonova observations show that the typical velocity of the ejecta is ∼0.1–0.3c and r-process elements of a mass of ∼0.05 M_⊙ are synthesized in the ejecta if the radioactive decay powers the emission (e.g., Andreoni et al. 2017; Arcavi et al. 2017; Coulter et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Evans et al. 2017; Kasliwal et al. 2017; Kilpatrick et al. 2017; Lipunov et al. 2017; Pian et al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Utsumi et al. 2017; Valenti et al. 2017, and see also the modelings, e.g., Kasen et al. 2017; Perego et al. 2017; Tanaka et al. 2017; Villar et al. 2017; Waxman et al. 2017). This mass estimate together with the rate of GW170817 supports early predictions (Lattimer & Schramm 1974; Eichler et al. 1989) that r-process elements in our Galaxy are predominately produced by mergers (e.g., Côté et al. 2017; Rosswog et al. 2017; Hotokezaka et al. 2018). However, the mechanism that ejects such a large amount of material still remains an open question (e.g., Shibata et al. 2017; Metzger et al. 2018).

X-ray and radio signals were discovered at nine and 16 days after the merger, respectively (Haggard et al. 2017; Hallinan et al. 2017; Troja et al. 2017). These signals are most likely explained by synchrotron radiation arising from the shock formed between the merger outflow and interstellar medium (ISM). Both off-axis and on-axis emission models are consistent with the observed data up to ∼30 days (e.g., Gottlieb et al. 2018b; Hallinan et al. 2017; Margutti et al. 2017; Troja et al. 2017. However, radio observations up to ∼100 days show that the flux density continued to rise as $\sim {t}^{0.8}$ (Mooley et al. 2018). X-ray and optical observations subsequently confirmed the increase in the flux up to 100 days and showed that the light curve peaked around ∼150 days (Ruan et al. 2017; Lyman et al. 2018; Margutti et al. 2018).

Recently, Nakar & Piran (2018) have shown that this moderately rising light curve implies that the emitting matter is moving toward us at the time of the emission; namely, the matter is moving within an angle θ < 1/Γ(Γis the Lorentz factor of the matter) from the line of sight toward us, otherwise the light curve would have risen much faster (F_ν ∝ t^α with α > 3). They have also shown a continuous energy injection into the blast wave. The isotropic equivalent energy increases like ${E}_{\mathrm{iso}}(\gt {\rm{\Gamma }})\propto {{\rm{\Gamma }}}^{-8\alpha +6(p-1)/(3-\alpha )}$ (where p is the electron's spectral index), otherwise the light curve would have declined.

This energy injection implies that the outflow must have a structure, either radial or angular or both (Lamb & Kobayashi 2017; Nakar & Piran 2018). Such a structure can arise naturally from the interaction of the jet with the surrounding ejecta that forms a cocoon (see e.g., Gottlieb et al. 2018a, 2018b; Kasliwal et al. 2017; Lazzati et al. 2017; Nakar & Piran 2018). This cocoon can arise in the case that the jet is choked or it emerged and generated a GRB pointing elsewhere.

Here, we explore a third possibility in which the observed synchrotron emission arises due to the fast tail of the dynamical ejecta and it has nothing to do with the question whether a jet existed or not and how it evolved. Therefore this model does not provide any clear connection between the nonthermal afterglow and the prompt γ-ray emission, GRB 170817A. In the dynamical ejecta scenario, suggested long ago by Nakar & Piran (2011) and elaborated by Piran et al. (2013), Hotokezaka & Piran (2015) and others, the synchrotron emission arises from the interaction of the mergers' dynamical ejecta with the surrounding ISM. These earlier studies focused on the late-time emission that would arise from the subrelativistic component of the outflow, but stressed that a very strong early signal is expected if a faster component exists.

Indeed, a high-velocity (mildly relativistic) tail of the ejecta can also explain the observed emission that followed GW170817 (Mooley et al. 2018). Such a high-velocity component was expected, but it was very difficult to estimate because only a very small amount of matter moves at these high velocities. The profile of the high-velocity tail was poorly known because of the complexity of hydrodynamics of mergers. Previous attempt to calculate this involved both analytic considerations (Kyutoku et al. 2014) and numerical simulations (Bauswein et al. 2013; Hotokezaka et al. 2013). The latter were limited by their resolution. Here we use the results of a new series of numerical simulation (Kiuchi et al. 2017) with a much higher spatial resolution that allows us to explore in details the velocity profile. These simulations find, in all cases considered, a light fast component with very steep energy profiles (see Figures 1–3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The structure of the paper is as follows. We show in Section 2 that recent high-resolution numerical relativity simulations reveal such high-velocity components and discuss their properties. In Section 3 we use these profiles to calculate the synchrotron radiation light curves arising from the shock between the ejecta and ISM and compare them with the observed data. Finally, in Section 4, we discuss the generic features and future expectations of the dynamical ejecta scenario and summarize our conclusions in Section 5.

A small fraction of the dynamical merger ejecta accelerates to mildly relativistic velocities when the shock formed between the colliding neutron stars emerges from the surface (Kyutoku et al. 2014). This component is the high-velocity tail of the dynamical ejecta. The flux of the synchrotron radiation arising from the ISM-ejecta shock can be significantly enhanced by even a small amount of a fast ejecta (Piran et al. 2013; Hotokezaka & Piran 2015). While fast ejecta have been found in previous numerical merger simulations (Bauswein et al. 2013; Hotokezaka et al. 2013; Metzger et al. 2015), the amount was not clear due to the low numerical resolution of those earlier simulations.

We begin by describing the ejecta profile obtained from recent high-resolution numerical relativity simulations (Kiuchi et al. 2017), whose the finest grid spacing, Δx, satisfies Δx ∼ 70 m. We discuss five models with a total mass of ≈2.7 M_⊙, a mass ratio of 1 or ≈0.8, and three different neutron star equations of state: models B, HB, and H (Kiuchi et al. 2017). These models are consistent with the mass and tidal deformability constraints of GW170817 (Abbott et al. 2017a). The parameters of these models and the resulting total kinetic energy and mass of the ejecta are listed in Table 1. Hereafter, we refer to equations of state that give small (large) radius neutron stars as "soft" ("stiff") equations of state.⁶ Here, we obtain the ejecta profile from snapshots of the simulations around 10 ms after the merger, corresponding a spatial scale of ∼1000 km (see, e.g., Hotokezaka et al. 2013 for the method). The radial velocity of the ejecta at each radius is converted to the one at infinity by using the gravitational potential of 2.7 M_⊙. Note that an artificial atmosphere.⁷

Table 1.  Parameters

Model	${m}_{1},{m}_{2}$ [M⊙]	${R}_{\mathrm{ns}}\,[\mathrm{km}]$	${M}_{\mathrm{ej}}\,[{M}_{\odot }]$	${E}_{\mathrm{ej}}$ [erg]
B	$1.35,1.35$	10.96	5 · 10−4	$7\cdot {10}^{49}$
HB	1.35, 1.35	11.61	3 · 10−3	1050
H	1.35, 1.35	12.27	3 · 10−3	9 · 1049
B uneq	1.51, 1.21	10.96	10−2	8 · 1050
HB uneq	1.51, 1.21	11.61	5 · 10−3	3 · 1050

Note. R_ns is the radius of a non-rotating cold neutron star of 1.35 M_⊙.

Download table as:  ASCIITypeset image

Note that, as in practically all other numerical relativistic merger simulations (Bauswein et al. 2013; Hotokezaka et al. 2013; Sekiguchi et al. 2015; Radice et al. 2016; Dietrich et al. 2017; Shibata et al. 2017) in all models, the total mass of dynamical ejecta is lower than the one needed to power the observed macronova/kilonova. This additional mass can be driven by other mass ejection mechanisms, such as winds from the surrounding accretion disk or from the hypermassive neutron star that have not been included in these merger simulations (see, e.g., Dessart et al. 2009; Fernández & Metzger 2013; Metzger & Fernández 2014; Perego et al. 2014; Fernández et al. 2015; Just et al. 2015; Fujibayashi et al. 2017; Siegel & Metzger 2017). However, the typical velocity of this wind ejecta is much slower than that of the dynamical ejecta, hence this ejecta is not relevant for the early synchrotron emission focused in this paper.

Figures 1–3 depict the kinetic energy distributions as a function of Γβ obtained from the simulations (see also Figures 6–8 for the ejecta mass distributions). Here, Γ is the Lorentz factor, and β is the velocity in units of the speed of light. The total kinetic energy of the ejecta is in the range from ∼7 · 10⁴⁹ to 8 · 10⁵⁰ erg. The ejecta of the softest model, B, is faster than that of models HB and H. This can be explained if the shock formed at the collision is stronger for more compact neutron stars, thus the ejected material is faster. For unequal mass cases, the ejecta are predominantly produced through the tidal disruption of the lower-mass neutron star, and it contains a larger amount of slow material and a smaller amount of fast material compared with the equal mass cases. Note also that the ejecta mass of model B is much smaller than the other models because the merging neutron stars promptly collapse into a black hole, thus the mass ejection occurs on only a short timescale.

The isotropic equivalent energy distributions at different polar angles are also shown in Figures 1–3, demonstrating the anisotropy of the ejecta. For the equal mass cases, the material ejected at 30°–60° has more kinetic energy. For the unequal mass cases, more material is ejected on the equatorial plane. Another important feature is that the ejecta moving faster than Γβ ∼ 1 is somewhat isotropic for all the cases. Therefore, the early light curve arising from this mildly relativistic tail is expected to be rather isotropic.

Also shown in Figures 1–3, for a comparison, is the quasi-spherical slow ejecta model used by Mooley et al. (2018):

$$\begin{eqnarray}&&E{(\gt {\rm{\Gamma }}\beta )=5\cdot {10}^{50}\mathrm{erg}({\rm{\Gamma }}\beta /0.4)}^{-5}\ (\mathrm{for}\,\beta \lt 0.8).\end{eqnarray} \tag{ 1 }$$

This distribution can fit the observed data up to ∼100 days after GW170817. Note that the nonthermal afterglow up to ∼100 days is produced by a component faster than β ≈ 0.6 (Mooley et al. 2018) and, of course, this profile cannot be extended to lower velocities β ≪ 0.4. The total kinetic energy distributions for model B is the closest to this distribution.

While we do not discuss any jet and ejecta interaction and the production mechanism of the prompt emission of GRB 170817A in this work, here we point out implications of the fast ejecta tail to short GRBs and in particular to GRB 170817A. Recently, Gottlieb et al. (2018b) proposed that a cocoon shock breakout can explain the observed γ-rays of GRB 170817A and found that a high-velocity component moving at velocities of ∼0.7–0.8c is needed to explain these observations. It is worth noting that the amount and the profile of the ejecta shown in Figures 1–3 are compatible with those needed to explain GRB 170817A according to this model. In the context of regular sGRB models, this fast-velocity ejecta leads to a longer delay than previously estimated between the GRB and the merge. A jet launched by the merger remnant must emerge from the ejecta surface that is moving at 0.7–0.8c to successfully produce a GRB. This velocity of the ejecta surface is faster by a factor of ∼2 than those adopted in the previous studies of the jet-ejecta interaction (e.g., Murguia-Berthier et al. 2014; Nagakura et al. 2014; Duffell et al. 2015; Gottlieb et al. 2018a). Therefore, the jet emergence time from the ejecta surface can be longer than those found in these previous works.

The interaction between the expanding merger ejecta and the surrounding ISM results in a shock in which particle acceleration and magnetic field amplification occur. As a result, this blast wave emits multiwavelength synchrotron radiation (Nakar & Piran 2011). Here, we numerically calculate the synchrotron radiation emitted by the accelerated electrons, as we are interested in mildly relativistic blast waves for which either limits, β ≪ 1, or Γ ≫ 1, cannot be used. To follow the evolution of the ejecta expanding in the ISM with a constant density, n, we solve the adiabatic radial expansion of the ejecta at each polar angle using the kinetic energy distributions given in the previous section. The synchrotron flux from the blast wave element at each solid angle at a given observer time is calculated assuming a power-law distribution of electrons, ${dN}/d{\gamma }_{e}\propto {\gamma }_{e}^{-p}$, and the standard equipartition parameters _e and _B, which reflect the conversion efficiency of the internal energy of the shocked ISM to the accelerated electrons' energy and magnetic energy respectively (e.g., Sari et al. 1998; Nakar & Piran 2011 and see Hotokezaka & Piran 2015 for details). We fix the viewing angle to be 30° (Abbott et al. 2017b) and set _e to be 0.1. Therefore, there are in total three free parameters, n, _B, and p, with which we fit the observed data in our afterglow modeling.

Figure 4 shows the calculated light curves for the different models as well as the observations of GW170817. The light curves match the observed data, for n = 0.008–0.04 cm⁻³, _B = 0.05–0.1, and p = 2.1. Note that even the higher values of the ISM density are consistent with the upper limit on the mean ISM density of NGC 4993,⁸ $n\lesssim 0.04\,{\mathrm{cm}}^{-3}$, inferred from the upper limit on the H i mass (Hallinan et al. 2017). For the equal mass cases, the radio and optical light curves continuously rise up to ∼100 days and then have a plateau phase lasting a few hundred days. This early component is produced by the tail of the ejecta faster than Γβ ∼ 0.6. Model H contains a lower kinetic energy at Γβ ∼ 1. As this fastest component dominates the early light curve, the flux density up to 20 days is lower than those of model B and HB.

Zoom In Zoom Out Reset image size

For unequal masses, the flux densities continuously rise until 1000 days because of the larger amount of material at low velocities. However, the slope of the calculated light curves is steeper than the observed one because the kinetic energy distribution of the these models declines more steeply at high velocities.

An important feature of the X-ray light curves is that they peak on a timescale of a few months to a year, and then their temporal evolution is different from that of the radio and optical light curves once the synchrotron cooling frequency becomes lower than 1 keV. For model H, in particular, the cooling frequency falls below 1 keV at early times ∼10 days because of the required relatively large values of the ISM density and _B. Note that the exact location of the cooling frequency can be higher, due to several uncertainties (see the discussion in the next section). However, it is generally expected that the cooling frequency crosses the X-ray band on a timescale of a few months to a year in this scenario.

Cooling frequency. An important feature of the synchrotron light curves obtained in Section 3 is the evolution of the synchrotron cooling frequency, which crosses the X-ray band on a timescale of a few months to a year. After this time, the X-ray flux density declines with time faster than the radio and optical flux densities. In the Newtonian limit (β ≪ 1), the cooling frequency at a given time, t, is approximately estimated as (e.g., Sari et al. 1998)

$$\begin{eqnarray}\begin{array}{rcl}{\nu }_{c} & \sim & 0.5\,\mathrm{keV}{\left(\displaystyle \frac{\beta }{0.6}\right)}^{-3}{\left(\displaystyle \frac{n}{0.01{\mathrm{cm}}^{-3}}\right)}^{-3/2}\\ & & \times \,{\left(\displaystyle \frac{{\epsilon }_{B}}{0.05}\right)}^{-3/2}{\left(\displaystyle \frac{t}{100\mathrm{days}}\right)}^{-2},\end{array}\end{eqnarray} \tag{ 2 }$$

or the cooling frequency for a given flux density is

$$\begin{eqnarray}&&{\nu }_{c}\propto {\beta }^{\tfrac{5p-6}{3}}{n}^{\tfrac{p-4}{6}}{\epsilon }_{B}^{\tfrac{p-8}{6}}{\epsilon }_{e}^{\tfrac{2(p+1)}{3}}{F}_{\nu }^{-\tfrac{2}{3}},\end{eqnarray} \tag{ 3 }$$

where F_ν is the flux density at a frequency in the range of ν_a, ν_m < ν < ν_c. Note that for a given flux density, the cooling frequency increases by decreasing n and _B and by increasing the velocity for the expected range of the electrons' power-law index, 2 < p < 3. In other words, it depends sensitively on the kinetic energy of the fast tail. The cooling frequency can be higher than those of our models if the kinetic energy at velocities faster than β ∼ 0.6 is only slightly larger. For instance, the kinetic energy distribution of Equation (1) results in ν_c ≳ 10 keV at 100 days (Mooley et al. 2018). Note also that the cooling frequency is much higher for the cocoon models (Lyman et al. 2018; Margutti et al. 2018; Mooley et al. 2018; Nakar & Piran 2018), which involve an outflow much faster than the fast tail of the dynamical ejecta (see, e.g., Gottlieb et al. 2018b; Lazzati et al. 2017; Nakar & Piran 2018).

The X-ray observations around 100 days show that the X-ray flux density and the photon index in the X-ray band are consistent with a single power-law spectrum from the radio to the X-ray bands (Ruan et al. 2017; D'Avanzo et al. 2018; Margutti et al. 2018; Mooley et al. 2018; Troja et al. 2018). This suggests that the cooling frequency is above the X-ray band, ruling out models with much lower cooling frequency, e.g., model H. However, it should be noted that the estimate of Equation (2) and the cooling frequency used in our light-curve modeling are correct only within an order of magnitude. In fact, Granot & Sari (2002) obtained the spectral breaks of the afterglows for the Blandford–McKee solution and found that the cooling frequency was higher by a factor of ≳3 than the simple order-of-magnitude estimate (Sari et al. 1998). Therefore, models B and HB for the equal mass cases can be considered as consistent with the current observed data. However, a generic expectation of the dynamical ejecta scenario is that the cooling frequency will cross the X-ray band on a timescale not much longer than a few months. Therefore, further multifrequency observations will easily confirm or rule out this model.

Note that model B results in a direct collapse to a black hole, i.e., a hypermassive neutron star does not form after the merger. The remnant is a hypermassive neutron star for model HB and H. Model B contradicts the kilonova/macronova observation because the total mass of the ejecta in model B is expected to be smaller than that inferred from the observed luminosity. Moreover, the early blue emission seems to require the existence of a hypermassive neutron star (e.g., Perego et al. 2017; Shibata et al. 2017). However, the condition for the formation of a hypermassive neutron star is sensitive to the binary's total mass and the maximal neutron star mass. Thus, we expect that a kinetic energy distribution similar to model B with a hypermassive neutron star remnant can be obtained with slightly different parameters of model B, e.g., a different maximal neutron star mass and/or a different binary total mass.

Late-time signal. In addition to the fast moving dynamical ejecta, there should be a main, subrelativistic ejecta component that have produced the UV/optical/IR macronova/kilonova signal. As already mentioned the dynamical ejecta masses of our models, 5 · 10⁻⁴–10⁻² M_⊙, are much smaller than those estimated from the macronova observations, $\sim 0.05\,{M}_{\odot }$ (e.g., Kasliwal et al. 2017). This component has the photospheric velocities of ∼0.1–0.3c and a kinetic energy of $\sim {10}^{51}$ erg, which is also larger than the kinetic energy of the dynamical ejecta in our dynamical ejecta models.⁹ Some other processes, not taken into account in the simulations used in this work, must be responsible for the ejection of this additional mass.

Considering now the total ejected mass as observed from the UV/optical/IR macronova, we estimate the peak time and flux density of this slow component (Nakar & Piran 2011):

$$\begin{eqnarray}&&{t}_{{\rm{p}}}\approx {10}^{4}\,\mathrm{day}{\left(\displaystyle \frac{E}{{10}^{51}\mathrm{erg}}\right)}^{\tfrac{1}{3}}{\left(\displaystyle \frac{n}{0.01{\mathrm{cm}}^{-3}}\right)}^{-\tfrac{1}{3}}{\left(\displaystyle \frac{\beta }{0.2}\right)}^{-\tfrac{5}{3}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\nu ,{\rm{p}}} & \approx & 2\,\mathrm{mJy}\,{\epsilon }_{B,-1}^{0.775}{\epsilon }_{e,-1}^{1.1}\left(\displaystyle \frac{E}{{10}^{51}\,\mathrm{erg}}\right)\\ & & \times \,{\left(\displaystyle \frac{n}{0.01{\mathrm{cm}}^{-3}}\right)}^{0.775}{\left(\displaystyle \frac{\beta }{0.2}\right)}^{1.75}{\left(\displaystyle \frac{\nu }{3\mathrm{GHz}}\right)}^{-0.55},\end{array}\end{eqnarray} \tag{ 5 }$$

where ${\epsilon }_{B,-1}$ and ${\epsilon }_{e,-1}$ are normalized by 0.1 and we set p = 2.1. Figure 5 depicts the flux densities from the dynamical ejecta and the macronova ejecta. Here, we assume that the macronova ejecta has a mass of $0.04\,{M}_{\odot }$, and a single velocity of $0.2c$ (top panel) and 0.1c (bottom panel), and set the ISM density and the same microphysics parameters as those used for the dynamical ejecta models (see Figure 4). The bumps around 10⁴–10⁵ days in the light curves are produced by the macronova component. The actual shape of the peak of the light curves is expected to be broader, due to the velocity structure of the macronova ejecta. Note that the cooling frequency of the late emission is around the uv/optical bands so that the estimate of Equation (5) is applicable up to the optical band and the X-ray flux density of the macronova component is much fainter than that estimated by Equation (5).

Zoom In Zoom Out Reset image size

This late-time emission is one of the notable differences between the dynamical ejecta and cocoon/structure jet scenarios. For the cocoon/structure jet scenario, in which the ISM density is much lower than that of the dynamical ejecta scenario, the flux density from the macronova ejecta is expected to be much fainter ≲50 μJy and the peak time is much longer, ≳100 years.

Angular size. The velocity of the fast tail that produces the afterglow at 100 days is ∼0.6c. Therefore, the angular size of the afterglow is $\lesssim 500\,\mu \mathrm{arcsec}(t/100\,\mathrm{day})$ for $\gt 100$ days. This is smaller by a factor a few compared with that expected in the cocoon and structured jet models. Recently, Mooley et al. (2018) found that the radio emission region of GW170817 exhibits a superluminal motion. This observation provided us direct evidence for the existence of a narrowly collimated jet and that the afterglow light curve around and after the peak is predominantly produced by the jet. Thus, the dynamical ejecta scenario is ruled out for explaining the observed light curve of GW170817.

The fast tail (Γβ ∼ 1) of the dynamical ejecta of binary neutron star mergers, calculated in recent high-resolution numerical simulations by Kiuchi et al. (2017) contains kinetic energy of 10⁴⁷–10⁴⁹ erg, depending on the neutron star equation of state and on the binary masses. Mergers with a softer neutron star equation of state, which gives smaller radius neutron stars, and with a mass ratio close to unity eject larger amounts of the fast tail. The fast tail has somewhat isotropic shape even for models in which the bulk of the material is large aspherical.

The interaction of this fast tail with the surrounding ISM produces a blast wave whose synchrotron emission (Nakar & Piran 2011) is observed as the radio to X-ray signals that followed GW170817. We compute this synchrotron emission arising from the ejecta profile obtained from the high-resolution numerical relativity simulations and compare them with the observed data of the nonthermal radiation that followed GW170817. We find that the multifrequency observed data can be reproduced well for the equal mass binary models with a relatively soft equation of state. In all cases, an ISM density of $\sim 0.01\,{\mathrm{cm}}^{-3}$ is required to obtain the observed flux level at ∼100 days after the merger. For unequal mass cases, the velocity gradient of the ejecta profile is steeper and the light curves rise more steeply than the observed one so that nearly equal mass mergers are favored in this scenario.

The dynamical ejecta scenario has three generic predictions.

• 1.  
  The cooling frequency crosses the X-ray band on a timescale of a few months to a year leading to a cooling break in the X-ray light curve.
• 2.  
  The outer ejecta velocity is ≲0.6c at t ≳ 100 days so that the angular size of the ejecta will be ≲500 μarcsec(t/100 day).
• 3.  
  The emission of the subrelativistic macronova/kilonova component of a mass of $\sim 0.05\,{M}_{\odot }$ and a velocity of ∼0.1c will be continuously observable with flux densities of 0.1–1 mJy on timescales of 10³–10⁵ days (Nakar & Piran 2011; Piran et al. 2013; Hotokezaka & Piran 2015).

A recent study of the very long baseline interferometry observation of GW170817 ruled out the dynamical ejecta scenario for the afterglow of GW170817 (Mooley et al. 2018). While this model does not explain the afterglow observed in GW170817 up to 300 days, the model may be relevant for the afterglows of future merger events.

We thank Paolo D'Avanzo, Gregg Hallinan, Raffaella Margutti, Kunal Mooley, David Radice, and Eli Waxman for useful discussions. We are also grateful to the anonymous referee for improving the manuscript. K.H. was supported by Lyman Spitzer, Jr. fellowship at Department of Astrophysical sciences at Princeton University. T.P. was partially supported by an advanced ERC grant TReX and by a grant from the Templeton foundation. T.P. and E.N. were supported by the I-Core center of excellence of the Israeli Science Foundation. Numerical simulations were performed on K computer at AICS (project numbers hp160211, hp170230, and hp170313), on Cray XC30 at cfca of National Astronomical Observatory of Japan, FX10 at Information Technology Center of the University of Tokyo, and on Cray XC40 at Yukawa Institute for Theoretical Physics, Kyoto University. This work was supported by Grant- in-Aid for Scientific Research (16H02183, 17H06361, and 15K05077) of JSPS and by a post-K computer project (Priority issue No. 9) of Japanese MEXT.

Figures 6–8 show the ejecta mass distributions as a function of velocity and isotropic mass distribution at different polar angles. Figure 9 shows the dependence of the ejecta distribution on the numerical resolution. The mass in the fast-velocity tail increases with the numerical resolution. Therefore, we consider the result of the simulations with the highest resolution as our current best estimate of the fast-velocity tail profile.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size