SHOCK BREAKOUT DRIVEN BY THE REMNANT OF A NEUTRON STAR BINARY MERGER: AN X-RAY PRECURSOR OF MERGERNOVA EMISSION

A merger would inevitably happen in an NS-NS binary once the gravity between the NSs could no longer be supported by angular momentum because of GW radiation release. After the merger, in some situations, a supra-massive NS rather than a black hole could be left with an extremely rapid differential rotation. This is supported indirectly by many emission features of short GRBs and afterglows, including extended gamma-ray emission, X-ray flares, and plateaus (Dai et al. 2006; Fan & Xu 2006; Rowlinson et al. 2010, 2013; Giacomazzo & Perna 2013; Zhang 2013). The presence of such a supra-massive NS is permitted by both theoretical simulations (Hotokezaka et al. 2011; Bauswein et al. 2013b) and observational constraints (e.g., the present lower limit of the maximum mass of Galactic NSs is precisely set by PSR J0348+0432 to $2.01\pm 0.04{M}_{\odot };\;$Antoniadis et al. 2013). On one hand, during the first few seconds after the birth of a remnant NS, a great amount of neutrinos can be emitted from the very hot NS. On the other hand, differential rotation of the NS can generate multipolar magnetic fields through some dynamo mechanisms (Duncan & Thompson 1992; Price & Rosswog 2006; Cheng & Yu 2014), making the NS a magnetar. Consequently, during the very early stage, the neutrino emission and ultra-strong magnetic fields together could drive a continuous baryon outflow (Dessart et al. 2009; Siegel et al. 2014; Siegel & Ciolfi 2015a), which provides an important contribution to form an isotropic merger ejecta. A few seconds later, the NS eventually enters into a uniform rotation stage and meanwhile a stable dipole structure of magnetic fields can form. From then on, the NS starts to lose its rotational energy via a Poynting flux-dominated wind. The spin-down luminosity carried by the wind can be estimated with the magnetic dipole radiation formula as ${L}_{{\rm{sd}}}={L}_{{\rm{sd,i}}}{(1+t/{t}_{{\rm{md}}})}^{-2}$, where ${L}_{{\rm{sd,i}}}={10}^{47}{{ \mathcal R }}_{{\rm{s,6}}}^{6}{{ \mathcal B }}_{14}^{2}$ ${{ \mathcal P }}_{{\rm{i}},-3}^{-4}\;\mathrm{erg}\ {{\rm{s}}}^{-1}$, ${t}_{{\rm{md}}}=2\times {10}^{5}{{ \mathcal R }}_{{\rm{s,6}}}^{-6}{{ \mathcal B }}_{14}^{-2}$ ${{ \mathcal P }}_{{\rm{i}},-3}^{2}\;{\rm{s}}$, and the zero point of time t is set at the beginning of the magnetic dipole radiation, which is somewhat later than the NS formation by several seconds. Here ${{ \mathcal R }}_{{\rm{s}}}$, ${ \mathcal B }$, and ${{ \mathcal P }}_{{\rm{i}}}$ are the radius, dipolar magnetic filed strength, and initial spin periods of the NS, respectively, and the convention ${Q}_{x}=Q/{10}^{x}$ is adopted in cgs units.

