MODELS OF KILONOVA/MACRONOVA EMISSION FROM BLACK HOLE–NEUTRON STAR MERGERS

The detection of gravitational waves from a binary-black hole (BH) merger by Advanced LIGO in event GW150914 (Abbott et al. 2016) marked the beginning of the era of gravitational-wave astronomy. As well as binary-BH mergers, BH–neutron star (NS) mergers are among the most promising gravitational-wave sources for ground-based detectors, such as Advanced LIGO, Advanced VIRGO (Acernese et al. 2015), and KAGRA (Somiya 2012). The detection of gravitational waves from BH–NS mergers is expected to be achieved in the next few years as statistical studies predict that more than once detection per year could be achieved by the above detectors (Abadie et al. 2010; Dominik et al. 2015; Kim et al. 2015).

Since a BH–NS binary contains a NS, there are many possible electromagnetic counterparts for the merger event (e.g., Li & Paczyński 1998; Roberts et al. 2011; Metzger & Berger 2012; Piran et al. 2013; Rosswog et al. 2013; Tanaka & Hotokezaka 2013; Takami et al. 2014; Tanaka et al. 2014; Kisaka et al. 2015). The so-called kilonova/macronova is one of the candidate electromagnetic counterparts for BH–NS mergers as well as a NS–NS mergers (Li & Paczyński 1998). The kilonova/macronova is induced by the NS material ejected during the merger (Rosswog 2005; Kyutoku et al. 2013, 2015; Kawaguchi et al. 2015). Since the NS consists of highly neutron-rich matter, r-process nucleosynthesis is expected to take place in the ejecta (Lattimer & Schramm 1974; Metzger et al. 2010), and emission powered by decay of the radioactive nuclei would occur. Simultaneous detection of a kilonova/macronova and gravitational waves would provide rich information about a merger event. It could be useful for determining the host galaxy of the source. As its lightcurve reflects the binary parameters, it could also be useful for extracting physical information about the binary.

Recently, a variety of numerical-relativity simulations have been performed for BH–NS mergers, and the quantitative dependence of the ejecta on a wide range of binary parameters has been revealed (e.g., Foucart et al. 2014; Kawaguchi et al. 2015; Kyutoku et al. 2015). Several groups have also studied the lightcurve of the kilonova/macronova by performing radiation-transfer simulations (Roberts et al. 2011; Kasen et al. 2013; Tanaka & Hotokezaka 2013; Tanaka et al. 2014). Specifically, Tanaka et al. (2014) performed simulations focusing on BH–NS binaries based on the hydrodynamical evolution obtained from numerical-relativity simulations. These radiation-transfer simulations also consider a detailed heating rate and opacity. However, the study only showed the results for a few ejecta models, and the dependence of kilonova/macronova emission on a wide range of binary parameters has not yet been studied.

In this paper, we introduce a semi-analytic model for BH–NS kilonova/macronova emission using fitting formulas for the mass and velocity of dynamical ejecta from BH–NS mergers. We calibrate these fitting formulas by comparing them with the results of recent numerical-relativity simulations for BH–NS mergers performed by the Kyoto group. The semi-analytic model for BH–NS kilonova/macronova emission is checked and calibrated by comparing it with the results of a multi-frequency radiation-transfer simulation performed by Tanaka et al. (2014).

This paper is organized as follows. In Section 2, we derive fitting formulas for the mass and the velocity of dynamical ejecta and a semi-analytic model for a kilonova/macronova lightcurve. In Section 3, we explore the possible range of kilonova/macronova magnitudes and its dependence on the binary parameters by using the model derived in Section 2. In Section 4, we apply our model to GRB 130603B as an illustration. Finally, a summary and remarks on this work are presented in Section 5. Our notation for physically important quantities is as follows: the BH mass, M_BH; the NS mass, M_NS; the ejecta mass, M_ej; the mass-weighted root mean square of the ejecta velocity, v_ave; the mass ratio, Q = M_BH/M_NS; the dimensionless spin parameter of the BH, $\chi ={{cS}}_{{\rm{BH}}}/{{GM}}_{{\rm{BH}}}^{2}$; the angle between the BH spin and the orbital angular momentum, i_tilt; and the compactness of the NS, ${ \mathcal C }={{GM}}_{{\rm{NS}}}/{c}^{2}{R}_{{\rm{NS}}}$. G and c denote the gravitational constant and the speed of light, respectively. The BH and the NS masses are the Arnowitt-Deser-Misner (ADM) masses at infinite separation.

In this section we derive fitting formulas for the mass and velocity of dynamical ejecta from BH–NS mergers, and an analytic model for the lightcurve of a kilonova/macronova. We only consider the dynamical ejecta which are ejected during the merger process. The ejecta from the BH-accretion disk system, which could be formed in the post-merger phase, are not taken into account. We remark on the effect of this additional ejecta component in Section 5.

2.1. Ejecta Mass

Foucart (2012) introduced a fitting formula for the mass remaining outside the remnant BH after BH–NS mergers. This remaining mass includes both the remnant disk mass and ejecta mass. As the ejecta mass is defined as the gravitationally unbound component of the remaining mass, we expect that a similar form of fitting model can also be useful for the ejecta mass. Therefore, by referring to Foucart (2012), we propose a fitting model for the ejecta mass as

$$\begin{eqnarray}\displaystyle \frac{{M}_{{\rm{ej}}}}{{M}_{{\rm{NS,\ast }}}} & = & {\rm{Max}}\{{a}_{1}{Q}^{{n}_{1}}(1-2\;{ \mathcal C }){{ \mathcal C }}^{-1}\\ & & -{a}_{2}\;{Q}^{{n}_{2}}\;{\tilde{r}}_{{\rm{ISCO}}}({\chi }_{{\rm{eff}}})\\ & & +\left.{a}_{3}\left(1-\displaystyle \frac{{M}_{{\rm{NS}}}}{{M}_{{\rm{NS,\ast }}}}\right)+{a}_{4},0\right\},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\chi }_{{\rm{eff}}}=\chi \;{\rm{cos}}\;{i}_{{\rm{tilt}}},\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}{\tilde{r}}_{{\rm{ISCO}}}(\chi ) & = & 3+{Z}_{2}\\ & & -{\rm{sign}}(\chi )\sqrt{(3-{Z}_{1})(3+{Z}_{1}+2{Z}_{2})},\\ {Z}_{1} & = & 1+{(1-{\chi }^{2})}^{1/3}\\ & & \times \;\{{(1+\chi )}^{1/3}+{(1-\chi )}^{1/3}\},\\ {Z}_{2} & = & \sqrt{3{\chi }^{2}+{Z}_{1}^{2}}.\end{eqnarray} \tag{ 3 }$$

