Relativistic envelopes and gamma-rays from neutron star mergers

TL;DR

This paper proposes a model where neutron star mergers eject an ultra-relativistic envelope, which explains the gamma-ray burst observed in GW170817 through a delayed shock breakout mechanism involving a magnetar and subsequent jet activity.

Contribution

The paper introduces a novel model of gamma-ray burst emission from neutron star mergers involving an ultra-relativistic envelope and shock breakout, explaining observed burst properties.

Findings

The model reproduces the luminosity and duration of GW170817's gamma-ray burst.

The Lorentz factor profile of the envelope matches observational constraints.

The proposed shock breakout mechanism provides a broad-angle gamma-ray emission explanation.

Abstract

We suggest that neutron star mergers eject an ultra-relativistic envelope of mass $m\sim 10^{-7}M_\odot$, which helps explain the gamma-ray burst from GW170817. One ejection mechanism is the ablation of the neutron star surface by the burst of neutrinos in the first $30\,\mu$s of the merger. Another, more efficient, mechanism for inflating the ultra-relativistic envelope is an internal shock in the massive ejecta from the merger. A strong shock is expected if the merger product is a magnetar, which emits a centrifugally accelerated wind. The shock propagates outward through the ejecta and accelerates in its outer layers at radii $r\sim 10^9-10^{10}$cm, launching an ultra-relativistic opaque envelope filled with $\sim 10^4$ photons per nucleon. The Lorentz factor profile of the envelope rises outward and determines its homologous expansion, which adiabatically cools the trapped photons. Once the magnetar loses its differential rotation and collapses into a black hole, a powerful jet forms. It drives a blast wave into the envelope, chasing its outer layers and eventually catching up with the envelope photosphere at $r\sim 10^{12}$cm. The ultra-relativistic photospheric breakout of the delayed blast wave emits a gamma-ray burst in a broad solid angle around the merger axis. This model explains the gamma-ray pulse from merger GW170817 with luminosity $L_\gamma\sim 10^{47}$erg/s, duration $\Delta t_{\rm obs}\sim 0.5$s, and characteristic photon energy $\sim 100$keV. The blast wave Lorentz factor at the envelope photosphere is consistent with $\Gamma\geq 5$ that we derive from the observed light curve of the burst. We suggest future tests of the model.

Equations

$${\tau }_{\mathrm{T}}=\frac{\zeta {\sigma }_{\mathrm{T}}{n}_{p}r}{\Gamma }=\frac{\zeta {\sigma }_{\mathrm{T}}L{Y}_e}{4\pi r{m}_p{c}^3{\Gamma }^3},$$

$$\Gamma >{\Gamma }_{\star }\approx 4.9r_{12}^{-1/3}{\left(\frac{\zeta Y_{e}L}{L_{\gamma }}\right)}^{1/3}{L}_{\gamma ,47}^{1/3}.$$

$$r\text{raise}2.15277pt<to0.0pt\text{hss}\text{lower}2.15277pt\sim 2{\Gamma }^{2}c\Delta t_{\mathrm{obs}}=3\times 10^{10}{\Gamma }^2\left(\frac{\Delta t_{\mathrm{obs}}}{0.5\mathrm{s}}\right)\mathrm{cm}.$$

$$\Gamma >5.7{\left(\frac{L_{\gamma }}{1.6\times 10^{47}erg/s}\right)}^{1/5}{\left(\frac{\Delta t_{\mathrm{obs}}}{0.5\mathrm{s}}\right)}^{-1/5}{\left(\frac{\zeta Y_{e}L}{L_{\gamma }}\right)}^{1/5}.$$

$$\hat{E}\approx \frac{v_{\star }^2}{2c^2}.$$

$$P$$

$$k{T}_{\max }\equiv \frac{2}{3}\hat{E}{m}_n{c}^2\approx 20\left(\frac{\hat{E}}{0.03}\right)\mathrm{MeV}(\mathrm{no}\mathrm{cooling}),$$

$$F\sim \frac{\rho v_{\star }^3}{2}(\mathrm{efficient}\mathrm{cooling}).$$

$$kT\approx 15{\rho }_{12}^{1/4}{v}_{\star ,10}^{3/4}\mathrm{MeV}(\mathrm{efficient}\mathrm{cooling}).$$

$${ϵ}_{\nu }\approx 3\Theta \equiv 3\frac{kT}{m_e{c}^2}(\mathrm{thermalized}).$$

$$l_{\nu }=\frac{1}{{\kappa }_{\nu }\rho }\approx 3.6\times 10^5{\rho }_{11}^{-1}{\left(\frac{\Theta }{30}\right)}^{-2}\mathrm{cm}.$$

$$w\equiv \frac{U+P}{\rho c^2}$$

$${\sigma }_{\nu \bar{\nu }}\approx 0.2{\sigma }_0\frac{c^4({\mathbf{p}}_{\nu }\cdot {\mathbf{p}}_{\bar{\nu }}{)}^2}{(E_{\nu }{E}_{\bar{\nu }}{)}^2}=0.2{\sigma }_0{ϵ}_{\nu }{ϵ}_{\bar{\nu }}(1-\cos \delta {)}^2,$$

$$\dot{Q}\sim {\sigma }_{\nu \bar{\nu }}c{n}_{\nu }{n}_{\bar{\nu }}(E_{\nu }+E_{\bar{\nu }})\sim \frac{{\sigma }_{\nu \bar{\nu }}}{c}\frac{F_{\nu }}{E_{\nu }}\frac{F_{\bar{\nu }}}{E_{\bar{\nu }}}(E_{\nu }+E_{\bar{\nu }}).$$

$$\dot{Q}\sim {\sigma }_0\overline{(1-\cos \delta {)}^2}\frac{{\Theta }_{\star }{F}^2}{m_e{c}^3}\sim 10^{35}{F}_{41}^{9/4}{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{cm}}^{-3},$$

$$F_h=\int \dot{Q}ds\sim H\dot{Q},$$

$$U\sim \dot{Q}{t}_a\sim 10^{30}{F}_{41}^{9/4}{t}_{-5}{\mathrm{erg}\mathrm{cm}}^{-3}$$

$$m\sim bV\frac{\dot{Q}{t}_a}{c^2}\sim 10^{-7}{M}_{\odot }{b}_1{V}_{16}{F}_{41}^{9/4}.$$

$$\dot{Q}(t,s)={\dot{Q}}_0(t){\left(\frac{a+s}{a}\right)}^{-4}{\left[1-\frac{s}{(s^2+a^2{)}^{1/2}}\right]}^2.$$

$${\dot{Q}}_0(t)=10^{35}{\left(\frac{t}{t_0}\right)}^{\zeta }{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{cm}}^{-3},t<t_0,$$

$$\mathcal{E}\sim \frac{M_{\mathrm{ej}}(\delta v{)}^2}{2}\approx 10^{51}\left(\frac{M_{\mathrm{ej}}}{10^{-2}{M}_{\odot }}\right){\left(\frac{\delta v}{10^{10}cm/s}\right)}^2\mathrm{erg}.$$

$$\tilde{\mathcal{E}}\equiv \frac{\mathcal{E}}{M_{\mathrm{ej}}{c}^2}\sim 0.01-0.1,$$

$${\gamma }_s{\beta }_s\approx A{\tilde{\mathcal{E}}}^{1/2}{\left(\frac{M_{\mathrm{ej}}}{\rho r^3}\right)}^{\alpha },$$

$$m_1\equiv m(\beta \gamma >1)\sim 10^{-6}{M}_{\odot }{\left(\frac{\tilde{\mathcal{E}}}{0.03}\right)}^4\left(\frac{M_{\mathrm{ej}}}{0.01M_{\odot }}\right),$$

$$\frac{m(\beta \gamma >3)}{m(\beta \gamma >1)}\sim 3\times 10^{-3}.$$

$$\gamma \beta \approx {\left(\frac{m}{m_1}\right)}^{-\psi }.$$

$$\frac{n_{\gamma }}{n_b}\approx \frac{(3/4)w_0{m}_p{c}^2}{2.7k{T}_0}\sim 10v_{10}^{7/4}{\dot{M}}_{31}^{-1/4}(\mathrm{without}\mathrm{shock}).$$

$$\frac{n_{\gamma }}{n_b}\sim 2Y_e\frac{7\mathrm{MeV}}{2.7k{T}_{\mathrm{rec}}}\sim 10,$$

$$\rho (m,t)\approx \frac{\xi m}{4\pi r^3(t)},\xi =\frac{r}{\Delta r}\text{raise}2.15277pt>to0.0pt\text{hss}\text{lower}2.15277pt\sim 1.$$

$$U(m)\sim 2{\gamma }_s^2{v}_s^2(m)\rho (m)\sim \frac{\tilde{\mathcal{E}}{c}^2}{r^3}{\left(\frac{\xi m}{4\pi }\right)}^{0.6}{M}_{\mathrm{ej}}^{0.4}.$$

Full text

1. Introduction

1.1. Ejecta from neutron star merger GW 170817

The recent detection of gravitational waves from neutron star merger GW 170817 and its electromagnetic counterpart opens a new window for the studies of neutron stars, cosmological gamma-ray bursts (GRBs), and the origin of heavy nuclei (Abbott et al., 2017a, b; Goldstein et al., 2017; Savchenko et al., 2017; Coulter et al., 2017; Evans et al., 2017; Soares-Santos et al., 2017). The electromagnetic radiation was emitted by ejecta from the merger, viewed at an angle of $\theta\sim 20^{\circ}-30^{\circ}$ from the rotation axis of the binary. The viewing angle and the distance to the merger $d\sim 40$ Mpc are both estimated from the observed gravitational wave signal, and its host galaxy was found at $d\approx 40$ Mpc.

The gamma-ray counterpart, GRB 170817A, had luminosity $L\sim 10^{47}$ erg/s and was emitted with a delay of $\sim 1.7$ s following the merger. It could be powered by a delayed jet from the central remnant when it breaks out from a massive cloud around the remnant (Kasliwal et al., 2017; Gottlieb et al., 2018; Bromberg et al., 2018; Pozanenko et al., 2018).

A day later, the cloud of expanding ejecta with mass $M_{\rm ej}\sim 5\times 10^{-2}M_{\odot}$ and speed $v/c\sim 0.1-0.3$ emitted optical radiation with luminosity $L\sim 10^{41}$ erg/s — the “kilonova.” Its light curve was consistent with being powered by the decay of r-process nuclei (Kasen et al., 2017; Drout et al., 2017; Tanaka et al., 2017; Tanvir et al., 2017).

At yet later times (weeks), X-ray and radio afterglow was observed (Troja et al., 2017; Margutti et al., 2017; Hallinan et al., 2017). The unusually late rising afterglow was proposed to result from deceleration of quasi-isotropic, moderately relativistic ejecta from the merger in an external medium (Nakar et al., 2018). It was also found consistent with a decelerating ultra-relativistic narrow jet with energy $E_{\rm jet}\sim 10^{50}$ erg (Lazzati et al., 2018; Granot et al., 2018; Lamb et al., 2018; Xie et al., 2018), and further evidence for a collimated jet came from radio imaging (Mooley et al., 2018b, a). The jet is initially invisible to off-axis observers, because of its strong collimation and Doppler-beaming with a very high Lorentz factor. Its emission comes into view after significant deceleration, long after the merger.