During a merger event, although the majority of matter finally falls into the central remnant NS, there is still a small fraction of matter ejected outward, e.g., the baryon wind blown from the NS mentioned above. Besides that component, a quasi-isotropic merger ejecta can also be contributed by a wind from a short-lived disk surrounding the NS, by an outflow from the colliding interface between the two progenitor NSs, and by a tidal tail due to the gravitational and hydrodynamical interactions. The latter two components are usually called dynamical components. These ejecta components differ from each other in mass, electron fraction, and entropy. It is difficult and nearly impossible to precisely describe the specific constitutions and distributions of a merger ejecta, which depend on the relative sizes of the two progenitor NSs, the equations of state of the NSs' matter, and magnetic field structures. In any case, according to the numerical simulations in the literature, on one hand, the mass of the dynamical components can range from $\sim {10}^{-4}{M}_{\odot }$ to a few times $\sim 0.01{M}_{\odot }$ (Oechslin et al. 2007; Bauswein et al. 2013a; Hotokezaka et al. 2013; Rosswog 2013). On the other hand, in the presence of a remnant supra-massive NS, the mass of a neutrino-driven wind is found to be at least higher than $3.5\times {10}^{-3}\;{M}_{\odot }$ (Perego et al. 2014), and the mass-loss rate due to ultra-high magnetic fields is about ${10}^{-3}-{10}^{-2}{M}_{\odot }\;{{\rm{s}}}^{-1}$ during the first 1–10 s (Siegel et al. 2014). Therefore, by combining all of these contributions (see Rosswog 2015 for a brief review), we can take ${M}_{{\rm{ej}}}=0.01{M}_{\odot }$ as a reference value for the total mass of an ejecta. Furthermore, we can adopt power-law density and velocity distributions of this mass as follows (Nagakura et al. 2014):

$$\begin{eqnarray}&&{\rho }_{{\rm{ej}}}(r,t)=\displaystyle \frac{(\delta -3){M}_{{\rm{ej}}}}{4\pi {r}_{{\rm{max}}}^{3}}{\left[{\left(\displaystyle \frac{{r}_{{\rm{min}}}}{{r}_{{\rm{max}}}}\right)}^{3-\delta }-1\right]}^{-1}{\left(\displaystyle \frac{r}{{r}_{\mathrm{max}}}\right)}^{-\delta },\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&{v}_{{\rm{ej}}}(r,t)={v}_{{\rm{max}}}\displaystyle \frac{r}{{r}_{{\rm{max}}}(t)},\quad \quad \mathrm{for}\quad \quad r\leqslant {r}_{\mathrm{max}}(t),\end{eqnarray} \tag{ 2 }$$

where r is the radius to the central NS, ${v}_{{\rm{max}}}$ is the maximum velocity of the head of the ejecta, which is probably on the order of $\sim 0.1c$, and the slope δ ranges from 3 to 4 according to the numerical simulation of Hotokezaka et al. (2013). We fix $\delta =3.5$ as in Nagakura et al. (2014). In fact, the variation of δ within a wide range cannot significantly affect the primary results of this paper except for an extremely high value (e.g., $\delta \gt 5$). The maximum radius of ejecta can be calculated by ${r}_{\mathrm{max}}(t)\approx {r}_{{\rm{max,i}}}+{v}_{\mathrm{max}}t$, where ${r}_{{\rm{max,i}}}\approx {v}_{\mathrm{max}}{\rm{\Delta }}t$, with ${\rm{\Delta }}t$ being the time on which the dipolar magnetic field is stabilized. Correspondingly, the minimum radius reads ${r}_{{\rm{min}}}(t)={r}_{{\rm{min,i}}}+{v}_{{\rm{min}}}t$ and its initial value can be determined by an escape radius as ${r}_{{\rm{min,i}}}={(2{{GM}}_{c}{r}_{{\rm{max,i}}}^{2}/{v}_{\mathrm{max}}^{2})}^{1/3}$, where G is the gravitational constant and ${M}_{c}$ is the mass of the remnant NS.

The huge energy released from a remnant NS (i.e., millisecond magnetar) could eventually drive an ultra-relativistic wind mixing Poynting flux and leptonic plasma, which catches up with the preceding merger ejecta very quickly. On one hand, if this wind always remains Poynting-flux-dominated, even until it collides with the ejecta, the material at the bottom of ejecta could be heated by absorbing low-frequency electromagnetic waves from the Poynting component. On the other hand, more likely, some internal dissipations (e.g., the ICMART processes; Zhang & Yan 2011) could take place in the NS wind to produce non-thermal high-energy emission. Subsequently, a termination shock could be formed at the interface between the wind and ejecta, if the wind magnetization has become sufficiently low there (Mao et al. 2010). As a result, the bottom of the ejecta can be heated by absorbing high-energy photons from the emitting wind region and/or by transmitting heat from the neighbor hot termination shock region. Additionally, even if we arbitrarily assume an extreme situation in which all of the wind energy is completely reflected from the interface, the material at the bottom of the ejecta would also be heated due to adiabatic compression by the high pressure of the wind. Therefore, in any case, the energy carried by the wind can always be mostly injected into the bottom of ejecta and heat the material there.