We generalize the model of Foucart (2012) as follows: (i) we set the exponents of Q, i.e., n₁ and n₂, as fitting parameters and (ii) we add a term proportional to the specific binding energy of the NS, $1-{M}_{{\rm{NS}}}/{M}_{{\rm{NS,\ast }}}$, where ${M}_{{\rm{NS,\ast }}}$ is the total baryon mass of the NS.

We use the results of recent numerical-relativity simulations performed by the Kyoto group (Kawaguchi et al. 2015; Kyutoku et al. 2015) to determine the fitting parameters, which are summarized in Table 2. We use the data for the BH–NS mergers with various mass ratios, Q, BH spin magnitudes, χ, BH spin orientations, i_tilt, and NS equations of state (EOSs). In the simulations, four phenomenological EOSs described in Read et al. (2009) are employed. We list the key quantities of a NS with the EOS models which are employed in the numerical-relativity simulations in Table 1. Table 2 lists the new results of numerical-relativity simulations. These new simulations and the computations for initial data are performed using the same method as in Kawaguchi et al. (2015). We note that the NS mass is fixed to be 1.35 M_⊙ for all the models. Lovelace et al. (2013) show that very massive ejecta (M_ej ≳ 0.2 M_⊙) can be produced in the case where the BH spin is extremely high (χ_eff ∼ 0.97). However, because Lovelace et al. (2013) is, so far, the only study for such a rapidly spinning case, and we do not use its data for the calibration, we only apply the fitting formulas for the case where χ_eff ≤ 0.9.

Table 1.  The Key Quantities of a NS with the Piecewise Polytropic EOSs (Read et al. 2009) Employed in the Numerical-relativity Simulations

Model	R1.35(km)	M∗,1.35(M⊙)	${{ \mathcal C }}_{1.35}$
APR4	11.1	1.50	0.180
ALF2	12.4	1.49	0.161
H4	13.6	1.47	0.147
MS1	14.4	1.46	0.138

Note. R_1.35, M_∗,1.35, and ${{ \mathcal C }}_{1.35}$ are the radius, the Baryon rest mass, and the compactness parameter for the isolated NS with M_NS = 1.35 M_⊙, respectively.

Download table as:  ASCIITypeset image

Table 2.  List of Results Obtained from Recent Numerical-relativity Simulations

ID	Q	χ	itilt [°]	EOS	Mej	vave [c]	Ref
1	3	0	0	APR4	<1 × 10−3	0.20	(1)
2	3	0	0	ALF2	0.003	0.22	(1)
3	3	0	0	H4	0.006	0.22	(1)
4	3	0	0	MS1	0.02	0.23	(1)
5	3	0.5	0	APR4	0.002	0.21	(1)
6	3	0.5	0	ALF2	0.02	0.24	(1)
7	3	0.5	0	H4	0.03	0.23	(1)
8	3	0.5	0	MS1	0.05	0.24	(1)
9	3	0.75	0	APR4	0.01	0.23	(1)
10	3	0.75	0	ALF2	0.05	0.25	(1)
11	3	0.75	0	H4	0.05	0.24	(1)
12	3	0.75	0	MS1	0.07	0.25	(1)
13	3	0.75	31	H4	0.03	0.22	New
14	3	0.75	62	H4	0.02	0.24	New
15	3	0.75	93	APR4	<1 × 10−3	0.21	New
16	3	0.75	93	H4	0.006	0.22	New
17	5	0.5	0	APR4	<1 × 10−3	0.23	(1)
18	5	0.5	0	ALF2	0.01	0.27	(1)
19	5	0.5	0	H4	0.02	0.26	(1)
20	5	0.5	0	MS1	0.05	0.27	(1)
21	5	0.75	0	APR4	0.008	0.25	(1)
22	5	0.75	0	ALF2	0.05	0.28	(1)
23	5	0.75	0	H4	0.05	0.27	(1)
24	5	0.75	0	MS1	0.08	0.28	(1)
25	5	0.75	33	APR4	0.005	0.30	(2)
26	5	0.75	33	ALF2	0.03	0.27	(2)
27	5	0.75	33	H4	0.04	0.27	(2)
28	5	0.75	32	MS1	0.07	0.28	(2)
29	5	0.75	63	APR4	0.001	0.27	(2)
30	5	0.75	63	ALF2	0.007	0.28	(2)
31	5	0.75	63	H4	0.01	0.25	(2)
32	5	0.75	63	MS1	0.01	0.27	(2)
33	5	0.75	94	APR4	<1 × 10−3	0.24	(2)
34	5	0.75	94	ALF2	<1 × 10−3	0.26	(2)
35	5	0.75	94	H4	0.001	0.28	(2)
36	5	0.75	93	MS1	0.01	0.27	(2)
37	7	0.5	0	APR4	<1 × 10−3	0.23	(1)
38	7	0.5	0	ALF2	<1 × 10−3	0.27	(1)
39	7	0.5	0	H4	0.003	0.29	(1)
40	7	0.5	0	MS1	0.02	0.30	(1)
41	7	0.75	0	APR4	<1 × 10−3	0.27	(1)
42	7	0.75	0	ALF2	0.02	0.29	(1)
43	7	0.75	0	H4	0.04	0.29	(1)
44	7	0.75	0	MS1	0.07	0.30	(1)
45	7	0.75	33	H4	0.03	0.28	New