1.2. The puzzling GRB 170817A

Generally, emission from GRB jets was expected to be weak and soft when viewed off-axis (e.g. Lazzati et al. 2017a, b). By contrast, the gamma-ray pulse of GRB 170817A is not soft; its spectrum peaks at $100-200$ keV (Goldstein et al. 2017), and an even harder spectrum was detected in a short time interval of \mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}0.1 s (Veres et al., 2018). Furthermore, GRB jets display a strong correlation between the apparent brightness and hardness of their emission; GRB 170817A is far off this correlation (Pozanenko et al. 2018).

In addition, the simple light curve of GRB 170817A favors a single emission event, such as a blast wave from a jet, rather than variable internal dissipation typical for GRB jets. The canonical GRBs viewed on-axis are extremely bright and have diverse (usually multi-peaked) light curves (Nakar, 2007). The off-axis GRB 170817A is dominated by a single weak pulse of width $\Delta t_{\rm obs}\sim 0.5$ s, smaller than but comparable with its delay $t_{\rm obs}\sim 1.7$ s. Goldstein et al. (2017) also reported an unusual transition from the gamma-ray pulse to a quasi-blackbody X-ray tail, although the tail has a low signal-to-noise ratio and its detailed spectral shape is uncertain.

A natural mechanism for a single, hard, gamma-ray pulse followed by a soft thermal tail could be the breakout of a shock wave from the massive cloud around the merger. Following previous theoretical calculations (e.g. Nakar & Sari 2012), the shock breakout model proposed for GRB 170817A (Kasliwal et al., 2017; Gottlieb et al., 2018; Bromberg et al., 2018) posits the shock temperature $kT\approx 50$ keV. It predicts the observed (Doppler-shifted) spectral peak at $3kT\Gamma\sim 150$ keV if the Lorentz factor of the shock-heated plasma is $\Gamma\sim 1$. The low $\Gamma$ is, however, in conflict with observations. In Section 2 we show that the light curve of GRB 170817A requires the gamma-ray source to have \Gamma\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}5.

The Lorentz factor $\Gamma\sim 5$ could be consistent with an off-axis component (“cocoon”) of an ultra-relativistic jet after its breakout from the massive cloud. However, this outflow is not expected to emit gamma-rays. Simulations of jet breakout and cocoon expansion show that the off-axis outflow is heated by internal shocks too early, before it becomes transparent to radiation, and this leads to reprocessed X-ray emission rather than gamma-ray emission (Lazzati et al., 2017a, b).

1.3. This paper

After estimating the lower limit on $\Gamma$ in GRB 170817A in Section 2, we turn to the theory of relativistic ejecta from neutron star mergers. We find that, before jets are launched and the GRB is emitted, the merger is likely to eject an ultra-relativistic opaque envelope. The envelope quickly expands around from the central massive cloud and thus greatly inflates the photospheric radius of the merger ejecta. In Sections 3 and 4 we describe two mechanisms for the envelope ejection, and calculate its expected self-similar structure. In both cases, we find a stratified structure with four-velocity $\gamma\beta$ growing outward and extending to $\gamma\beta\gg 1$. We estimate the expected mass of the ultra-relativistic envelope and its photon-to-baryon ratio.

The first mechanism of the envelope ejection is the ablation of neutron star surface at the very beginning of the merger, when it suddenly (in $30\,\mu{\rm s}$) becomes a powerful source of neutrinos (Section 3). The reaction $\nu\bar{\nu}\rightarrow e^{\pm}$ injects heat near the surface of the merging stars, resulting in huge energy per baryon and accelerating the surface layers to ultra-relativistic speeds. This effect is missed by the existing merger simulations (e.g. Dessart et al. 2009; Bauswein et al. 2013; Hotokezaka et al. 2013, 2018; Radice et al. 2016, 2018; Kiuchi et al. 2018), because they do not resolve the heating and dynamics of low-density surface layers. We find that an ultra-relativistic ablated mass m\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}10^{-8}M_{\odot} promptly escapes the vicinity of the merger, before the outflow becomes heavily polluted by baryons forming the massive cloud of slower ejecta.

The second mechanism is described in Section 4. The expanding cloud of massive ejecta can develop a strong internal shock, which propagates outward and accelerates to ultra-relativistic speeds in the outermost, low-density layers of the cloud. One appealing scenario invokes the formation of a rapidly spinning, short-lived magnetar following the merger. The fast outflow from the magnetar drives a shock into the cloud, which appears favorable for production of the “blue” kilonova (Metzger et al., 2018). We show that after the shock crosses the cloud and accelerates in its outer layers, an ultra-relativistic envelope is inflated. We estimate its mass m\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}10^{-7}M_{\odot} and describe its self-similar structure.

The presence of the envelope weakly affects the bright beamed GRBs seen by on-axis observers, however it can play a key role in off-axis GRB production. In Section 5 we discuss how a delayed launch of a jet inside the envelope leads to the production of an off-axis, single-pulse GRB, and compare our model predictions with GRB 170817A. We find that the blast wave from the jet in the envelope can explain the observed luminosity, hardness, and light curve of the burst. Comparison with previous models and observational implications are discussed in Section 6.

2. Relativistic motion in GRB 170817A

Kasliwal et al. (2017) found a lower limit for the Lorentz factor of the GRB source $\Gamma>2.5$. They used the considerations of photon-photon opacity due to collisions $\gamma\gamma\rightarrow e^{+}e^{-}$, which occur for photons of energies above $E_{\rm thr}\sim\Gamma$ MeV (e.g. Lithwick & Sari 2001). It is, however, difficult to derive robust limits from photon-photon collisions. In principle, the source is allowed to be completely opaque to MeV photons, as no such photons were observed, and their number can only be guessed from extrapolations of the observed spectrum.

Instead, a simple lower bound on $\Gamma$ can be obtained from the standard consideration of the scattering opacity in the source (e.g. Lithwick & Sari 2001). Let $L\equiv 4\pi r^{2}\Gamma^{2}n_{b}m_{p}c^{3}$ be the isotropic equivalent of the kinetic power of the relativistic outflow emitting the observed gamma-rays with luminosity $L_{\gamma}$ (note that $L$ does not include $L_{\gamma}$ and in general may be smaller or larger than $L_{\gamma}$). Here $r$ is the emission radius, $n_{b}$ is the baryon number density in the flow rest frame, and $m_{p}$ is the proton mass. A characteristic optical depth to Thomson scattering is given by

[TABLE]

where $n_{p}=Y_{e}n_{b}$ is the proton density, $Y_{e}$ is the proton-to-nucleon ratio, and \zeta\mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1 is a numerical coefficient, which depends on the radial density profile of the outflow.

If $L$ is approximately uniform on radial scales $\Delta r>r/\Gamma^{2}$ then one can show that $\zeta\sim 1$. This situation is expected if the outflow is launched on a timescale $\Delta t>r/\Gamma^{2}c$. In particular, a ballistic outflow with $\Gamma\approx const$ and $L\approx const$ has a density profile $\rho\propto r^{-2}$. The characteristic optical depth for photons (emitted isotropically in the fluid frame) can be found by integrating the scattering coefficient $\alpha_{\rm sc}\propto\rho$ along the photon trajectory. This calculation gives Equation (1) with $\zeta\sim 1$ (Abramowicz et al., 1991; Beloborodov, 2011).

In the opposite limit, one can consider an outflow ejected impulsively, within $\Delta t\approx 0$. Then its density profile is set by radial spreading during the outflow expansion to the radius of GRB emission. This radial spreading is controlled by the distribution of four-velocity, which has a significant width for any realistic ejection mechanism. In particular, the homologous envelope described later in this paper can be idealized as an impulsive ballistic ejection with a power-law distribution of four-velocity (see Section 5). A photon propagating in the homologous envelope will see a steeply decreasing density. In this case $\zeta$ can be as low as $1/4$.

The GRB radiation can escape if \tau_{\rm T}\mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1, which gives a lower limit on $\Gamma$,

[TABLE]

Another constraint is set by the minimum dispersion of the photon arrival times $\Delta t_{\rm disp}\sim r/2\Gamma^{2}c$. It applies as long as the GRB-emitting shell has a minimum angular size $\Delta\theta_{\min}\sim\Gamma^{-1}$, which is valid for any expanding relativistic cloud accelerated by its internal pressure.111The angular size of the emitting shell in GRB 170817A is likely to significantly exceed $\Gamma^{-1}$ (below we argue that \Gamma\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}5). The fact that we receive radiation at the polar angle $\theta\sim 0.5$ rad, as inferred from the gravitational wave signal, suggests \Delta\theta\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}0.5>\Gamma^{-1}. The minimum dispersion $\Delta t_{\rm disp}$ should not exceed the observed duration of the main peak of GRB 170817A $\Delta t_{\rm obs}\sim 0.5$ s, which gives the constraint

[TABLE]

Note that only the duration of the peak is relevant for this constraint; it is not affected by the arrival time of the peak $t_{\rm obs}\sim 1.7$ s.

Combining the two constraints in Equations (2) and (3), and using the observed $L_{\gamma}\approx 1.6\times 10^{47}$ erg/s and $\Delta t_{\rm obs}\approx 0.5$ s (Goldstein et al., 2017), we find

[TABLE]

Figure 1 shows the constraints on the GRB source in the $r$-$\Gamma$ plane. If most of the observed luminosity comes from the photosphere of the explosion, the source must be located near the line of $\tau_{\rm T}=1$. If it is produced by a photospheric shock breakout then it should also be located near the line of $r/2\Gamma^{2}c=0.5$ s, and so in this case it will be near the intersection of the two lines.

A similar minimum $\Gamma\sim 5$ for GRB 170817A may be obtained by considering the scattering opacity of $e^{\pm}$ pairs created by MeV photons, using a high-energy extrapolation of the observed spectrum (Matsumoto et al., 2019).

One could also consider the possibility that the source is magnetically dominated, i.e. powered by Poynting flux from the central engine and carries a negligible amount of baryonic matter. Then the plasma emitting the GRB would be made entirely of $e^{\pm}$ pairs created in photon-photon collisions. The source becomes transparent and releases radiation when its compactness parameter $\ell=L_{\gamma}\sigma_{\rm T}/m_{e}c^{3}R\Gamma^{3}$ decreases to $\sim 10$, which requires $\Gamma>10\,r_{11}^{-1/3}$.

3. Ablation of the neutron star surface

3.1. Hot sandwich at the collision interface

The two merging stars are strongly deformed by tidal forces, forming cusps pointing approximately toward each other (but still significantly misaligned, because of the orbital rotation of the binary). The merger begins with the tangential collision of the cusps, forming a growing interface between them (see e.g. Figure 4 in Bauswein et al. 2013). The surface layers at the interface are shocked, compressed, and heated, forming a thin “sandwich.” Local thermodynamic equilibrium and nuclear statistical equilibrium are quickly established in the sandwich, with pressure contributions from nuclei, electrons, positrons, and Planckian radiation, at a high temperature $T$.