When the bottom of a merger ejecta is heated by an injected NS wind, a pressure balance can naturally be built between the wind and the ejecta bottom. On one hand, the pressure balance can gradually extend to larger radii through thermal diffusion, which, however, happens very slowly for an extremely optical thick ejecta. On the other hand, the high pressure of the bottom material can lead itself to expand adiabatically and to reach a high speed. This speed could result in the formation of a forward shock propagating outward into the ejecta.

By denoting the radius and speed of a shock front by ${r}_{{\rm{sh}}}$ and ${v}_{{\rm{sh}}}$, respectively, the increase rate of the mass of shocked ejecta can be calculated by

$$\begin{eqnarray}&&\displaystyle \frac{{dM}}{{dt}}=4\pi {r}_{{\rm{sh}}}^{2}[{v}_{{\rm{sh}}}-{v}_{{\rm{ej}}}({r}_{{\rm{sh}}},t)]{\rho }_{{\rm{ej}}}({r}_{{\rm{sh}}},t),\end{eqnarray} \tag{ 3 }$$

where ${v}_{{\rm{ej}}}({r}_{{\rm{sh}}},t)$ and ${\rho }_{{\rm{ej}}}({r}_{{\rm{sh}}},t)$ are the velocity and density of the upstream material just in front of the shock. Obviously, we also have ${{dr}}_{{\rm{sh}}}={v}_{{\rm{sh}}}{dt}$. As the shock propagates, its bulk kinetic energy, which has been previously gained from adiabatic acceleration, can again be partly converted into the internal energy of the newly shocked material. The rate of this shock heating effect can be written as

$$\begin{eqnarray}&&{H}_{{\rm{sh}}}=\displaystyle \frac{1}{2}{[{v}_{{\rm{sh}}}-{v}_{{\rm{ej}}}({r}_{{\rm{sh}}},t)]}^{2}\displaystyle \frac{{dM}}{{dt}}.\end{eqnarray} \tag{ 4 }$$

Then the total internal energy accumulated by the shock, ${U}_{{\rm{sh}}}$, can be derived from

$$\begin{eqnarray}&&\displaystyle \frac{{{dU}}_{{\rm{sh}}}}{{dt}}={H}_{{\rm{sh}}}-{P}_{{\rm{sh}}}\displaystyle \frac{d(\epsilon V)}{{dt}}-{L}_{{\rm{sh}}}.\end{eqnarray} \tag{ 5 }$$

Here ${P}_{{\rm{sh}}}={U}_{{\rm{sh}}}/(3\epsilon V)$ is an average pressure, $V\sim (4/3)\pi {r}_{{\rm{sh}}}^{3}$ is the volume of the whole shocked region experiencing an adiabatic expansion, and the fraction is introduced by considering that the shock-accumulated heat is mostly deposited at a small volume immediately behind the shock front. The product of this average pressure and the corresponding volume represents the cooling effect due to adiabatic expansion. ${L}_{{\rm{sh}}}$ is the luminosity of shock thermal emission, which is caused by the diffusion and the escape of the shock heat.

Approximately following a steady diffusion equation $L=(4\pi {r}^{2}/3\kappa \rho )(\partial u/\partial r)c$, where κ is opacity, and ρ and u are densities of mass and internal energy, respectively, we roughly estimate the luminosity of shock thermal emission by

$$\begin{eqnarray}&&{L}_{{\rm{sh}}}\approx \displaystyle \frac{{r}_{{\rm{max}}}^{2}{U}_{{\rm{sh}}}c}{\epsilon {r}_{{\rm{sh}}}^{3}+({r}_{{\rm{max}}}^{3}-{r}_{{\rm{sh}}}^{3})}\left[\displaystyle \frac{1-{e}^{-(\epsilon {\tau }_{{\rm{sh}}}+{\tau }_{{\rm{un}}})}}{\epsilon {\tau }_{{\rm{sh}}}+{\tau }_{{\rm{un}}}}\right],\end{eqnarray} \tag{ 6 }$$