Note. "New" in the "Ref" Column denotes the data obtained by our new numerical simulations. Q, χ, i_tilt, M_ej, and v_ave are the mass ratio of the binary, the dimensionless spin parameter of the BH, the initial misalignment angle between the BH spin and the orbital angular momentum, the ejected mass, and the mass-weighted root mean square of the ejecta velocity distribution, respectively.

Download table as:  ASCIITypeset image

We determine the fitting parameters from the simulation data using the least squares method. The best-fit values for the parameters were obtained as follows:

$$\begin{eqnarray}{a}_{1} & = & 4.464\times {10}^{-2},\\ {a}_{2} & = & 2.269\times {10}^{-3},\\ {a}_{3} & = & 2.431,\\ {a}_{4} & = & -0.4159,\\ {n}_{1} & = & 0.2497,\\ {n}_{2} & = & 1.352.\end{eqnarray} \tag{ 4 }$$

Figure 1 plots the comparison of the ejecta mass fitting formula with the results of numerical simulations and the relative fitting error of the data as a function of the ejecta mass. The data for M_ej ≥ 0.05 M_⊙ are fitted within ∼20%, while the data with 0.04 M_⊙ and 0.02 M_⊙ are fitted within ∼30% and ∼40%, respectively. Because the results obtained by the simulations include errors due to numerical discretization, some dispersion is unavoidable even if the fitting model is appropriate. Since only limited data in Table 2 were published with explicit error measurements, we assume the estimated numerical error of the simulation data using

$$\begin{eqnarray}&&{\rm{\Delta }}{M}_{{\rm{ej,NR}}}=\sqrt{{(0.1{M}_{{\rm{ej}}})}^{2}+{(0.005{M}_{\odot })}^{2}}\end{eqnarray} \tag{ 5 }$$

referring to the estimated numerical error discussed in the appendix of Kyutoku et al. (2015) and Kawaguchi et al. (2015). Figure 1 shows that the errors of the fitting are consistent with these estimated errors. The error of the ejecta mass induces the relative error in the peak luminosity only by about a half of its relative error because the luminosity of the kilonova/macronova is approximately proportional to ${M}_{{\rm{ej}}}^{1/2}$ in the early phase. This error in the peak luminosity is comparable to or even smaller than the systematic error of model of kilonova/macronova lightcurves, as we will see below. The reduced χ² for this fit is defined as,

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{{N}_{{\rm{NR}}}-{N}_{{\rm{p}}}-1}\displaystyle \sum _{i=1}^{{N}_{{\rm{NR}}}}\;{\left(\displaystyle \frac{{M}_{{\rm{ej}}}-{M}_{{\rm{ej,fit}}}}{{\rm{\Delta }}{M}_{{\rm{ej,NR}}}}\right)}^{2},\end{eqnarray} \tag{ 6 }$$

where N_NR = 45 and N_p = 6 are the number of data points and parameters. χ² for the model of Equation (1) is 0.85 and χ² becomes larger than 1 when we reduced the number of parameters, such as a₃, a₄, n₁, and n₂. Thus, we refrain from increasing the number of parameters in this paper to improve the fitting accuracy, in order to avoid the simulation data being over-fitted. We note that this fitting formula could have systematic errors due to the choice of NS EOSs that are used for the fitting. For example, we should check whether our fitting formula can appropriately predict the ejecta mass for two EOSs that give the same NS compactness but give different ${M}_{{\rm{NS,\ast }}}$ with the same M_NS. This can only be checked by testing the fitting model with the data using an EOS which is not used in this paper. We keep this as a future task.

Zoom In Zoom Out Reset image size

We also derive a fitting model for the averaged velocity of the ejecta as a simple linear model of Q:

$$\begin{eqnarray}&&{v}_{{\rm{ave}}}=(0.01533\;Q+0.1907)\;c.\end{eqnarray} \tag{ 7 }$$

The relative error of the fitting is always within 10%. However, as is discussed in Kyutoku et al. (2015), we expect that the ejecta velocity measured in the numerical-relativity simulation can be overestimated by ∼20% and, thus, the relative error in the velocity fitting formula can be ∼30%. This error in ejecta velocity can cause an ∼15% relative error in t_c and the bolometric luminosity in the lightcurve model introduced below, which only weakly affects the following discussion.

2.2. The Kilonova/Macronova

Here, we derive a model for the kilonova/macronova from the anisotropic ejecta with a velocity distribution resulting from BH–NS mergers with reference to the kilonova/macronova model introduced in Piran et al. (2013) and Kisaka et al. (2015). Rosswog et al. (2013), Kyutoku et al. (2013, 2015), and Kawaguchi et al. (2015) showed that the ejecta expand homologously and exhibit a crescent-like shape in most cases for BH–NS mergers (see Figure 2). We describe this morphology of the ejecta by modeling the ejecta shape as a partial sphere in the longitudinal and latitudinal directions as shown in Figure 3. We employ spherical coordinates, setting the ejecta on the equatorial plane. The latitudinal coordinate θ is measured from the equatorial plane. Due to the homologous expansion, each shell with the same radius has a velocity v = r/t, where t is the elapsed time since the merger.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Kyutoku et al. (2013, 2015) showed that, the opening angle of the ejecta arc is typically φ_ej ≈ π, while the half thickness of the ejecta in the latitudinal direction is typically θ_ej ≈ 1/5. The ejecta have an approximately flat distribution in their expanding velocity. Assuming a homogeneous mass distribution in the directions of θ and φ, and a homologous expansion of the mass shell for each velocity, the density of the ejecta is given by

$$\begin{eqnarray}&&\rho (v,t)=\displaystyle \frac{{M}_{{\rm{ej}}}}{2{\varphi }_{{\rm{ej}}}{\theta }_{{\rm{ej}}}({v}_{{\rm{max}}}-{v}_{{\rm{min}}})}{v}^{-2}{t}^{-3},\end{eqnarray} \tag{ 8 }$$

where v_max and v_min are the highest and lowest values of the radial velocity of the ejecta, respectively. Here, we assume that sin θ_ej ≈ θ_ej ≪ 1.