The opposite tangential velocities $v_{\parallel}$ of the colliding stars imply a huge velocity shear at the interface, which immediately leads to Kelvin-Helmholtz instability (Price & Rosswog, 2006; Kiuchi et al., 2018), with a growth rate comparable to the shear rate. The limited numerical resolution of the global merger simulations makes it difficult to observe the fast shear damping that develops on smallest scales, and local simulations (Zrake & MacFadyen, 2013) show more details of the Kelvin-Helmholtz instability. The efficient damping of the tangential motion suggests that the collision is “sticky,” releasing its entire specific energy $v_{\star}^{2}/2=(v_{\parallel}^{2}+v_{\perp}^{2})/2$, not just the normal component $v_{\perp}^{2}/2$. Then the energy released in the sandwich per unit baryon rest mass is given by

[TABLE]

It may exceed 0.03 in mergers of massive neutron stars.

The hot compressed sandwich is bounded by two shocks propagating into the stars with speed $v_{\rm sh}\sim v_{\star}$. The energy density in the sandwich is $U\approx\xi\rho c^{2}\hat{E}$, where $\rho$ is the upstream (pre-shock) density, $\xi=(\hat{\gamma}+1)/(\hat{\gamma}-1)$ is the shock compression factor, and $\hat{\gamma}$ is the adiabatic index of the post-shock matter. The sandwich pressure is

[TABLE]

As the two shocks propagate into denser subsurface layers of the colliding stars, the sandwich pressure grows, $P\propto\rho$.

The approximate pressure balance between the shocks implies $P(x)\approx const$ across the sandwich, where the $x$-axis is chosen normal to the collision interface.222Variation of $P$ in the $x$-direction is small because the shocks are in hydrodynamical causal contact and not far from pressure equilibrium. At the same time, $P(y,z)$ strongly varies along the collision interface, decreasing from the center (the initial touch point of the colliding stars) to the outer parts of the sandwich, where shocks form later (Figure 2). Initially the collision interface area grows superluminally, faster than the shocked matter could be squeezed out from the sandwich. At a later stage, the pressure gradient in the $y$-$z$ plane begins to drive a fan-like “fountain” from the sandwich. However, the nucleon density is strongly non-uniform in the $x$ direction — the sandwich is made of layers of stratified density. The older layers in the middle were shocked at a low pressure and later pressurized through compression. As $P$ grows proportionally to the upstream density $\rho$, the old shocked layers are strongly compressed to stay in the approximate pressure equilibrium with the propagating shocks. This compression implies strong adiabatic heating, which produces a low-mass, thin layer of ultra-relativistic material. The compressional heating is discussed in more detail in Appendix.

In principle, the compressed layer heated to specific enthalpy $W>c^{2}$ could be partially ejected with an ultra-relativistic speed. The basic effect may be illustrated by a vessel of hot gas compressed to a small volume by external pressure. The work performed to compress the vessel is stored in the gas internal energy, and a sufficiently strong compression makes the gas relativistically hot, $W\gg c^{2}$. When the external pressure is eventually removed, the gas will explode and achieve ultra-relativistic speeds. This toy model does not, however, capture additional effects expected at the interface of the colliding stars. In particular, mixing and transport effects should suppress the relativistic ejection (see discussion in Appendix).

The magnetospheres of the colliding neutron stars will also be strongly compressed in collision. In an idealized model, this would create a magnetic “pillow” at the interface between the stars, with magnetic pressure in an approximate balance with the ram pressure of the two shocks propagating into the stars $B^{2}/8\pi\sim P$. As the ram pressure grows up to $P\sim 10^{34}\rho_{14}$ erg cm*-3* (see Equation 6), the magnetic field in the pillow is amplified up to $B\sim 5\times 10^{17}\,\rho_{14}^{1/2}$ G, where $\rho$ is the matter density upstream of the shocks. This implies compression of the magnetosphere by a huge factor $B/B_{0}$ for reasonable pre-merger magnetic fields $B_{0}$. The resulting pillow thickness $\delta\sim(B_{0}/B)R$ is many orders of magnitude smaller than the stellar radius $R\sim 10$ km.

However, the idealized picture of a compressed magnetic pillow is destroyed by instabilities of the Raleigh-Taylor type, which will tend to mix the compressed magnetic field into dense stellar material. Note also that the magnetic field is amplified in a much thicker layer as a result of Kelvin-Helmholtz instability.

3.2. Neutrino emission from the sandwich

When the two stars have just touched, the sandwich pressure and temperature are initially modest; at this earliest stage most of the collision energy converts to Planckian radiation. When the shocks propagate into deep and dense layers, \rho\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}10^{12} g cm*-3*, the pressure becomes dominated by nucleons rather than radiation. The simplest estimate for the sandwich temperature is given by the upper limit,

[TABLE]

which neglects any contributions to pressure other than nucleons.

The rate of neutrino and antineutrino production by the sandwich is quickly increasing with temperature $T$, and becomes significant before $T$ approaches $T_{\max}$. The cooling becomes significant when the shocks reach crustal layers with densities $\rho\gg 10^{10}$ g cm*-3*; then the cooling timescale becomes shorter than the shock age (the details will be described elsewhere). Neutrinos are mainly produced by the $e^{\pm}$ capture reactions $e^{+}+n\rightarrow p+\bar{\nu}$ and $e^{-}+p\rightarrow n+\nu$. If the neutrinos escape, their energy flux can be estimated as

[TABLE]

This simply states that most of the energy released in the shock converts to the neutrino flux. The effective temperature $T$ of escaping neutrinos is approximately related to their energy flux by $F\approx\sigma T^{4}$, where $\sigma=ac/4$ is the Stefan-Boltzmann constant.333The emitted neutrino flux is not exactly $\sigma T^{4}$ for two reasons: (i) their spectrum is not exactly thermal, and (ii) even for completely thermalized neutrinos, their (fermion) statistics is different from photon statistics. The numerical factor resulting from these corrections is $\sim 1$ and weakly affects the estimate in Equation (9). This gives an estimate,

[TABLE]

It may be viewed as a lower bound on $T$ at large densities. The upstream $\rho$ appearing in Equations (8) and (9) is lower near the edges of the sandwich, where the shocks formed later and therefore had less time to propagate into deep subsurface layers (Figure 2). The sandwich size measured along the collision interface grows from the initial contact point to $\sim 10$ km on a timescale $t\sim 3\times 10^{-5}$ s.

As the two shocks bounding the sandwich propagate into the subsurface layers of increasing density $\rho$, the post-shock temperature grows and so do the energies of emitted neutrinos $\epsilon_{\nu}$ (in units of $m_{e}c^{2}$). The cross section for neutrino interaction with matter grows as $\epsilon_{\nu}^{2}$, and the neutrinos eventually become absorbed near (or inside) the sandwich. This occurs when the shocks propagate sufficiently deep below the stellar surface, $\rho\gg 10^{11}$ g cm*-3*. Then neutrino transport occurs in a diffusive regime, with the local neutrino density $U_{\nu}$ close to local thermodynamic equilibrium. The mean energy of the thermalized neutrinos (in units of $m_{e}c^{2}$) is

[TABLE]

Cross sections for neutrino interactions with nucleons and leptons are summarized e.g. in Chen & Beloborodov (2007). For our estimates it is sufficient to include one process — neutrino absorption by nucleons, which has the largest opacity $\kappa_{\nu}\approx 3\times 10^{-20}\,\Theta^{2}$ cm2 g*-1*. The mean free path of neutrinos is

[TABLE]

One can see that at densities $\rho\gg 10^{11}$ g cm*-3* the neutrinos are self-absorbed and thermalized.

At the advanced stage of collision, most of the produced neutrinos are trapped in the middle of the sandwich, and the escaping neutrinos diffuse out from its outer parts (Figure 2), where density \rho\mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}10^{12} g cm*-3*. The simple estimate in Equation (8) then suggests a characteristic value for the escaping flux $F\sim 10^{41}$ erg s*-1* cm*-2*. This estimate is, however, crude. Accurate spatial distribution of $F$ could be found with expensive three-dimensional simulations of the collision dynamics and neutrino transport. The result will depend on the masses of the colliding stars. Massive mergers produce high $F$ because they have a high collision speed $v_{\star}$.

The flux $F\sim 10^{41}$ erg s*-1* cm*-2* is emitted from the “neutrino-sphere” — the surface where the neutrinos decouple from matter and begin to stream freely. The temperature of the neutrino-sphere $T_{\star}$ satisfies the approximate relation $F\sim\sigma T_{\star}^{4}$ which gives $\Theta_{\star}\sim 30$.

3.3. Heating of surface layers by reaction $\nu\bar{\nu}\rightarrow e^{+}e^{-}$

Some of the neutrinos and anti-neutrinos escaping from the neutrino-sphere can collide with each other and convert to $e^{\pm}$ pairs, depositing heat. The rate of this “neutrino heating” is independent of the local matter density and in the low-density regions near the stellar surface it injects huge energy per nucleon, rising the local specific enthalpy,

[TABLE]

to relativistic values w\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1. The heated surface layers will expand with relativistic speeds and leave the star. Thus a fraction of the stellar crust will be ablated by neutrino heating.

Neutrinos collide and turn into $e^{\pm}$ pairs with a significant cross section when there is a significant angle between their directions, \delta\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1. The neutrino heating wave propagates with a speed $v\sim c\cos\delta$, comparable to $c/2$. The wave is faster than the shocks, so ablation of surface layers occurs before the shock arrival (Figure 2).

We wish to obtain a rough estimate for the mass ablated to highly relativistic speeds. The first step is to evaluate the heating rate $\dot{Q}$ due to reaction $\nu\bar{\nu}\rightarrow e^{+}e^{-}$. The cross section for this reaction is given by (Goodman et al., 1987),

[TABLE]

where $\sigma_{0}\approx 1.7\times 10^{-44}$ cm2, $\epsilon=E/m_{e}c^{2}$ (with subscripts $\nu$ and $\bar{\nu}$ corresponding to the colliding neutrino and anti-neutrino), ${\mathbf{p}}=(E/c,p_{x},p_{y},p_{z})$ is the four-momentum of the colliding neutrino/anti-neutrino, and $\delta$ is the angle between their directions.

Each reaction $\nu\bar{\nu}\rightarrow e^{+}e^{-}$ creates an $e^{\pm}$ pair with total energy $E_{\nu}+E_{\bar{\nu}}$. The corresponding energy deposition rate per unit volume $\dot{Q}$ may be estimated by replacing the quasi-thermal spectrum of the neutrino-sphere by a delta-function at the average neutrino energy, $\epsilon_{\nu}\approx\epsilon_{\bar{\nu}}\approx 3\Theta_{\star}$. Then one finds

[TABLE]

We will also use the estimate $F_{\nu}\approx F_{\bar{\nu}}=F\approx\sigma T_{\star}^{4}$. This gives

[TABLE]

where \overline{(1-\cos\delta)^{2}}\mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1 is a numerical factor obtained after averaging over the directions of the colliding neutrinos.