where ${\tau }_{{\rm{sh}}}=\kappa M/4\pi {r}_{{\rm{sh}}}^{2}$ and ${\tau }_{{\rm{un}}}={\int }_{{r}_{{\rm{sh}}}}^{{r}_{\mathrm{max}}}\kappa {\rho }_{{\rm{ej}}}{dr}$ are the optical depths of the shocked and unshocked ejecta, respectively. Here, the former optical depth, which only influences the decrease of the shock emission after breakout, is calculated by considering that most of the shocked material is concentrated within a thin shell behind the shock front (Kasen et al. 2015b). The value of the parameter can be fixed by equating the shock luminosity during the breakout to the simultaneous heating rate, because after that moment freshly injected shock heat can escape from the ejecta almost freely. The opacity of the merger ejecta is predicted to be on the order of magnitude of $\sim (10\mbox{--}100)\;{\mathrm{cm}}^{2}\;{{\rm{g}}}^{-1}$, which results from the bound–bound, bound–free, and free–free transitions of lanthanides synthesized in the ejecta (Kasen et al. 2013). This value is much higher than the typical value of $\kappa =0.2\;{\mathrm{cm}}^{2}\;{{\rm{g}}}^{-1}$ for normal supernova ejecta. In this paper we take $\kappa =10\;{\mathrm{cm}}^{2}\;{{\rm{g}}}^{-1}$. Fairly speaking, some reducing effects on the opacity could exist, e.g., (1) the lanthanide synthesis in the wind components of ejecta could be blocked by neutrino irradiation from the remnant NS by enhancing electron fractions (Metzger & Fernández 2014) and (2) the lanthanides in the dynamical components could be ionized by the X-ray emission from the NS wind (Metzger & Piro 2014).

The temporal evolution of shock thermal emission, as presented in Equation (6), is obviously dependent on the dynamical evolution of the shock. Due to the slow thermal transmission in the optical thick ejecta, a significant pressure/temperature gradient must exist in the ejecta during shock propagation, which could lead the shock to be continuously accelerated. Therefore, the dynamical evolution of such a radiation-mediated shock is completely different from the internal and external shocks in GRB situations. For an optical thin GRB ejecta, a pressure balance can be built throughout the whole shocked region simultaneously with the shock propagation. In that case, the shock velocity can simply be derived from shock jump conditions (Dai 2004; Yu & Dai 2007; Mao et al. 2010). On the contrary, in the present optical thick case, a detailed dynamical calculation of the shock in principle requires an elaborate description of the energy and mass distributions of the ejecta (see Kasen et al. 2015b for a one-dimensional hydrodynamical simulation for a similar process), which is beyond the scope of an analytical model. Nevertheless, a simplified and effective dynamical equation can still be obtained according to the energy conservation of the system.

The total energy of the shocked region can be written as¹ $E=\frac{1}{2}{{Mv}}_{{\rm{sh}}}^{2}+U$, where $U$ is the total internal energy of the shocked region. For a radiation-mediated shock, the value of U should be much higher than ${U}_{{\rm{sh}}}$. In principle, the concept of a "shocked region" that is used here can generally include the NS wind regions, because the total mass of wind leptons is drastically smaller than that of shocked ejecta, and more importantly, the energy released from the NS is continuously distributed in both the wind and ejecta through thermal diffusion, which makes them behave like a whole. By ignoring the relatively weak energy supply from radioactivity, the variation of the total energy can be written as

$$\begin{eqnarray}&&{dE}=(\xi {L}_{{\rm{sd}}}-{L}_{{\rm{e}}}){dt}+\displaystyle \frac{1}{2}{v}_{{\rm{ej}}}^{2}({r}_{{\rm{sh}}},t){dM},\end{eqnarray} \tag{ 7 }$$