We use a diffusion approximation for the radiation transfer. We also assume a gray-opacity with κ = 10 cm² g⁻¹, which is shown in Kasen et al. (2013) and Tanaka & Hotokezaka (2013) to be a good approximation for determining the bolometric luminosity of lathanoid-rich kilonovae/macronovae. Because the density of the ejecta decreases in time due to the free expansion, the optical depth of the ejecta decreases with time and hence emission from inner ejecta eventually becomes visible to the observer. Here, in order to calculate the luminosity, we assume that the photons do not diffuse from the radial edge or the longitudinal edge of the ejecta, but only from the latitudinal edge. This assumption is justified by the fact that θ_ej is small for the ejecta from BH–NS mergers and the areas of the radial edge and the longitudinal edge are smaller by ∼θ_ej and ∼θ_ej/φ_ej than the area of the latitudinal edge, respectively. Thus, the contribution of diffused photons to the luminosity is dominated by the photons diffused in the latitudinal direction until the photons start to escape from the whole ejecta region. This assumption is also consistent with the results of previous radiation-transfer simulations, which show that the emission from the radial edge is smaller than that from the latitudinal edge by a factor of about 3 until t ∼ 1 day (∼t_c below; see Figure 3 in Roberts et al. 2011 or Figure 8 in Tanaka et al. 2014). At later times, the photons diffuse isotropically since the whole region of the ejecta becomes visible (t > t_c). Considering the random walk of photons, the depth of the visible mass is determined by the condition that the distance to the latitudinal edge is comparable to the distance that a photon diffuses, namely ${vt}({\theta }_{{\rm{ej}}}-\theta )\approx {ct}/\tau$. Here, $\tau \approx \kappa \rho {vt}({\theta }_{{\rm{ej}}}-\theta )$ is an optical depth measured from the latitudinal edge. From this condition, we obtain the depth of the visible mass θ_obs as

$$\begin{eqnarray}{\theta }_{{\rm{obs}}}(t) & = & {\theta }_{{\rm{ej}}}-{\left[\displaystyle \frac{2{\varphi }_{\mathrm{ej}}{\theta }_{{\rm{ej}}}({v}_{{\rm{max}}}-{v}_{{\rm{min}}})c}{\kappa {M}_{{\rm{ej}}}}\right]}^{1/2}t\\ & =: & {\theta }_{{\rm{ej}}}\left(1-\displaystyle \frac{t}{{t}_{{\rm{c}}}}\right),\end{eqnarray} \tag{ 9 }$$

where we set

$$\begin{eqnarray}&&{t}_{{\rm{c}}}={\left[\displaystyle \frac{{\theta }_{{\rm{ej}}}\kappa {M}_{{\rm{ej}}}}{2{\varphi }_{\mathrm{ej}}({v}_{{\rm{max}}}-{v}_{{\rm{min}}})c}\right]}^{1/2}.\end{eqnarray} \tag{ 10 }$$

Using this result, the mass of the photon-escaping region is given by,

$$\begin{eqnarray}{M}_{{\rm{obs}}}(t) & = & {M}_{{\rm{ej}}}\\ & & -{\displaystyle \int }_{0}^{{\varphi }_{{\rm{ej}}}}d\varphi {\displaystyle \int }_{-{\theta }_{{\rm{obs}}}(t)}^{{\theta }_{{\rm{obs}}}(t)}{\rm{sin}}\theta d\theta {\displaystyle \int }_{{v}_{{\rm{min}}}}^{{v}_{{\rm{max}}}}{{dvv}}^{2}{t}^{3}\rho \\ & = & {M}_{{\rm{ej}}}\displaystyle \frac{t}{{t}_{{\rm{c}}}}\;\;(t\lt {t}_{{\rm{c}}}).\end{eqnarray} \tag{ 11 }$$

At t = t_c, the whole region of the ejecta becomes visible. Thus for t > t_c, we set M_obs(t) = M_ej.

Following Piran et al. (2013), we assume that the observed luminosity is dominated by the energy release via radioactive decay. Korobkin et al. (2012) and Wanajo et al. (2014) showed that the specific heating rate is given approximately by a power law $\dot{\epsilon }(t)\approx {\dot{\epsilon }}_{0}{(t/{\rm{day}})}^{-\alpha }$, and we set ${\dot{\epsilon }}_{0}\;=1.58\times {10}^{10}\;{\rm{erg}}\;{{\rm{g}}}^{-1}\;{{\rm{s}}}^{-1}$ and α = 1.2 following Tanaka et al. (2014). The resulting bolometric luminosity is given by