It may be instructive to compare the estimate (15) with previous numerical simulations of steady heating by the reaction $\nu\bar{\nu}\rightarrow e^{+}e^{-}$ by Birkl et al. (2007) and Zalamea & Beloborodov (2011), who focused on accretion disks around accreting black holes, in Kerr metric. Birkl et al. (2007) also calculated steady heating in spherical geometry, with and without gravitational bending of neutrino trajectories. These simulations gave the efficiency of converting neutrino flux to heat, $F_{h}/F$, where

[TABLE]

$s>0$ is the altitude above the neutrino source, $F_{h}$ is the vertically integrated heating rate, and $H$ is the characteristic scale-height of the heating region. In the simplest spherical model, the thermal neutrino source is described by its surface flux $F$ and radius $R$. The characteristic $H$ is comparable to $0.1R$. We have checked that the heating rate calculated by Birkl et al. (2007), scaled to $R\sim 10$ km and $F\sim 10^{41}$ erg s*-1* cm*-2*), is approximately consistent with the estimate in Equation (15). It gives a rough estimate of the heating efficiency $F_{h}/F\sim H\dot{Q}/F$ approaching 0.1.

In a real merger the neutrino source geometry is neither spherical nor axisymmetric, and the neutrino-sphere is not parallel to the stellar surface (Figure 2). Furthermore, an essential difference from the previous work is that here we deal with an initial-value problem rather than a steady state. The ablation of surface layers is triggered by the suddenly arising burst of neutrinos from the sandwich.

It takes a very short time $t\sim x/v_{\rm sh}\sim 10^{-5}$ s for the sandwich density and temperature to reach high values so that its neutrino flux approaches $F\sim 10^{41}$ erg s*-1* cm*-2*. As the sandwich size grows, the area of the neutrino-sphere grows on a similar timescale. Thus, the local $\dot{Q}(t)$ measured at the stellar surface is a steeply increasing function of time, shaped by the evolution of the neutrino-sphere and the emitted flux $F$. This function is also slightly affected by a propagation delay: the wave of heat injection propagates with a speed of v_{h}\sim c\cos\delta\mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}c from the neutrino-sphere to the stellar surface.

3.4. Estimate for relativistically ablated mass

Only the uppermost ablated layers, which have low densities, reach highly relativistic speeds. The relativistic ablation ends after a short time $t_{a}$ when sufficient amount of matter is lifted from the NS surface and fills the main heating region $H\sim 1$ km. Later the neutrino-driven outflow becomes a relatively slow quasi-steady wind, which was studied previously in detail (Qian & Woosley, 1996; Thompson et al., 2001, 2004; Metzger et al., 2007; Dessart et al., 2009). The quasi-steady wind at $t\gg t_{a}$ is mainly heated through neutrino absorption by baryons. By contrast, heating at the initial ablation stage $t<t_{a}$ is dominated by the process $\nu\bar{\nu}\rightarrow e^{+}e^{-}$. There is much less matter in the initial relativistic outflow, and so neutrinos are mainly absorbed in $\nu\bar{\nu}$ collisions rather than by baryons.

Ablation may be roughly described as a two-step process: (1) enthalpy $w$ is deposited in the heating zone s\mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}H\sim 1 km, where matter begins its acceleration, and (2) adiabatic expansion at $s>H$ converts enthalpy to bulk kinetic energy. For instance, layers with asymptotic $\gamma\beta\sim 1$ are still relatively slow in the heating zone s\mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}H. Their modest characteristic speed $v_{a}\sim(0.1-0.2)c$ gives them time t_{a}\sim H/v_{a}\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}10^{-5} s to accumulate enthalpy $w\sim 1$ before leaving the heating zone.

One can estimate the mass of ultra-fast ablated layers from the condition w\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1 in the heating zone. The deposited energy density

[TABLE]

will give dimensionless enthalpy $w>1$ in layers of density $\rho<(U+P)/c^{2}\sim 10^{9}$ g cm*-3*. These layers occupy an initial volume V\sim Ah\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}10^{16} cm3, where $A\sim 10^{11}-10^{12}$ cm2 is the area of the neutrino-sphere near the sandwich edge, where ablation occurs, and $h\sim 0.3-1$ km is the characteristic hydrostatic scale-height of the stratified crust in the colliding tidal cusps of the stars. This gives the relativistically ablated mass $\sim V\,\dot{Q}t_{a}/c^{2}$.

The relativistic ablation is not isotropic and may peak in a solid angle $\Omega_{a}\sim 1$, which corresponds to a beaming factor of $b\sim 4\pi/\Omega_{a}$. The isotropic equivalent of the ablated mass $m$ viewed within $\Omega_{a}$ is given by

[TABLE]

We conclude that the observed isotropic equivalent of the ablated mass with $\gamma\beta>1$ may exceed $10^{26}$ g, depending on the viewing angle and the precise values of the beaming factor $b$ and the neutrino-sphere flux $F$.

3.5. Numerical simulation

The features of relativistic ablation discussed above are illustrated by the following simplified numerical model. Let us replace the merging stars by a single sphere of radius $R$ and choose the surface heating rate $\dot{Q}$ in the form

[TABLE]

Here $s=r-R_{\star}>0$ is the distance from the neutrino-sphere $R_{\star}$, which is somewhat below the stellar surface, $R_{\star}<R$. Our simulation assumes the neutrino-sphere radius $R_{\star}=10$ km and the star radius $R=11$ km.

A true spherical source of neutrinos would have $a=R_{\star}$ in Equation (19). However, our spherical model is designed as a proxy for the colliding stars, and we use $a<R_{\star}$ to parameterize the relatively small size of the ablation region near the sandwich edge (Figure 2). The dependence of $\dot{Q}$ on $s$ has two parts: the power law $(1+s/a)^{-4}$ describes the reduction of $n_{\nu}n_{\bar{\nu}}$ with distance from the neutrino-sphere, and the term in square brackets roughly describes the dependence $\sigma_{\nu\bar{\nu}}\propto\overline{(1-\cos\delta)^{2}}$, where $\tan\delta\sim a/s$ represents the characteristic angle between the colliding neutrinos. Equation (19) implies a steep decline of $\dot{Q}$ with $s$. This leads to a modest characteristic thickness of the main heating region, comparable to or smaller than 1 km.

The $\dot{Q}_{0}(t)$ in Equation (19) describes the time dependence of surface heating by the neutrino wave from the sandwich. At early times, $\dot{Q}_{0}(t)$ is a steep function, because it is proportional to $F^{9/4}$, and $F$ is quickly increasing as the sandwich pressure grows, $P\propto\rho$ (Equation 6). We do not know the exact shape of this function and replace it with a simple power law,

[TABLE]

followed by constant $\dot{Q}_{0}(t)=\dot{Q}_{0}(t_{0})$ at $t>t_{0}$. Our sample model will have $\zeta=9$ and $t_{0}=10^{-5}$ s. This crude description of the heating onset captures its main feature: the steep rise on a short timescale, which will lead to relativistic ablation of the surface layers.

The simulation tracks the dynamics of the outer crust with total mass $M_{\rm sim}=10^{28}$ g. The crust initially occupies a spherical shell with the outer radius $R=11$ km. It is initially static and stratified in hydrostatic equilibrium with a power-law index $q$: $\rho\propto x^{q}$ where $x=R-r$ is the depth below the stellar surface at a radius $r<R$. We have run simulations with $q=5$ and $q=3$. Our simulations are performed using a relativistic Lagrangian hydrodynamic code described in Lundman et al. (2018), with some modifications. In particular, we use a non-uniform discretization in the mass coordinate $0<m<M_{\rm sim}$, which allows us to resolve well the dynamics of the low-density layers near the stellar surface. The simulated mass $M_{\rm sim}$ is discretized into $10^{4}$ subshells.

The results are shown in Figures 3-5. One can see that the heat is deposited in the surface layers on the timescale $t\sim(3-7)\times 10^{-5}$ s. Then the layers leave the main heating zone, expand and cool, converting the accumulated heat to bulk kinetic energy. The relativistically ablated matter approaches its final (asymptotic) momentum $\gamma\beta\gg 1$ after a longer time, when it has lost its enthalpy through adiabatic cooling.

Figure 5 shows how much mass escapes with $\beta\gamma$ larger than a given value. The result is sensitive to the geometric parameter $a$. For instance, in the most “optimistic” case of $a=4$ km, mass $m\approx 10^{26}$ g is ejected with $\gamma\beta>2$, and $m\approx 10^{25}$ g is ejected with $\gamma\beta>4$. One can also see that the distribution $\gamma\beta(m)$ extends to very high values of $\gamma\beta\sim 10^{3}$. This is expected, as at the onset of ablation the neutrino wave deposits comparable energy everywhere in the heating zone in the upper crust (or even above it), regardless of the local density of matter. As a result, the outermost layers of the neutrino-driven outflow form an ultra-relativistic fireball that freely expands with acceleration by the radiation pressure $\gamma\propto r$. The fireball Lorentz factor is limited to $\gamma\sim 10^{3}$, because at higher $\gamma$ the fireball becomes transparent to radiation and acceleration becomes inefficient.

4. Early internal shock in the merger ejecta

This section describes the second mechanism for generating a self-similar ultra-relativistic envelope. It invokes a shock wave accelerating in the outer layers of the massive cloud around the merger remnant.

This mechanism is similar to shock breakout in supernovae, see e.g. detailed calculations in Tan et al. (2001) for spherical shocks, which are parallel to the surface of the supernova progenitor. These calculations were applied to merging neutron stars by Kyutoku et al. (2014). We note that the merger shocks propagate obliquely to the neutron star surface. When the oblique shock approaches the upper crust, the sudden density drop will make the shock perpendicular to the surface, rather than parallel, and this effect introduces an upper cutoff on the velocity of the ejected material (Matzner et al., 2013). Therefore, it is unclear if a sufficient amount of ultra-relativistic ejecta can be produced by the merger shocks launched inside the neutron stars. Instead, below we focus on internal shocks in the large cloud around the merger remnant. Shock acceleration in the expanding cloud is more capable of ejecting a relativistic envelope.

Clouds around the merger remnants were predicted to have masses $M_{\rm ej}\sim 10^{-2}M_{\odot}$ and expansion speeds $v_{\rm ej}/c\sim 0.1-0.3$ (e.g. Bauswein et al. 2013; Siegel & Metzger 2017). The cloud inferred from the kilonova in GW 170817 is somewhat more massive, $M_{\rm ej}\sim 0.05M_{\odot}$. An internal shock must develop in the cloud if matter ejected at later times has a higher speed. In particular, this is expected if the merger produces a magnetar — a massive, differentially rotating neutron star which generates ultrastrong magnetic fields. Although the field dynamics in the merger is not fully understood, it is plausible that an ultrastrong field is generated and then rises to the stellar surface due to magnetic buoyancy (Kluźniak & Ruderman, 1998). Then a strong magnetosphere should form, and the neutrino-driven wind from the massive neutron star becomes centrifugally accelerated (Mestel & Spruit, 1987; Metzger et al., 2018). This fast wind will drive a shock wave into the earlier, slower ejecta. Metzger et al. (2018) argued that such a wind could help explain the “blue” part of the kilonova emission in GW170817.

A shock with a velocity jump $\delta v$ crossing the cloud of mass $M_{\rm ej}$ will dissipate energy