where ξ represents the energy injection efficiency from the NS wind, which could be much smaller than one if the NS wind is highly anisotropic, and ${L}_{{\rm{e}}}$ is the total luminosity of the thermal emission of merger ejecta. As a general expression, the specific form of the energy injection is not taken into account. Substituting the expression of E into Equation (7), we can find the dynamical equation of the forward shock to be

$$\begin{eqnarray}&&\displaystyle \frac{{{dv}}_{{\rm{sh}}}}{{dt}}=\displaystyle \frac{1}{{{Mv}}_{{\rm{sh}}}}\left[(\xi {L}_{{\rm{sd}}}-{L}_{{\rm{e}}})-\displaystyle \frac{1}{2}({v}_{{\rm{sh}}}^{2}-{v}_{{\rm{ej}}}^{2})\displaystyle \frac{{dM}}{{dt}}-\displaystyle \frac{{dU}}{{dt}}\right].\end{eqnarray} \tag{ 8 }$$

In order to clarify the expressions of ${L}_{{\rm{e}}}$ and ${dU}/{dt}$, note that $\tilde{U}=U-{U}_{{\rm{sh}}}$, which represents the internal energy excluding the shock-accumulated part. The evolution of this internal energy component can be given by

$$\begin{eqnarray}&&\displaystyle \frac{d\tilde{U}}{{dt}}=\xi {L}_{{\rm{sd}}}-\tilde{P}\displaystyle \frac{{dV}}{{dt}}-{L}_{{\rm{mn}}},\end{eqnarray} \tag{ 9 }$$

where adiabatic cooling is also calculated with an average pressure as $\tilde{P}=\tilde{U}/3V$ and

$$\begin{eqnarray}&&{L}_{{\rm{mn}}}\approx \displaystyle \frac{\tilde{U}c}{{r}_{{\rm{max}}}}\left[\displaystyle \frac{1-{e}^{-({\tau }_{{\rm{sh}}}+{\tau }_{{\rm{un}}})}}{{\tau }_{{\rm{sh}}}+{\tau }_{{\rm{un}}}}\right].\end{eqnarray} \tag{ 10 }$$

The above expression is different from Equation (6) because the majority of the internal energy of the shocked region is deposited in the innermost part of the region, which is much deeper than the shock front. Then we have ${L}_{{\rm{e}}}={L}_{{\rm{sh}}}+{L}_{{\rm{mn}}}$. The emission component represented by ${L}_{{\rm{mn}}}$ actually is the mergernova emission discussed in Yu et al. (2013).

By substituting Equations (4), (5), and (9) into (8), we can obtain another form of the dynamical equation as

$$\begin{eqnarray}&&\displaystyle \frac{{{dv}}_{{\rm{sh}}}}{{dt}}=\displaystyle \frac{1}{{{Mv}}_{{\rm{sh}}}}\left(\displaystyle \frac{U}{3V}\right)\displaystyle \frac{{dV}}{{dt}},\end{eqnarray} \tag{ 11 }$$

which is just the expression adopted in Kasen et al. (2015b) to calculate the shock breakout for super-luminous supernovae. This equation can be easily understood in the framework of adiabatic acceleration of a "fireball." As discussed above, what is accelerating the ejecta material is actually the local internal pressures at different radii, which are much lower than the pressure of the NS wind, due to the significant delay of pressure transmission. The work performed by these varying pressures can be effectively estimated with an average internal pressure of $(U/3V)$ with respect to a volume variation of ${dV}=4\pi {r}_{{\rm{sh}}}^{2}{v}_{{\rm{sh}}}{dt}$. With the above dynamical equation, the energy conservation and assignments of the system can be described reliably. By considering that the internal energy at the time t is on the order of magnitude $U\sim {L}_{{\rm{sd}}}t$, Equation (11) can naturally determine a kinetic energy that is also on the same order of magnitude $\frac{1}{2}{{Mv}}_{{\rm{sh}}}^{2}\sim {L}_{{\rm{sd}}}t$.