$$\begin{eqnarray}L(t) & = & (1+{\theta }_{{\rm{ej}}}){\epsilon }_{{\rm{th}}}{\dot{\epsilon }}_{0}{M}_{{\rm{ej}}}\\ & & \times \;\left\{\begin{array}{ll}\displaystyle \frac{t}{{t}_{{\rm{c}}}}{\left(\displaystyle \frac{t}{{\rm{day}}}\right)}^{-\alpha } & t\leqslant {t}_{{\rm{c}}},\\ {\left(\displaystyle \frac{t}{{\rm{day}}}\right)}^{-\alpha } & t\gt {t}_{{\rm{c}}}.\end{array}\right.\end{eqnarray} \tag{ 12 }$$

Here, _th is the efficiency of thermalization introduced in Metzger et al. (2010). The factor $(1+{\theta }_{{\rm{ej}}})$ is introduced to include the contribution from the radial edge effectively.

In order to check the validity of the analytic model obtained above, we compare Equation (12) with the results obtained from radiation-transfer simulations in Tanaka et al. (2014)⁵ using the ejecta mass and ejecta velocity profiles obtained from numerical-relativity simulations. Specifically, we focus on the cases for "APR4Q3a75," "H4Q3a75," and "MS1Q3a75" in Tanaka et al. (2014) and "MS1Q7a75" in Hotokezaka et al. (2013). We plot the lightcurves predicted by our analytic model in Figure 4. Here, we set φ_ej = π, θ_ej = 1/5, v_min = 0.02c, and _th = 0.5 according to Kyutoku et al. (2015) and Korobkin et al. (2012). Using the result of the numerical simulations, (M_ej and v_max) are set to be (1 × 10⁻² M_⊙, 0.41c), (5 × 10⁻² M_⊙, 0.35c), (7 × 10⁻² M_⊙, 0.42c), and (7 × 10⁻² M_⊙, 0.51c) for "APR4Q3a75," "H4Q3a75," "MS1Q3a75," and "MS1Q7a75," respectively. These values of v_max are chosen so that mass-weighted root mean squares of the velocity distribution agree with the values of v_ch in Tanaka et al. (2014). As shown in Figure 3, we found that the lightcurves of the analytic model agree with the results of the radiation-transfer simulations within a factor of ∼1.4.

Zoom In Zoom Out Reset image size

The same heating rate and thermalization factor are adopted for both the analytic model and the radiation-transfer simulation. While the wavelength-dependent opacity is taken into account in the radiation-transfer simulations, we employ gray-opacity with κ = 10 cm² g⁻¹ in the analytic model, which is expected to be good approximation for the bolometric lightcurves for lanthanoid-rich ejecta, as shown in Kasen et al. (2013) and Tanaka & Hotokezaka (2013). A much more simplified morphology for the ejecta is used in the analytic model than in Tanaka et al. (2014). In particular, the parameters of the ejecta morphology, namely θ_ej and φ_ej, are fixed for different binary parameters. Kyutoku et al. (2015) and Kawaguchi et al. (2015) showed that the ejecta exhibit a similar shape in most cases, particularly if the ejecta mass is larger than 0.01 M_⊙. Kyutoku et al. (2015) also showed that there is the exception that φ_ej becomes as large as 2π with a substantial mass ≳0.01 M_⊙ (model ALF2-Q7a75). Moreover, θ_ej has a variation up to factor of ≲2 among the models with substantial mass ejection. Changes of θ_ej and φ_ej by a factor of 2 change L by ∼60% and ∼40%, respectively, and t_c by ∼40%.

2.3. Bolometric Correction

To predict the magnitude at each wavelength we need to impose an additional assumption because we cannot know the spectra from the analytic model derived in Section 2.2. However, the spectra of the kilonovae/macronovae are determined by very complex frequency-dependent radiative processes of lanthanoids and, hence, it is not easy to model them analytically. Thus, instead, we introduce a phenomenological approach to reproduce the spectra of kilonova/macronova models in the following.

We define a bolometric correction of a specific band (referred to as the X-band) ΔM_X (X = u, g, r, i, z, J, H, K) as

$$\begin{eqnarray}&&{\rm{\Delta }}{M}_{X}(t):= {M}_{{\rm{bol}}}(L(t))-{M}_{X}(t),\end{eqnarray} \tag{ 13 }$$

where M_X and M_bol are the X-band AB magnitude and the bolometric magnitude, respectively. As far as the photospheric emission dominates the total luminosity, the dominant factor determining the color temperature of the emission, or the bolometric correction, is the temperature near the photosphere. Thus, we need to obtain the temperature at the photosphere to estimate M_X(t). One possible estimator for the temperature at the photosphere is the effective temperature of the surface emission, which is given from Equation (12) by the Stefan–Boltzmann law, i.e.,

$$\begin{eqnarray}{T}_{{\rm{eff}}} & = & {\left(\displaystyle \frac{L(t)}{\sigma S(t)}\right)}^{1/4}\\ & \propto & \left\{\begin{array}{cc}{(t{/M}_{{\rm{ej}}}^{1/(2+2\alpha )})}^{-1-\alpha } & \;\;t\leqslant {t}_{{\rm{c}}},\\ {(t{/M}_{{\rm{ej}}}^{1/(2+\alpha )})}^{-2-\alpha } & \;\;t\gt {t}_{{\rm{c}}},\end{array}\right.\end{eqnarray} \tag{ 14 }$$

where S(t) ∝ t² is the surface area of the latitudinal edge.

Another simple guess for the temperature at the photosphere is the local temperature in the limit that the radiative cooling is negligible. In this limit, the local internal energy density is proportional to ${M}_{{\rm{ej}}}\dot{\epsilon }(t){t}^{-2}$, and assuming that the radiation pressure is dominant there, the local temperature becomes proportional to

$$\begin{eqnarray}&&T\propto {({M}_{{\rm{ej}}}{t}^{-2-\alpha })}^{1/4}={(t{/M}_{{\rm{ej}}}^{1/(2+\alpha )})}^{-(2+\alpha )/4}.\end{eqnarray} \tag{ 15 }$$

For both cases, the estimate of the temperature at the photosphere is written as a function ${t}^{\prime }=t/{M}_{{\rm{ej}}}^{1/n}$, where $n=2+\alpha$ or $2+2\alpha$. This suggests that the temperature at the photosphere and, thus, the bolometric corrections are approximately the same for all the models once the time is rescaled by ${t}^{\prime }=t/{M}_{{\rm{ej}}}^{1/n}$ with $n=2+\alpha$ or $2+2\alpha$, i.e., n = 3.2 or 4.4.

In fact, as shown in Figure 5, we found that the evolution of the bolometric corrections of each band filter obtained using the radiation-transfer simulations of Tanaka et al. (2014) agree approximately with each other by rescaling the elapsed time by $t/{M}_{{\rm{ej}}}^{1/3.2}$. This result shows that this scaling rule holds in the mass range M_ej = 0.01–0.07 M_⊙. This fact allows us to estimate the bolometric corrections for kilonovae/macronovae with arbitrary ejected mass from the results in each band with a single mass. However, here, the agreement of the bolometric correction is tested only with a few models with similar binary parameters, for example the spin parameter is fixed to 0.75. Thus, changing the parameter of the binary could introduce error in the agreement. The specific effects of binary parameters on the bolometric correction can only be checked by performing the radiation-transfer simulation while systematically varying the binary parameters, and we keep this for future work.

Zoom In Zoom Out Reset image size

As we can see from Equation (14), there is some ambiguity in choosing the value of n. Moreover, the emission from the photosphere may not be dominant after the system becomes optically thin. However, the time-rescaled bolometric corrections between models with different values of ${M}_{{\rm{ej}}}$ agree approximately with the numerical results in Tanaka et al. (2014), and choosing different value of n leads to only small differences (<0.5 mag) in the resulting band magnitudes in t ≤ 14 days as long as n = 3.2–4.4. Thus, in this paper, we employ the bolometric correction obtained by rescaling the result of a radiation-transfer simulation, specifically the result for APR4Q3a75 in Tanaka et al. (2014), by employing n = 3.2 to predict each band magnitude of the kilonova/macronova models. The bolometric corrections for the ugrizJHK-band magnitudes calculated from the results for APR4Q3a75 in Tanaka et al. (2014) are summarized in Table 3.

Table 3.  The Bolometric Corrections for ugrizJHK-band Magnitudes Calculated from the Results of the Model APR4Q3a75 in Tanaka et al. (2014)

Rescaled Time	Bolometric Correction ΔMX (mag) (X: Band Filter)
$\left(\tfrac{t}{{\rm{day}}}\right){\left(\tfrac{0.01{M}_{\odot }}{{M}_{{\rm{ej}}}}\right)}^{1/3.2}$	X =	u	g	r	i	z	J	H	K
1.5		−0.28	−0.45	−0.47	−0.71	−0.97	−1.61	−2.37	−4.55
2.0		−0.19	−0.57	−0.59	−0.32	−0.48	−0.07	1.18	0.41
2.5		0.25	−0.13	−0.20	0.05	−0.28	−1.50	0.86	2.22
3.0		0.93	0.93	0.60	0.56	0.37	−1.98	−4.29	−3.78
3.5		0.47	1.25	1.36	1.26	1.04	−1.38	−4.65	−6.36
4.0		−0.34	1.29	1.20	1.60	1.44	−0.74	−3.31	−6.12
4.5		0.03	0.74	1.15	1.51	1.65	−0.29	−2.33	−5.25
5.0		−0.39	0.37	0.89	0.99	1.76	0.42	−1.73	−4.10
5.5		−0.69	0.68	0.74	0.44	1.13	0.94	−0.69	−3.57
6.0		−1.21	0.41	0.74	0.50	1.15	1.05	−0.66	−3.55
6.5		−3.65	−0.58	1.05	0.94	1.28	1.66	−0.67	−3.76
7.0		−4.40	−0.87	−0.16	1.04	1.82	1.46	−0.07	−3.68
7.5		−1.74	−2.55	0.14	1.30	2.09	1.39	0.01	−3.24
8.0		−1.36	−2.87	0.15	1.31	2.11	1.38	0.04	−3.11
8.5		−3.75	−1.90	−0.35	1.16	2.26	1.53	0.15	−2.40
9.0		−4.70	−1.57	−0.57	1.08	2.29	1.59	0.21	−2.10
9.5		−3.93	−0.36	0.62	1.24	2.50	2.31	0.31	−1.34
10.0		−4.00	−0.44	0.56	1.20	2.47	2.29	0.32	−1.30
10.5		3.18	−0.52	0.84	1.34	2.74	2.43	0.30	−1.11
11.0		⋯	−0.60	1.97	1.92	3.79	2.93	0.20	−0.54
11.5		⋯	−0.68	1.89	1.87	3.73	2.90	0.21	−0.51
12.0		⋯	0.29	1.59	2.14	3.74	2.93	0.43	0.33
12.5		⋯	0.53	1.44	2.18	3.70	2.91	0.50	0.61
13.0		⋯	0.43	1.36	2.12	3.64	2.89	0.51	0.63
13.5		⋯	⋯	2.05	3.59	3.00	3.00	1.27	1.04
14.0		⋯	⋯	2.19	3.97	2.77	3.02	1.49	1.17

Download table as:  ASCIITypeset image

Figure 6 shows a comparison of the ugrizJHK-band magnitudes between the model obtained in Section 2.2 and the results of radiation-transfer simulations in Tanaka et al. (2014). The band magnitudes are calculated with a combination of the mass and the velocity fitting formulas, the bolometric luminosity model, and the bolometric correction model (we call this model "the analytic model" in the following). We set the hypothetical distance as 400 Mpc, which is the typical effective distance to a BH–NS merger event expected to be detected by gravitational-wave detectors (Abadie et al. 2010). For each band magnitude, we find that the analytic model always agrees with the previous results of the radiation-transfer simulations within ∼1 mag.⁶ We note that the variation in the ejecta geometry, and the uncertainties in the opacity and the heating rate may be sources of error in the lightcurve model and the radiation-transfer simulations, as we discuss in Section 5.

Zoom In Zoom Out Reset image size

Once the detection of gravitational waves from a compact binary merger is achieved, the chirp mass ${{ \mathcal M }}_{{\rm{chirp}}}={[{Q}^{3}/(1+Q)]}^{1/5}{M}_{{\rm{NS}}}$, the symmetric mass ratio $\nu =Q/{(1+Q)}^{2}$, and the effective spin parameter χ_eff = χ cos i_tilt of the binary will be estimated (e.g., O'Shaughnessy et al. 2014a, 2014b). Assuming a BH–NS merger event, the analytic model we constructed in this paper predicts the brightness and the duration of the kilonova/macronova for detected events. Tanaka et al. (2014) pointed out that observations in the i-band filter of wide-field 8 m class telescopes, such as Subaru/Hyper Suprime Cam (Miyazaki et al. 2006) and the Large Synoptic Survey Telescope (Ivezic et al. 2008), will be useful for follow-up observations of kilonovae/macronovae. It has also been pointed out that near-infrared observations by wide-field space telescopes, such as the Wide-Field Infrared Survey Telescope (WFIRST; Green et al. 2012), are promising. Thus, as an illustration, we show the evolutions of i- and H-band magnitudes, and the parameter dependence of BH–NS kilonovae/macronovae in this section.

In Figure 7, we plot the evolution of the i-band magnitude for several binary parameters. The label for the model denotes the EOS name, the mass ratio Q, and the effective spin parameter χ_eff of the binary. Specifically, "Q3," "Q5," and "Q6" denote the models with Q = 3, 5, and 6, respectively. "a0," "a5," and "a75" denote the models with χ_eff = 0, 0.5, and 0.75, respectively. We set M_NS = 1.35 M_⊙. We can see that the kilonova/macronova model becomes brighter as Q becomes smaller and χ_eff gets larger. This dependence reflects the fact that the ejecta mass increases with the decrease of the mass ratio of the binary and with the increase of the effective spin of the BH. We note that, exceptionally, the dependence of the luminosity on the mass ratio is not monotonic for the rapidly spinning case χ_eff ≳ 0.75, as seen in Figure 8. The kilonova/macronova model is always dimmer for APR4 than for H4, which also reflects the fact that the a NS with a larger radius produces more massive ejecta. This dependence of the ejecta mass has been shown by previous numerical simulations (e.g., Kawaguchi et al. 2015; Kyutoku et al. 2015).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 8, we show the contour plots of i- and H-band magnitudes at t = 3, 7, and 14 day after the merger for APR4 and H4 as a function of chirp mass ${{ \mathcal M }}_{{\rm{chirp}}}={[{Q}^{3}/(1+Q)]}^{1/5}{M}_{{\rm{NS}}}$ and ${\chi }_{{\rm{eff}}}=\chi \;{\rm{cos}}\;{i}_{{\rm{tilt}}}$. The magnitudes and the chirp masses are calculated assuming M_NS = 1.35 M_⊙. The distance is set to be 400 Mpc. We cut the regions from the plot in which the ejecta mass is below 10⁻³ M_⊙ because the data of the radiation-transfer simulation are not sufficient to calculate the bolometric correction 14 days after the merger for the case where the ejecta mass is below 10⁻³ M_⊙. We note that the i- and H-band magnitudes are always dimmer than 26 mag in those regions for both APR4 and H4.

For the case in which the NS EOS is H4, the i-band magnitude can reach ∼23 mag for a rapidly spinning BH, χ_eff ≥ 0.5 for a wide range of Q. The emission in the i-band with ≤26 mag can last for ∼7 day after the merger. On the other hand, a BH–NS merger with a small BH spin and a large chirp mass, particularly with ${\chi }_{{\rm{eff}}}\lesssim 0.67\;({{ \mathcal M }}_{{\rm{ch}}}/{M}_{\odot })-1.45$, produces no ejecta or ejecta with a mass below 10⁻³ M_⊙, and leads to a kilonova/macronova dimmer than ∼26 mag in the i-band. In contrast to the i-band magnitude, the H-band magnitude becomes brighter as time elapses. This time dependence reflects the fact that the color temperature decreases with time and, thus, the bolometric correction in the H-band increases with time. Due to this time dependence, the H-band magnitude is brighter than 25 mag until ∼14 day if the BH has a sufficiently large spin parameter, i.e., if ${\chi }_{{\rm{eff}}}\gtrsim 0.46\;({{ \mathcal M }}_{{\rm{ch}}}/{M}_{\odot })-0.72$.

For the case in which the NS EOS is APR4, ejecta heavier than 10⁻³ M_⊙ are produced when ${\chi }_{{\rm{eff}}}\gtrsim 0.23({{ \mathcal M }}_{{\rm{ch}}}/{M}_{\odot })-0.18$. There is no region where the i-band magnitude reaches 23 mag as long as χ_eff ≤ 0.9, and only a very narrow region with a very high BH spin (χ_eff ≳ 0.8), which can reach 24 mag only for the first ∼3 day. At 7 day after the merger, the i-band magnitude is dimmer than 26 mag in most of the region shown in the plot. The i-band magnitude reaches 26 mag for the case in which χ_eff ≳ 0.8 and ${{ \mathcal M }}_{{\rm{ch}}}\lesssim 2.8\;{M}_{\odot }$. Similarly, the H-band magnitude reaches 25 mag ∼7 day after the merger for the case in which χ_eff ≳ 0.85 and ${{ \mathcal M }}_{{\rm{ch}}}\lesssim 2.8\;{M}_{\odot }$, and is always dimmer than 25 mag for other cases.

Our results provide a guide for electromagnetic follow-up observations. In optical wavelengths, in order to maximize the possibility of detecting electromagnetic counterparts, follow-up observations using 8 m class telescopes within ∼3 day after the merger are crucial. A typical limiting magnitude for 8 m class telescopes is ∼26 mag for a 10 minute exposure. Therefore, rapid follow-up observations with 8 m class telescopes will open the possibility of detecting kilonova/macronova emission even for faint models with soft EOS (APR4).

At near-infrared wavelengths, the emission brighter than 25 mag lasts for t ∼ 14 day. Since a limiting magnitude for infrared telescopes, such as WFIRST, is ∼25 mag, infrared satellite observations will provide a greater chance of finding kilonovae/macronovae from BH–NS mergers.

The analytic model we introduced in Section 2 can also be used to constrain the binary parameters of kilonova/macronova emission. GRB 130603B is a gamma-ray burst which is possibly associated with a kilonova/macronova (Berger et al. 2013; Tanvir et al. 2013). Yang et al. (2015) and Jin et al. (2015) showed that GRB 060614 could also be a gamma-ray burst with a kilonova/macronova association. We apply our model to a possible kilonova/macronova candidate associated with GRB 130603B in this section.

A near-infrared excess from an afterglow model associated with GRB 130603B was observed by the Hubble Space Telescope (HST) in the H-band (Berger et al. 2013; Tanvir et al. 2013). According to their discussion, the redshift of the host galaxy is z = 0.356, and the near-infrared emission corresponds to an emission with a J-band absolute magnitude ${M}_{{\text{}}J,{\rm{abs}}}=-15.7\;{\rm{mag}}$ at ≈7 day after the prompt emission in the gamma-ray burst rest frame.

We calculate the J-band magnitudes at 7 day after the merger using the analytic model, ${M}_{{\text{}}J,{\rm{model}}}$, and plot the difference between these magnitudes and the observed magnitude ${\rm{\Delta }}{M}_{{\text{}}J}={M}_{{\text{}}J,{\rm{model}}}-{M}_{{\text{}}J,{\rm{abs}}}$ as a function of ${{ \mathcal M }}_{{\rm{chirp}}}$ and χ_eff in Figure 9. The magnitudes and the chirp masses are calculated assuming M_NS = 1.35 M_⊙. Four different cases employing APR4, ALF2, H4, and MS1 EOS are shown. The magnitude obtained by the analytic model is always larger than the observed magnitude. We note that our models tend to be fainter than the results of the radiation-transfer simulations, as we see Figure 6, and thus this comparison may be conservative.

Zoom In Zoom Out Reset image size

The blue and light-blue regions in Figure 9 show the regions where the differences of the J-band magnitude between the prediction and the observation are within 1 and 2 mag, respectively. Supposing that the error of the analytic model is within 1 mag, the binary parameters are constrained to this blue region. We note that the observational errors or uncertainties are smaller than the uncertainties in the model. The observational error of the HST H-band data is about 0.3 mag (Tanvir et al. 2013). The Galactic extinction in the direction to GRB 130603B is negligible (A_H ∼ 0.01 mag). The extinction in the host galaxy is estimated to be A_V ∼ 0.8–1.0 mag (de Ugarte Postigo et al. 2014; Fong et al. 2014) and, thus, according to the Galactic extinction law (A_J/A_V ∼ 0.282; Cardelli et al. 1989), the extinction in the J-band can be estimated as A_J ≈ 0.23–0.28 mag. Figure 9 shows that a BH–NS merger with a small spin χ_eff ≤ 0.3 is disfavored for all EOS cases. In particular, for the case in which the NS EOS is APR4, there is no region with χ_eff ≤ 0.9 in which the predicted magnitude is consistent with the observation within 1 mag. Thus, a NS with a small radius is disfavored for this observation (unless the BH is rapidly spinning with χ_eff > 0.9) if the ejecta are dominated by dynamical ejecta, which was also pointed out by Hotokezaka et al. (2013) using models in a smaller parameter region. Applying the same arguments to GRB 060614, we obtain similar results, and BH–NS mergers with a large NS radius, a small mass ratio, or a rapidly spinning BH are favored, which is consistent with the discussion in Yang et al. (2015).

We developed fitting formulas for the mass and the velocity of dynamical ejecta from a generic BH–NS merger based on recently performed numerical-relativity simulations. We combined these fitting formulas with a semi-analytic model for the BH–NS kilonova/macronova lightcurve that reproduces the results of radiation-transfer simulations with a small error. Specifically, the semi-analytic model reproduces the results of each band magnitude obtained by previous radiation-transfer simulations within ∼1 mag.

This semi-analytic model shows that a kilonova/macronova can be observed in optical wavelengths by 8 m class telescopes within 3 day from the merger if ${\chi }_{{\rm{eff}}}\gtrsim 0.23\;({{ \mathcal M }}_{{\rm{ch}}}/{M}_{\odot })-0.18$ for the case in which the NS EOS is stiffer than APR4, namely for the case in which the NS radius is larger than ∼11 km. On the other hand, if the optical-wavelength emission is observed by 8 m class telescopes for the case in which ${\chi }_{{\rm{eff}}}\lesssim 0.67({{ \mathcal M }}_{{\rm{ch}}}/{M}_{\odot })-1.45$, the origin of the emission may not be the kilonova/macronova from the dynamical ejecta but some other components unless the NS EOS is stiffer than H4. At near-infrared wavelengths, the follow-up observation of kilonova/macronova by space telescopes is useful because the emission can be bright and last long enough to be detected for the wide parameter space of a BH–NS merger. In particular, if the NS EOS is as stiff as H4, the emission can be detected even at ∼14 day after the merger for the case in which ${\chi }_{{\rm{eff}}}\gtrsim 0.46\;({{ \mathcal M }}_{{\rm{ch}}}/{M}_{\odot })-0.72$. We also applied our model to GRB 130603B as an illustration, and showed that a BH–NS merger with a rapidly spinning BH or a large NS radius is favored, which is consistent with previous studies (Hotokezaka et al. 2013).

There are several error sources for the results of the analytic model and the radiation-transfer simulations. Variation in ejecta geometry may be a source of error in the lightcurve model. It can change the lightcurve of the analytic model by ∼0.5 mag. Errors in the mass and the velocity fitting formulas can also induce errors in the lightcurves of ∼0.6 mag. Moreover, the uncertainties in the opacity and the heating rate can be additional factors which may affect the results. Since there are still some uncertainties in the opacity of lanthanoids, the real opacity could change from the one which is employed in both the analytic model and the simulations in Tanaka et al. (2014). For example, because the bolometric luminosity is proportional to κ^−1/2 in t < t_c, the increase of κ by a factor of 3 would make the bolometric magnitude larger by 0.6 mag. Uncertainties in the heating rate and the thermal efficiency can induce errors in the lightcurve linearly.

We should also note that there may exist some other ejecta components which are produced in the post-merger phase and could be additional sources of emission. For example, Kiuchi et al. (2015), Just et al. (2015), and Fernández et al. (2015) show that mass comparable to (or even more than) the dynamical ejecta may be ejected as a disk wind from a BH–NS merger remnant. While a screening effect of subsequent ejecta by preceding ejecta is argued for NS–NS mergers (Kasen et al. 2015), this may not take place for the kilonovae/macronovae from BH–NS mergers because the dynamical ejecta is confined in a narrow region around the equatorial plane; additional components of ejecta could enhance the electromagnetic emission. Thus, additional components may be visible and the emission may be brighter than we predict in this paper. This may make it difficult to constrain the binary parameters from kilonova/macronova observations. For example, additional components may allow a larger parameter region (which is consistent with GRB 130603B) to exist even for the case in which the NS EOS is as soft as APR4. Therefore, we should take the possibility of other emission into account to constrain the models, and we keep this as a future task.

We thank F. Foucart for helpful comments on the preprint version of the manuscript. This study was supported in part by the Grants-in-Aid for the Scientific Research of JSPS (14J02950, 24244028, 24740117, 15H02075, 15H06857), MEXT (25103515, 15H00788), and RIKEN iTHES project.