[TABLE]

As long as the shock propagates deep inside the cloud, its speed remains mildly relativistic. When it reaches the outer layers of the cloud, where density is lower, the shock will accelerate and bring the outer layers to ultra-relativistic speeds. This process of shock breakout will inflate a self-similar relativistic envelope around the merger remnant.

Below we estimate the mass of the relativistic envelope and its structure, and illustrate it with a hydrodynamical simulation. Then we estimate the photon-to-baryon ratio in the inflated envelope.

4.1. Previous results for shock breakout in static clouds

Shock acceleration in non-relativistic hydrodynamics was studied in detail six decades ago (see Zel’dovich & Raizer (1967) and references therein). Its relativistic version was proposed as a possible mechanism for outflows in GRBs (Paczyński, 1998).

Tan et al. (2001) provided analytical fits for relativistic mass ejection by shock breakout, which were tested against hydrodynamical simulations. In their simulations a shock emerges from an initially static star with density $\rho=0$ at the stellar surface, and a polytropic mass stratification with depth $x$ below the surface, $\rho\propto x^{n}$, with a typical $n=3$. The form of their analytical approximation is motivated by earlier results (obtained in the relativistic and non-relativistic limits), which are valid for more general density profiles. Therefore, similar ejecta are expected for shock breakout in clouds with different density distributions, as long as $\rho$ steeply drops in the outer layers of the cloud.

The main dimensionless parameter of the problem is the ratio ${\cal E}/M_{\rm ej}c^{2}$. For strong shocks in the merger ejecta, we expect

[TABLE]

which corresponds to $\delta v\sim(0.1-0.5)c$ during the shock propagation inside the cloud, before the breakout.

The breakout problem has two parts: shock dynamics and subsequent expansion of the shock-heated fluid, with adiabatic cooling and bulk acceleration. The growth of the shock speed $\beta_{s}=v_{s}/c$ with decreasing density $\rho\ll M_{\rm ej}/r^{3}$ is approximately described by

[TABLE]

where $\gamma_{s}=(1-\beta_{s}^{2})^{-1/2}$ and $A\sim 1$ is a numerical factor. The power-law index $\alpha\approx 0.2$; its more accurate value is 0.187 when $\gamma_{s}\beta_{s}\ll 1$ (Zel’dovich & Raizer, 1967) and $\alpha=\sqrt{3}-3/2\approx 0.23$ when $\gamma_{s}\beta_{s}\gg 1$ (Johnson & McKee, 1971; Pan & Sari, 2006). As $\beta_{s}$ grows with decreasing $\rho$, the dissipated energy per unit mass increases, however the energy density decreases. Most of the energy ${\cal E}$ is dissipated in the dense, heavy part of the ejecta, and only a fraction of ${\cal E}$ is delivered to the outer, low-density layers that eventually develop highly relativistic motion $\gamma\beta>1$.

Figure 6 in Tan et al. (2001) shows how the ejecta energy is distributed over the asymptotic $\gamma\beta$ for several choices of $\tilde{{\cal E}}$. In particular, for $\tilde{{\cal E}}\sim 0.03$ they find that the ejecta with asymptotic $\gamma\beta>1$ carry the energy of ${\cal E}_{1}\approx 6\times 10^{-5}M_{\rm ej}c^{2}$, and these ejecta have mass $m_{1}\sim 10^{-4}M_{\rm ej}$. When a similar estimate is applied to the merger cloud, it gives $m_{1}\sim 10^{-6}M_{\odot}(M_{\rm ej}/0.01M_{\odot})$. Ejecta with yet higher $\gamma\beta\gg 1$ have a significantly smaller mass, e.g. $m\sim 10^{-6}M_{\rm ej}$ for $\gamma\beta>3$.

These estimates are sensitive to $\tilde{{\cal E}}$. One can see from Figure 6 in Tan et al. (2001) that a change of $\tilde{{\cal E}}$ from 0.03 by a factor of 3 (in either direction) changes $m_{1}$ by approximately two orders of magnitude. A crude estimate in the relevant parameter range may be written as

[TABLE]

[TABLE]

While Tan et al. (2001) considered a static star with a certain density profile, similar order-of-magnitude estimates apply to shocks in expanding clouds. In Section 4.2 we perform a detailed calculation for a sample cloud model and find the accurate distribution of the asymptotic four-velocity in the ejected envelope.

4.2. Simulation of early shock breakout from an expanding cloud

At the start of the simulation (time $t=t_{0}$), we specify the cloud parameters as follows. We place a spherical shell of mass $M_{\rm ej}=10^{-2}M_{\odot}$ with the outer radius $R_{c}=2\times 10^{9}$ cm and the inner radius of $10^{9}$ cm. The outer half of the shell is expanding with a uniform speed $v_{0}=0.1c$, and the inner half is expanding with $v_{0}+\delta v$, where $\delta v/c=(2\tilde{{\cal E}})^{1/2}\approx 0.45$ corresponds to $\tilde{{\cal E}}=0.1$.

We assume that the cloud was adiabatically cooled as it expanded from the merger remnant of radius $\sim 10^{6}$ cm, and we give the shell a low (insignificant) enthalpy $w=(U+P)/\rho c^{2}=10^{-3}$. The density profile of the shell is flat except near its outer edge, where density falls off exponentially on a scale $\Delta r=0.3R_{c}$. We choose the moderate $\Delta r/R_{c}$ keeping in mind that ejecta from neutron star mergers are initially hot, and there is significant pressure in the cloud until it strongly expands and cools adiabatically. Even if the cloud was initially ejected with a sharp edge, the pressure drop in the outermost layers will accelerate them, creating a positive velocity gradient in the radial direction. It leads to stretching of $\Delta r$ and makes the density decline at the edge smooth and gradual. Therefore, freely expanding warm clouds in general cannot have sharp edges. They are also generally expected to have a positive gradient of $v_{0}(r)$, so our assumption of $v_{0}(r)=const$ in the outer layers is a rather crude simplification of the expansion velocity profile.

Our initial condition with the velocity jump $\delta v$ inside the cloud of radius $\sim 10^{9}$ cm roughly corresponds to shock launching from the central object at time $\sim 1/2$ s after the merger. In the simulation, the jump immediately launches a forward shock in the middle of the cloud. There is also a reverse shock, however we are mainly interested in the dynamics of the forward shock, which propagates outward, reaches the low-density layers, and accelerates. Similar to the ablation simulation in Section 3, we use the Lagrangian mass coordinate $m$ counted from outside, and employ the Lagrangian relativistic hydrodynamics code of Lundman et al. (2018). The non-uniform discretization in $m$ allows us to track the evolution of the entire massive shell $M_{\rm ej}=0.01M_{\odot}$ while resolving the dynamics of the outer, low-$m$ layers, where the forward shock reaches highly relativistic speeds.

Figure 6 shows snapshots of fluid density $\rho(m)$ and four-velocity $\gamma\beta(m)$ in our simulation. It demonstrates the shock evolution and its effect on the structure of the cloud. As the forward shock enters the outer low-density layers and accelerates, it loses causal contact with the inner massive part of the cloud. The shock completely crosses the outer half of the cloud in $\sim 1/4$ s, and after this the shocked layers continue to expand with acceleration, converting heat to bulk kinetic energy. Then the ejecta become cold and ballistic.

The final profile of $\gamma\beta$ as a function of the Lagrangian mass coordinate $m$ is shown in Figure 7. It may be approximately described by power laws with different slopes $\psi$ in the regions \gamma\beta\mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1 and $\gamma\beta\gg 1$,

[TABLE]

In our sample simulation with $\tilde{E}=0.1$, we find $m_{1}\approx 10^{-5}M_{\odot}$. The power-law slope is $\psi\approx 0.18$ at m\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}m_{1} and $\psi\approx 1/4$ at $m\ll m_{1}$.

In this sample simulation we assumed a flat pre-shock velocity profile $\beta_{0}=0.1$. More realistic expanding clouds have an increasing profile $v_{0}(r)$, shaped during the cloud formation near the central object. Changing the shape of $v_{0}(r)$ slightly changes the results, as long as $v_{0}\sim 0.1$. We also run models with a faster pre-shock speed $v_{0}=0.1-0.3$; then the profile of $v_{0}(r)$ significantly affects the final distribution of $\gamma\beta$ after shock breakout. Furthermore, the detailed shape of this distribution is affected by the initial density profile in the outer layers, and the initial enthalpy in the cloud. However, in all runs we found the final $\gamma\beta(m)$ qualitatively similar to that shown in Figure 7: a shallow power law at \gamma\beta\mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1, and a steeper power law at $\gamma\beta\gg 1$, with $\psi$ close to 1/4.

4.3. Photon-to-baryon ratio in the envelope

The early shock breakout described above occurs at a modest distance $R_{\rm sh}$ from the merger. This distance depends on the timescale of launching the shock, $t_{\rm sh}$, which may be related to the magnetar formation in the center; we assume $t_{\rm sh}<1$ s. During this time the ejecta expands to a radius $r\sim v_{0}t_{\rm sh}<10^{10}$ cm, and the shock breaks out at a similar radius. It emits a burst of radiation when it reaches the ejecta photosphere, however this burst is too weak to be observed, as discussed below. Instead, the main radiative effect of the early internal shock is an increase of the photon number trapped in the ejecta. The large photon number will play a role later, when a new explosion from the merger remnant energizes the envelope at large radii r\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}10^{12} cm, and a detectable GRB is emitted (Section 5).

First let us estimate photon number carried by the massive cloud without an internal shock. The ejection of mass $M_{\rm ej}$ in time $t_{\rm ej}$ implies a characteristic mass outflow rate $\dot{M}\sim M_{\rm ej}/t_{\rm ej}=2\times 10^{31}(M_{\rm ej}/0.01M_{\odot})(t_{\rm ej}/1{\rm~{}s})^{-1}$. The ejection radius $R_{0}\sim 20$ km (the size of the merger remnant) and expansion speed $\beta=v/c\sim 0.1-0.3$ determine a characteristic density $\rho_{0}\sim\dot{M}/4\pi R_{0}^{2}\,v\sim 10^{8}\dot{M}_{31}/\beta_{-1}$ g cm*-3*. Note that the expansion speed requires an initial enthalpy per unit rest mass $w_{0}\sim\beta^{2}/2$ and hence a thermal energy density at the base of the outflow $U_{0}=(3/4)w_{0}\rho_{0}c^{2}\sim 10^{27}\dot{M}_{27}\beta_{-1}$ erg cm*-3*. This rough, order-of-magnitude estimate is sufficient to evaluate a characteristic temperature from $aT_{0}^{4}\sim U_{0}$: $kT_{0}\sim 1\,\dot{M}_{31}^{1/4}\beta_{-1}^{1/4}$ MeV. A more accurate estimate takes into account that at temperatures $\sim 1$ MeV $e^{\pm}$ pairs make a contribution to $U_{0}$ comparable that of photons; then $T_{0}$ is reduced by the factor of $(11/4)^{-1/4}\approx 0.8$.

The photon number density is $n_{\gamma}=U_{0}/2.7kT_{0}$, and the baryon number density is $n_{b}=\rho_{0}/m_{p}$. The photon-to-baryon ratio in the massive cloud is then given by

[TABLE]

The internal energy and entropy of the cloud are dominated by radiation, and during its adiabatic expansion the photon number remains approximately constant, because it is proportional to entropy. The conservation of the photon and baryon numbers implies that their ratio remains unchanged as the cloud expands.

The simple model of adiabatic expansion can be refined by including two effects. First, photon number is increased by a moderate factor $\sim 2$ when the temperature decreases below $m_{e}c^{2}=511$ keV and the $e^{\pm}$ pairs annihilate into photons. Second, when the cloud temperature drops to $kT_{\rm rec}\sim 150-200$ keV, free nucleons recombine into $\alpha$-particles, releasing 28 MeV per $\alpha$-particle (see Beloborodov (2003) for a discussion of $T_{\rm rec}$ in GRBs). The number of photons generated by the recombination may be estimated as

[TABLE]

where $Y_{e}<0.5$ is the proton-to-baryon ratio in the cloud. The produced photon number is small, comparable to that in Equation (27) and orders of magnitude smaller than $n_{\gamma}/n_{b}$ generated by the shock discussed below.

Deep inside the massive cloud the photon number may be calculated taking into account the synthesis of heavy neutron-rich elements and their $\beta$-decay at a later stage. However, nucleosynthesis of elements heavier than helium is inefficient in the outermost, low density layers of main interest for us, and so we do not include this effect here.

Let us now consider the production of photon number by an internal shock in the cloud. The shock boosts $n_{\gamma}/n_{b}$, because it generates entropy and radiation at radii $r\sim v_{0}t_{\rm sh}\gg R_{0}$. The jump conditions for a shock with speed $v_{s}=\beta_{s}c$ give the generated energy density $U\approx 2\gamma_{s}^{2}\beta_{s}^{2}\rho c^{2}$, where $\rho$ is the fluid density ahead of the shock. We are interested in the outer layers of the cloud, $m\ll M_{\rm ej}$, where the shock accelerates to \gamma_{s}\beta_{s}\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1. The cloud density $\rho(m)$ may be estimated as follows. Each layer of the cloud at radii of interest expands ballistically, with some positive velocity gradient, and the expansion is homologous, with $\rho(m,t)\propto t^{-3}\propto r^{-3}$. The parameter $\xi=r/\Delta r$ (which describes the sharpness of the density decline at the outer edge) freezes for homologous ballistic expansion, and so the density of the pre-shock layers may be written as

[TABLE]

Using this equation for $\rho$ and Equation (23) for $\gamma_{s}\beta_{s}$, one can express $U\approx 2\gamma_{s}^{2}\beta_{s}^{2}\rho c^{2}$ as a function of the mass coordinate of the propagating shock, $m$. This gives the following estimate for the post-shock energy density,

[TABLE]

Here $r\sim R_{\rm sh}$ is the shock breakout radius. Assuming that the energy density $U$ is quickly thermalized, we can estimate the blackbody radiation density $aT^{4}\approx U$, which gives the temperature

[TABLE]

The generated photon number per nucleon in the region where the shock accelerates to $\gamma_{s}\gg 1$ is given by

[TABLE]

In this last estimate we have omitted the weak dependence on $M_{\rm ej}$ and $\tilde{{\cal E}}$.

The assumption of quick thermalization of radiation in a heated flow can be verified as follows (e.g. Levinson 2012; Beloborodov 2013). Radiation must relax to a Planckian spectrum with the photon density $n_{\gamma}=U/2.7kT$ if the plasma efficiently emits photons. The two main processes of photon production by the thermal plasma is double Compton scattering and bremsstrahlung. In particular, double Compton scattering occurs with rate $\dot{n}_{\rm DC}\approx 0.1\,n_{\gamma}n\,\sigma_{\rm T}c\,\Theta^{2}$, where $n$ is the electron/positron number density and $\Theta=kT/m_{e}c^{2}$. The number of photons produced during the expansion timescale for the post-shock plasma is $n_{\rm DC}\approx\dot{n}_{\rm DC}\,r/c$. It becomes exponentially large, ensuring thermalization, if

[TABLE]

Here $\tau_{\rm T}\approx n\sigma_{\rm T}\Delta r$ is the scattering optical depth of the ejecta outside the current shock radius. It is related to the mass coordinate $m$ by

[TABLE]

where $\sigma_{\rm T}$ is Thomson cross section, and $Y_{e}=n_{p}/n_{b}$ is the proton-to-baryon ratio. For simplicity, we neglected the possible presence of $e^{\pm}$ pairs (which could only increase the photon production rate). One can see that the condition (33) is satisfied, and so the shock-generated radiation is thermalized in the shells of interest $m>10^{25}$ g, as long as the shock crosses the cloud before it expands to $r\sim 10^{10}$ cm.

4.4. Free neutrons

An important feature of the merger ejecta is their neutron-rich composition $Y_{e}<0.5$. The baryons at the base of the outflow are initially free nucleons, predominantly neutrons. As the matter expands and cools, the nucleons recombine into $\alpha$-particles, and the neutron excess implies some leftover free neutrons. Deep inside the massive cloud most of the free neutrons become locked into heavy, neutron-rich nuclei after $\sim 1$ s of expansion (e.g. Metzger et al. 2010). However, this process is less efficient in the outer layers of main interest here, $m\sim 10^{25}-10^{27}$ g, because their density is well below the typical density inside the massive cloud.

The free neutrons and ions are still well coupled by frequent nuclear collisions, so to a first approximation one can treat them as a single fluid. However, this approximation is not valid on small scales comparable to the shock thickness, and the drift of neutrons relative to the ions changes the shock dissipation mechanism (Beloborodov, 2017). In the absence of free neutrons, the shock is mediated by radiation and has a thickness comparable to the photon free path. In the presence of free neutrons, the shock is partially mediated by neutrons, which have much longer free paths. In addition, the neutron-ion collisions around the shock cause spallation of $\alpha$-particles (Belyanin et al., 2001; Beloborodov, 2003).

When the shock becomes highly relativistic, \gamma_{s}\beta_{s}\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1, the neutron-ion collisions in the shock become inelastic and generate pions. The pions immediately decay into ultra-relativistic leptons and generate a nonthermal inverse Compton cascade. In the presence of a magnetic field, the cascade would be capable of producing a significant photon number through synchrotron emission (Beloborodov, 2010; Vurm & Beloborodov, 2016). However, the envelope magnetization is likely low — both mechanisms of the envelope ejection (Sections 3 and 4) suggest that it expands outside the magnetic fields generated by the merger. This suggests weak synchrotron emission by the cascade from neutron collisions.

4.5. Photospheric emergence of the early shock

We also note that the early shock breakout radiates little energy when it reaches the photosphere of the cloud. It does not produce a detectable burst, because of the modest emission radius $r\sim R_{\rm sh}<10^{10}$ cm. The photosphere is located in the outermost layers with mass

[TABLE]

which we estimated from Equation (34) by setting $\tau_{\rm T}\sim 1$.444Density of the photospheric layers is low compared with the inner parts of the cloud. Therefore heavy nuclei are not synthesized in these layers, and there is no bound-free absorption of photons, so we assume Thomson opacity. Since $m_{\star}$ is so small, the internal shock must accelerate to a high Lorentz factor $\gamma_{s}(m_{\star})\sim 10$ as it reaches the photosphere. The energy of the shocked photospheric layers is given by

[TABLE]

It is 2-3 orders of magnitude smaller than the energy of GRB 170817A, and so the early shock breakout is hardly capable of emitting detectable radiation. The large number of photons produced by the shock inside the cloud remain trapped by the huge optical depth and experience strong adiabatic cooling.

5. GRB production in the envelope

The main conclusion from the preceding sections is that the merger GW 170817 likely ejected a low-mass, opaque envelope expanding with a stratified Lorentz factor $\gamma(m)\gg 1$. The cold ballistic envelope becomes capable of emitting a GRB only if it is reheated by some dissipation process. A simple way to accomplish this is to drive a new shock wave. Therefore, below we consider a scenario where the merger remnant produces a delayed explosion (e.g. Gottlieb et al. 2018). In particular, if the remnant is a super-massive neutron star with a limited lifetime (Lipunova & Lipunov 1998), the explosion may be associated with its collapse. The collapse is promoted by the generation of ultra-strong magnetic fields and loss of differential rotation, as well as by cooling due to neutrino emission, on a timescale of a few seconds. Then the nascent spinning black hole launches powerful, ultra-relativistic, magnetized jets.

Compared with the pre-collapse massive neutron star (the magnetar), the black hole is much more capable of launching the jets. During the collapse, the source of the baryonic wind polluting the magnetosphere of the magnetar disappears behind the event horizon. At the same time, the accretion disk of the merger debris continues to sustain a strong magnetic field threading the black hole. The Poynting flux from the black hole (of radius $\sim 5$ km) may exceed that from the magnetar, because it is more compact than the magnetar and is spinning faster. These conditions are favorable for formation of an ultra-relativistic jet, which is collimated by the surrounding slower ejecta.

Our proposed model for GRB production is schematically summarized in Figure 8. The jets first propagate inside the massive cloud and then expand into the large ultra-relativistic envelope. The forward shock from the jet (and its cocoon in the cloud) forms a blast wave which initially expands forward and sideways around the jet. Simulations by Duffell et al. (2018) suggest that the blast wave will be launched into the outer envelope if the jet itself is successful, i.e. if it exits the massive cloud. At later times the blast wave shape becomes nearly spherical as it travels with almost speed of light and has the radius $r\approx ct$. The jet must be strongly collimated, as required by the late afterglow observations of GRB 170817A (Lazzati et al., 2018; Granot et al., 2018; Lamb et al., 2018; Mooley et al., 2018a). Therefore, the blast wave in the envelope has an anisotropic power, however it is less beamed than the jet.

When viewed at large angles from the polar axis, the explosion emission will be dominated by the blast wave in the envelope rather than the jet itself. By contrast, when viewed on-axis, the jet kinetic energy will strongly dominate over the energy of its forward shock in the envelope, and the observer will see a much brighter beamed GRB emitted by the jet plasma with \Gamma\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}10^{2}.

Below we focus on the off-axis GRB expected from the blast wave in the envelope, and compare it with GRB 170817A. We will assume that the jet launches the blast wave into the envelope at a time comparable to one second after the merger.

5.1. Homologous density profile

Let us first evaluate the radial structure of the relativistic envelope before the blast wave from the jet. The envelope structure can be calculated in the ballistic approximation. The picture becomes particularly simple at long times/large radii: the envelope may be thought of as a sequence of shells ejected at $r\approx 0$ and $t\approx 0$ with different speeds, so that a shell with mass coordinate $m$ and velocity $v(m)$ has the radius,

[TABLE]

A shell $dm$ has the thickness $dr=t\,|dv/dm|\,dm$, and hence the envelope has the following density distribution (measured in the fixed lab frame),

[TABLE]

The proper density $\rho=\rho_{\rm lab}/\gamma$ is then given by

[TABLE]

where we have used the relation $d(\gamma\beta)=\gamma^{3}d\beta$.

In particular, for the four-velocity distribution of the form $\gamma\beta=(m/m_{1})^{-\psi}$ we find

[TABLE]

The mass stratification of the ultra-relativistic homologous envelope is such that the radial thickness of an outer shell $m$ occupies a radial thickness of $\Delta r/r\sim\gamma^{-2}(m)$. The outer (smaller) $m$ occupies a progressively smaller $\Delta r\sim m\,(dm/dr)^{-1}$ because it moves faster and has a stronger relativistic compression factor $\gamma^{-2}$.

5.2. Blast wave emergence at the photosphere

The edge of the envelope has a diverging Lorentz factor $\gamma$, and so it is out of reach for the blast wave. Effectively, the blast wave propagates in an unbound medium, resembling the shocks in external media that produce GRB afterglows. Unlike the standard afterglow model, here the medium is opaque, and the shock is mediated by radiation. Furthermore, the blast wave is accelerating, because the envelope density is steeply decreasing with radius and time, and its Lorentz factor $\gamma$ increases with radius. The growth of the blast wave Lorentz factor $\Gamma$ gradually allows it to catch up with faster (and less massive) outer shells of the envelope. Eventually, the blast wave emerges at the photosphere of the envelope and produces a pulse of observed emission.

The photosphere of the homologously expanding envelope is located where the column number density of protons is $\sim\sigma_{\rm T}^{-1}$. The photosphere mass coordinate $m_{\star}$ and radius $R_{\star}$ are related by

[TABLE]

Radiation produced by the blast wave begins to escape to a distant observer when the shock mass coordinate $m_{s}(t)$ approaches $m_{\star}(t)$.

Relativistic shocks in a photon-rich medium are capable of creating copious $e^{\pm}$ pairs (Beloborodov 2017; Lundman et al. 2018; Ito et al. 2018). This effect tends to prolong the photospheric emergence of the shock and delay its transition to complete transparency. However, the shock in the envelope described above will have a mildly relativistic jump, and pair creation will not be so efficient (we leave this for future study). In any case, the optical depth of pairs created by the shock is limited to $\sim 1$ by $e^{\pm}$ annihilation (Beloborodov, 2017). Therefore, radiation will begin to escape the blast wave when $m_{s}=m_{\star}$ even if pair creation is strong.

The observed delay of the shock appearance at the photosphere $R_{\star}$ is given by

[TABLE]

Hereafter we assume that the relativistic envelope was inflated by the mechanism described in Section 4, and use Equation (26). Combining Equations (41) and (42) we find

[TABLE]

Note that $\psi/(1+4\psi)=(1/8)\pm 0.01$ in the relevant range of $0.2<\psi<0.3$, and so

[TABLE]

where the observed time of the photospheric shock breakout $t_{\star}$ is expressed in seconds. Then we also find $R_{\star}$ from Equation (42),

[TABLE]

and the mass of the photospheric layers $m_{\star}$ (Equation 41),

[TABLE]

The shock power depends on the jump of the fluid Lorentz factor,

[TABLE]

The shock jump conditions give the relative velocity between the downstream ($\Gamma$) and the upstream ($\gamma)$: $\beta_{\rm rel}=(\chi^{2}-1)/(\chi^{2}+1)$, and the corresponding Lorentz factor is $\gamma_{\rm rel}=(1-\beta_{\rm rel}^{2})^{-1/2}=(\chi^{2}+1)/2\chi$. The energy dissipated by the shock as it crosses the photospheric layers of mass $\sim m_{\star}$ is given by

[TABLE]

The photospheric shock breakout radiates a pulse of radiation with energy $\sim{\cal E}_{\star}$ and the characteristic peak duration $\Delta t_{\rm obs}\approx\Delta t_{\star}\approx R_{\star}/2\Gamma_{\star}^{2}c$ determined by the Lorentz factor $\Gamma_{\star}=\Gamma(m_{\star})$ of the post-shock plasma. The peak width can be compared with its arrival time $t_{\rm obs}\approx t_{\star}$,

[TABLE]

Thus, the relative width of the observed peak, $\Delta t_{\star}/t_{\star}$, can be used as a proxy for the shock strength at the photosphere.

Elsewhere we describe in more detail the blast wave propagation through the homologous relativistic envelope, and evaluate emission expected at $t_{\rm obs}\gg t_{\star}$, after the main peak of the off-axis GRB. It is produced by the deeper shells $m_{s}\gg m_{\star}$, which were heated by the blast wave at smaller radii $r_{s}\ll R_{\star}$ and optical depths $\tau_{s}\gg 1$. The heated opaque shells behind $m_{\star}$ release their radiation with a delay $\Delta t_{\rm obs}(m_{s})$, when they expand to transparency. This delayed emission is partially thermalized and adiabatically cooled, and so it is much softer than the GRB peak.

5.3. Comparison with GRB 170817A

The above predictions can be compared with observations of GRB 170817A. Its arrival time was $t_{\star}\approx 1.7$ s and the main pulse had a width $\sim 3$ shorter than $t_{\star}$, which implies $\chi_{\star}\sim\sqrt{3}$ according to Equation (49). Substituting these values to Equations (48) we find

[TABLE]

One can see that the envelope with $m_{1}\sim 3\times 10^{27}$ g (which is in the expected range for the envelope model in Section 4) is consistent with the observed energy of the main peak of GRB 170817A, ${\cal E}\sim 4\times 10^{46}$ erg (Goldstein et al., 2017). The corresponding mass of the photospheric layers is $m_{\star}\sim 6\times 10^{25}$ g.

Furthermore, from Equations (44) and (45) we find that the pre-shock ejecta at the photosphere had Lorentz factor

[TABLE]

and the blast wave broke out at radius

[TABLE]

The Lorentz factor of the radiating plasma immediately behind the shock is

[TABLE]

The observed spectrum of the initial pulse peaked at $E_{\rm pk}\sim 10^{2}$ keV, which roughly corresponds to the average photon energy $\bar{E}\sim 10^{2}$ keV (the detailed shape of the spectrum of GRB 170817A is uncertain, because of poor photon statistics). This should be compared with the average energy of photons emitted by the blast wave at $R_{\star}$,

[TABLE]

Thus, the observed $E_{\rm pk}$ is consistent with the photon-to-baryon ratio $n_{\gamma}/n_{b}\sim 10^{4}$ expected in the envelope described in Section 4, see Equation (32).

GRB 170817A was also reported to have a soft tail of emission after the main peak. In an accompanying paper we study soft emission expected after the blast wave breaks out of the opaque relativistic envelope. It has a decreasing luminosity $L_{\rm obs}(t_{\rm obs})$ and a decreasing average photon energy $\bar{E}(t_{\rm obs})$. However, quantitative tests of the tail prediction are difficult for GRB 170817A, because its tail is barely detected and its properties are poorly known and still debated (cf. Burgess17).

6. Discussion

The observed timing of GRB 170817A and its luminosity implies ultra-relativistic expansion of the gamma-ray source, \Gamma\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}5 (Figure 1). This constraint shows that neutron star mergers eject ultra-relativistic outflows at large angles from the rotation axis, at least up to $\theta\sim 20-30^{\circ}$ (the viewing angle for GW 170817).

We have argued that this broad ultra-relativistic outflow has the form of a self-similar envelope expanding from the center with a stratified Lorentz factor. This picture follows from our investigation of a possible mechanism for ultra-relativistic ejecta. Both mechanisms described in Sections 3 and 4 inflate a self-similar envelope with a Lorentz factor profile increasing outward. This envelope contains significant mass and is opaque. In particular, the ultra-relativistic envelope ejected by a magnetar shock (Section 4) can have mass exceeding $10^{-7}M_{\odot}$. It is sufficient for inflating the GRB photosphere to radii $r\sim 10^{12}$ cm.

A plausible scenario for producing gamma-rays invokes a delayed explosion from the merger remnant. The explosion launches a blast wave into the inflated envelope, which eventually emerges at its photosphere and emits a gamma-ray burst (Section 5, see Figure 8). The burst radius $r\approx 10^{12}$ cm and Lorentz factor $\Gamma\approx 5$ (the red circle in Figure 1) as well as the predicted burst luminosity $L_{\gamma}\approx 10^{47}$ erg/s, are in agreement with observations (Section 5.3). Furthermore, the expected average energy of the emitted photons $\bar{E}\sim 10^{2}$ keV is consistent with observations.

6.1. Comparison with previous work

Our model for GRB 170817A shares some features with the shock-breakout models of Gottlieb et al. (2018), Bromberg et al. (2018), Nakar et al. (2018). However, there are important differences.

(i) The previous models require a small Lorentz factor $\Gamma\sim 1$ in order to explain the observed photon energy $\bar{E}\sim 10^{2}$ keV. These models adopted the plasma temperature behind the shock $kT\approx 50$ keV, regulated by $e^{\pm}$ creation, as discussed in earlier papers (e.g. Nakar & Sari 2012). Then the observed average photon energy $\bar{E}\approx 3kT\,\Gamma\approx 150\,\Gamma$ keV implies $\Gamma\sim 1$. By contrast, we find that the observed light curve requires \Gamma\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}5. Therefore, we conclude that the shock model with $kT\approx 50$ keV is in tension with observations of GRB 170817A.

(ii) The previous models assumed that shock breakout occurs in a photon-poor cloud, $n_{\gamma}/n_{b}\ll 10^{3}$. By contrast, the envelope inflated by the mechanism described in Section 4 is photon-rich, $n_{\gamma}/n_{b}>10^{4}$. Then the delayed explosion in the envelope emits a spectrum with a reduced energy per photon in the jet rest frame, which is consistent with the high-$\Gamma$ Doppler boost giving the observed $\bar{E}\sim 10^{2}$ keV. We have not calculated yet the detailed GRB spectrum expected from shock breakout in the photon-rich envelope. A similar problem was studied by Levinson (2012), and we are currently working on complete, first-principle simulations that will give the emitted spectra with various $n_{\gamma}/n_{b}\gg 10^{3}$. The results will be reported in a future paper. Note that our photon-rich model neglected the photon number generated by the shock itself via downstream bremsstrahlung emission while the previous models relied on this emission. The self-generation of photons could reach the required $n_{\gamma}/n_{b}\sim 10^{4}$ if the shock is slower and/or the outer layers of the envelope manage to synthesize heavy nuclei (Nakar, 2019).

(iii) In the previous models, the shock acceleration and the production of gamma-rays occurred at a well defined characteristic radius — the cloud “edge” where density suddenly and steeply drops by many orders of magnitude. In this respect, the models were similar to the canonical shock breakout in a stellar explosion. By contrast, the envelope described in this paper is equivalent to an infinite medium. The fact that at any given time the envelope extends to a finite radius becomes irrelevant, since its leading edge has a diverging Lorentz factor and is out of reach for a blast wave.

The expanding envelope may be idealized as a flow ejected impulsively from the center with a self-similar (power-law) distribution of Lorentz factor $\gamma(m)$. Its density profile is determined by $\gamma(m)$ and is also self-similar (Section 5.1). The acceleration of a blast wave launched in such an envelope occurs over a few decades in radius rather than at an edge of a cloud. This qualitatively changes the dynamics and radiation of the blast wave. It has to “chase” each layer of the envelope, and catches up with layers of higher $\gamma$ at progressively larger radii $r\propto\gamma^{2}$ until finally reaching the photosphere.

6.2. Future observational tests

Our results suggest a few observational implications that may be tested in the future.

(1) Our model for off-axis short GRBs predicts that the relative width of the gamma-ray pulse $\Delta t_{\rm obs}/t_{\rm obs}$ reflects the shock strength at the photosphere (Equation 49). The blast wave power is expected to decrease with the viewing polar angle $\theta$. This suggests that the luminosity of the gamma-ray counterpart should decrease with $\theta$ while $\Delta t_{\rm obs}/t_{\rm obs}$ should increase. Our model also predicts an anti-correlation between the pulse hardness $\bar{E}$ and relative duration $\Delta t_{\rm obs}/t_{\rm obs}$, as both are controlled by the blast wave strength at the photosphere, $\Gamma_{\star}/\gamma_{\star}$. These correlations may be tested by future observations. The blast wave is fastest when viewed on-axis, directly in front of the powerful collimated jet with \Gamma\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}10^{2}. Its gamma-ray emission is also shortest when viewed on-axis, $\Delta t_{\rm obs}/t_{\rm obs}\ll 1$. However, the on-axis luminosity from the blast wave should be outshined by the extremely bright, beamed emission from the jet itself, and therefore the above correlations should break at small $\theta$ where the jet comes into view. The jet is expected to emit a canonical short GRB many orders of magnitude brighter than GRB 170817A.

(2) Relativistic ablation of the neutron star surface at the onset of the merger (Section 3) suggests an immediate gamma-ray burst, overlapping with the gravitational wave signal. Relativistic ablation creates an ultra-relativistic fireball in a short time $t\sim 3\times 10^{-5}$ s after the two stars touch. Its energy ${\cal E}_{f}$ has a flat distribution over $\gamma$, up to enormous $\gamma\sim 10^{3}$ (Figure 5). Its outermost, fastest layers become transparent while still being radiation dominated, w\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1, and hence radiate away most of their energy, similar to the fireball models of Paczynski (1986) and Goodman (1986). Thus, a significant fraction of the ablation fireball energy ${\cal E}_{f}$ is radiated away. The fireball energy is quite uncertain though. It is sensitive to both the neutrino luminosity from the collision sandwich and the geometric parameter $a$ that describes the effective area of the neutrino-sphere near the stellar surface. In the most optimistic simulation with $a=4$ km this energy is ${\cal E}_{f}\sim 3\times 10^{45}$ erg; then the ablation fireball emits a quasi-thermal burst with luminosity up to $10^{51}$ erg/s and duration $\sim 10^{-5}$ s. Its observed temperature is close to the temperature at the base of the outflow, $kT\sim 3-6$ MeV. Even in the optimistic model, this initial, hard “ablation burst” is weak and difficult to detect, however it might become detectable with future, more sensitive detectors.

(3) The ejecta acceleration by the internal shock simulated in Section 4 has observational consequences for the kilonova emission. Optical/IR emission from GW 170817 can be explained as a superposition of a “blue” and a “red” component, that came from material with different opacities, with and without synthesized lanthanides, see Kasen et al. (2017). Waxman et al. (2018) pointed out that the data could also be fitted by a model with a simple (fixed, gray) opacity, if the emitting material was ejected with a power-law velocity distribution $v\propto m^{-\psi}$ with $\psi\sim 0.6$. For the cloud simulated in Section 4, both opacity and velocity may be expected to vary. The early internal shock creates a monotonic four-velocity distribution with a changing slope $\psi=-d\ln(\gamma\beta)/d\ln m$ (Figure 7). In the present paper we focused on the outer layers ($m\ll M_{\rm ej}$), which are relevant for the GRB production. The kilonova is emitted by much deeper layers (the massive part of the cloud, $m$ comparable to $M_{\rm ej}$). Qualitatively, one may expect that the variation of speed and density across the post-shock cloud will lead to the emission of blue and red kilonova components. The faster parts of the cloud will have a lower density, fail to produce lanthanides, and emit the blue component. The slower and denser parts may synthesize the lanthanide material of high opacity and emit the red component. However, it is unclear if our simple spherically symmetric simulations are capable of giving enough mass at low velocities needed for the observed red kilonova. The presence of heavier and slower outflow at large polar angles may need to be invoked to increase red emission and explain the GW 170817 observations. We leave the detailed analysis of this topic for future work.

(4) The ejected envelope is eventually decelerated by an external medium and produces afterglow emission for a broad range of viewing angles. Such deceleration afterglow is in general expected for dynamical ejecta from mergers (Nakar & Piran, 2018; Hotokezaka et al., 2018), regardless of the presence or absence of a collimated jet. Furthermore, Nakar & Piran (2018) argue that the initial slow rise of the afterglow of GW 170817 comes from ejecta moving toward us, along the line of sight, rather than the collimated jet viewed from the side. At present it is unclear if/when the envelope described in this paper can dominate the observed afterglow. Its calculated stratification $\gamma\beta(m)$ (Figure 7) may be used to develop a detailed afterglow model and check if the decelerating envelope could overshine the off-axis emission from the decelerating jet at late stages of the explosion.

We thank Ehud Nakar, the referee of this paper, for his comments. A.M.B. is supported by NSF grant AST-1816484, NASA grant NNX15AE26G, a grant from the Simons Foundation #446228, and the Humboldt Foundation. Y.L. is supported by NSF grant AST-1816484. C.L. was supported by the Swedish National Space Board under grant number Dnr. 107/16.

Appendix A Compressional heating at the collision interface

The pressure growth in the sandwich between the colliding stars implies a strong compressional heating of their surface layers, as seen from the following consideration. Let $\rho_{0}$ be the initial, pre-shock density of an old layer in the sandwich, and $P_{0}$ — its pressure when it was just crossed by the shock. Pressure in low-density layers is strongly dominated by radiation (and $e^{\pm}$ pairs), so their adiabatic index is $\hat{\gamma}\approx 4/3$. As long as neutrino cooling is negligible, the layer compression by increasing $P$ to a higher density $\rho_{0}^{\prime}$ occurs adiabatically,

[TABLE]

where we used pressure balance across the sandwich, $P_{0}^{\prime}\approx P$. Energy per baryon in the compressed layer grows proportionally to $P_{0}^{\prime}/\rho_{0}^{\prime}\propto P^{1-1/\hat{\gamma}}$.

Instead of energy per baryon it is more convenient to consider dimensionless enthalpy per unit mass, $w$. The initial $w_{0}=(P_{0}+U_{0})/\rho_{0}c^{2}$ is related to $\hat{E}$ defined in Equation (5) by

[TABLE]

As the sandwich pressure grows to $P\gg P_{0}$, the dimensionless enthalpy of the layer is amplified as

[TABLE]

where $\rho\gg\rho_{0}$ is the present density of matter just upstream of the shocks. Note that $\rho$ and $\rho_{0}$ are the pre-shock densities of different layers; $\rho/\rho_{0}$ should not be confused with the compression factor of the old layer due to the increasing pressure, $\rho_{0}^{\prime}/\rho_{0}$. The last equality in Equation (A3) may be slightly changed (by a numerical factor close to unity) when the shock pressure becomes dominated by nucleons, which leads to the postshock adiabatic index $\hat{\gamma}\approx 5/3$ instead of $4/3$. However, the scaling $w_{0}^{\prime}/w_{0}\propto\rho^{1/4}$ is weakly affected by this transition, because the scaling applies to the old, low-density layers, which remain radiation-dominated with $\hat{\gamma}=4/3$.

The moderate initial $w_{0}<0.1$ increases with time to $w_{0}^{\prime}>1$ when the shocks propagate into layers of density $\rho>w_{0}^{-4}\rho_{0}$. The layers compressed to relativistic enthalpy $w_{0}^{\prime}\gg 1$ have the potential of being ejected with highly relativistic speeds, if their internal energy has a chance to convert to bulk kinetic energy via adiabatic expansion without losing it to the neighboring heavy and non-relativistic layers.

The maximum compression factor may be reached close to the moment when the squeezed matter begins to leak out from the sandwich. At the beginning of the collision, the size of the collision interface in the $y$-$z$ plane expands with a superluminal speed $dr_{A}/dt$. Here $r_{A}=(y_{A}^{2}+z_{A}^{2})^{1/2}$ represents the curve where the surfaces of the two stars intersect, which defines the edge of the sandwich (this curve is not a circle, as the tangential motion of the colliding stars breaks the axial symmetry of the interface). The interface area grows with time because of the converging motion of the two stars, which brings it into contact with more material. The two shocks bounding the collision sandwich intersect at $r_{A}$, which initially grows with rate $dr_{A}/dt>c$. Later $dr_{A}/dt$ is reduced (because of the curvature of the stellar surfaces) and eventually becomes subliminal and subsonic. Then the hot interface material pushes the shocks aside (so that they no longer intersect at the edge) and leaks out with the sound speed through the edges of the sandwich. This mass loss will buffer the growth of pressure in the sandwich.

The time at which the sandwich matter begins to be ejected from the edges may be roughly estimated as $t\sim R/c\approx 3\times 10^{-5}$ s. During this time, the two shocks propagate a significant distance x\sim tv_{\rm sh}\mathrel{\hbox{\raise 2.15277pt\hbox{>}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}1 km into the deep stellar interior, where the upstream density can be quite large, $\rho\sim 10^{13}-10^{14}$ g cm*-3*. Therefore, we estimate the maximum possible compression in the sandwich using the ram pressure that corresponds to the upstream $\rho\sim 10^{14}$ g cm*-3*. The corresponding $\rho_{0}$ of layers reaching $w_{0}^{\prime}>1$ is found from Equation (A3),

[TABLE]

Recalling that $\rho_{0}$ is the pre-shock density of sub-surface layers with a hydrostatic scale-height $h\sim 0.3-1$ km, one can roughly estimate their mass as $m\sim r^{2}h\rho_{0}$, where r\mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}R\sim 10 km is the characteristic extension of the old layers along the interface. The initial volume occupied by these layers is V_{0}\sim r^{2}h\mathrel{\hbox{\raise 2.15277pt\hbox{<}\hbox to0.0pt{\hss\lower 2.15277pt\hbox{\sim}}}}10^{16} cm*-3*, and their nucleon mass is

[TABLE]

This mass has the potential of being ejected with a highly relativistic speed. However, a more detailed analysis suggests that this mechanism is hindered by two processes in the sandwich between the colliding stars: (1) The idealized picture where mass $m$ estimated in Equation (A5) resides in a very thin (strongly compressed) layer is likely incorrect, because the sandwich is also the site of Kelvin-Helmholtz instability. The high-$w$ material is likely to be dispersed into small bubbles or filaments and mixed into dense, massive, low-$w$ material before escaping the sandwich. Only a small fraction of the bubbles might be able to escape with a highly relativistic momentum and avoid sharing it with the non-relativistic matter. (2) Neutrino emission and transport tends to steal energy from the layers with high $w$ and reduce their pressure. This amplifies their compression and limits enthalpy per unit rest mass, $w$.