Toward Rapid Transient Identification and Characterization of Kilonovae

As pointed out, to perform accurate NR and radiative transfer simulations remains a challenging task, and further work including a better microphysical treatment is needed to allow a detailed understanding of ejecta, r-processes, and EM emission. In the following, we give a brief overview about some approaches without guarantee of completeness.

There are two major classes of kilonova models, those that model dynamic ejecta and those that model winds. Historically, kilonovae from dynamical ejecta were proposed earlier than those from wind kilonova. In this paper, we use parameterized models that concentrate on dynamic ejecta (Kawaguchi et al. 2016; Dietrich & Ujevic 2017), although the technique is generic enough to use wind models as well. It will be important when using this method on a transient to either include a wind model in addition to a dynamic ejecta model or, perhaps better, include one capable of incorporating both, such as the toy model in Metzger (2017).

One model is driven by the merger of two neutron stars, where material ejected during or following the merger assembles into heavy elements by the r-process (Metzger et al. 2010). Kulkarni (2005) postulated that ⁵⁶ Ni or free neutrons provided the power to these events, although ⁵⁶ Ni cannot be produced in the neutron-rich environments (Metzger 2017). Li & Paczynski (1998) first pointed out that radioactive ejecta from compact binary mergers are a source of electromagnetic emission. Metzger et al. (2010) used radioactive heating rates derived from r-process nuclear network calculations to determine the correct luminosity scale for the corresponding light curves. Kasen et al. (2013) then pointed out the high opacities of the resulting Lanthanides, with Barnes & Kasen (2013) and Tanaka & Hotokezaka (2013) performing the first realistic simulations for kilonova light curves.

EM emission then occurs during the radioactive decay of the resulting nuclei. Barnes et al. (2016) explore the emission profiles of the radioactive decay products, which include nonthermal β-particles, α-particles, fission fragments, and γ-rays, and the efficiency with which their kinetic energy is absorbed by the ejecta. By determining the net thermalization efficiency for each particle type and implementing the results into detailed radiation transport simulations, they provide kilonova light-curve predictions. Metzger et al. (2015) also explore the β-decay of the ejecta mass powering a "precursor" to the main kilonova emission, which peaks on a timescale of a few hours in the blue. Rosswog et al. (2017) use semianalytical models based on nuclear network simulations studying in detail the effect of the nuclear heating rate and ejecta electron fraction. The work of Rosswog et al. (2017) shows in detail how light-curve predictions change significantly for different nuclear physics parameters, e.g., the usage of different mass models.

Another kilonova mode, also driven by the r-process, in which radioactively powered transients are produced by accretion disk winds after the compact object merger was proposed by Kasen et al. (2015). In this model, the light curves contain two distinct components consisting of a ≈2-day blue optical transient and ≈10-day infrared transient. For this model, mergers resulting in a longer-lived neutron star or a more rapidly spinning black hole result in a brighter and bluer transient.

In addition to the numerical work, a handful of analytical models have also been developed with the purpose of approximating kilonova light curves. Based on NR simulations, Kawaguchi et al. (2016) derive fitting formulas for the mass and the velocity of ejecta from a generic BHNS merger and combine this with an analytic model of the kilonova light curve based on the radiative Monte Carlo (MC) simulations of Tanaka et al. (2014). Dietrich & Ujevic (2017) expand this work by using a large set of NR simulations to explore the EM signals from BNSs. The NR fit estimating the ejecta mass, velocity, and morphology is extended by an analytical model also based on the radiative MC simulations of Tanaka et al. (2014).

Parameterized models as proposed in Kawaguchi et al. (2016) and Dietrich & Ujevic (2017) directly tie GW parameters to expectations about the potential kilonova counterpart. They do not require NR and radiative transfer simulations to be completed, which is an impossible task over the few days of observations. Assumptions about the equation of state (EOS) of neutron stars, as well as measurement of the mass of the compact objects involved, allow the computation of the luminosity and light curves of kilonovae.

Because the ejecta morphology, the thermalization efficiency, and the opacity are not well constrained, it is advantageous to use a variety of models that estimate these quantities in different ways. In general, the luminosity will depend on the thickness of the ejecta, which is one of the main differences between BNS and BHNS systems. The thinner the ejecta becomes, the higher the density and temperature become. This affects the color temperature of the spectrum and consequently has a large impact on the detected light curve.

There are two limiting cases: (i) the ejecta are geometrically thick and approximately spherical, and (ii) the ejecta are geometrically thin. In general, due to shock-driven ejecta, BNS mergers correspond mostly to the former case and BHNS systems to the latter case; however, a clear distinction is impossible. The morphology of ejecta affects the diffusion timescale and changes the evolution of the light curve before the system becomes optically thin. When the system is optically thin, the difference in morphology may not be important for the light-curve evolution anymore. Since information about ejecta velocity is primarily contained in the light curve during the optically thick phase, modeling of this phase is important to constrain ejecta velocity.

As a first comparison between different models, we consider spherical ejecta with ${M}_{\mathrm{ej}}\approx 5\times {10}^{-3}$ and ${v}_{\mathrm{ej}}\approx 0.2$ (see Barnes et al. 2016). Here and in the following, we give ${v}_{\mathrm{ej}}$ in fractions of the speed of light and masses in fractions of the mass of the Sun M_⊙. For the nonspherical parameterized models of Kawaguchi et al. (2016) and Dietrich & Ujevic (2017), we further assume $\theta =0.2\,\mathrm{rad}$. From Rosswog et al. (2017), we include a model with ${M}_{\mathrm{ej}}=0.0079$ and ${v}_{\mathrm{ej}}=0.12$, which is closest to our fiducial model.

Additionally, we include the approximant of Metzger et al. (2015), which focused on the blue transient produced at a time around merger, which uses a neutron mass cut ${m}_{{\rm{n}}}={10}^{-4}$, opacity of $\kappa =30\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, and electron fraction ${Y}_{{\rm{e}}}=0.05$.

Figure 1 shows the bolometric luminosity and the light curves in the g (dashed) and i (solid) bands. The kilonova models have significant short-term dynamics, with changes of more than a magnitude in less than a day. Both the Kawaguchi et al. (2016) and Dietrich & Ujevic (2017) models are based on the MC simulations of Tanaka et al. (2014), for which a constant thermal efficiency is assumed (${\epsilon }_{\mathrm{th}}=0.5$).

Zoom In Zoom Out Reset image size

The model of Barnes et al. (2016) includes a time-dependent efficiency, which leads to a faster decay of the bolometric luminosity and magnitude because after a few days after the merger the thermalization efficiency drops below the constant thermalization efficiency employed in the Tanaka et al. (2014) simulations. Rosswog et al. (2017) employ both time-dependent and constant efficiencies and use a more complex density profile. The model picked from Rosswog et al. (2017) shows a smaller bolometric luminosity than other models; notice, however, that as shown in Rosswog et al. (2017), the usage of different mass models affects the luminosity by about $\approx 600 \%$, i.e., all presented models come with large uncertainties and are crucially dependent on nuclear physics assumptions. The model of Metzger et al. (2015) describes the blue transient arising from a small fraction of the ejected mass that expands sufficiently rapidly such that the neutrons are not captured and instead β-decay, giving rise to a clear peak in the bolometric luminosity visible around the time of merger.

Comparing Dietrich & Ujevic (2017) and Kawaguchi et al. (2016), we see a clear difference in the g band. This has already been pointed out in Tanaka et al. (2014). The main difference seems to arise from the difference of employed bolometric corrections, which itself will depend on the ejecta morphology. Since BHNS ejecta are much more nonspherical and are concentrated in the equatorial plane, they have higher temperatures that make the spectrum bluer than BNS ejecta with the same mass.

We can take the opportunity of having a variety of kilonova models accessible to compare the light-curve colors. It is common in dedicated searches for kilonovae to make color cuts (Doctor et al. 2017). Figure 2 shows the difference between the g and i bands for the models presented in Figure 1. As expected, all of the kilonova models show differences of at least 2 mag, especially on later timescales. For this reason, independent of the employed kilonova model, the proposed analysis will optimize the strategy for the detection of GW optical counterparts. Given the relative consistency in color among the models, imaging the transients in both the blue/green and the near-infrared can help differentiate from other transients. Due to the high opacities of Lanthanide elements, most models predict emission in the near-infrared wavelengths. These observations are required within the first few days owing to the faint magnitudes involved. As explained above, the significant changes in magnitude over day timescales can also help differentiate them as compared to possible background transients such as SNe Ia.

Zoom In Zoom Out Reset image size

As shown in Figure 1 (left panel), the bolometric luminosity of the models from Kawaguchi et al. (2016), Dietrich & Ujevic (2017), Barnes et al. (2016), and Rosswog et al. (2017) can be significantly different (we do not include the blue transient proposed in Metzger et al. [2015] in the following analysis since it is powered by a different mechanism). While similar ejecta masses, velocities, and energy deposition rates are employed, the models use different density profiles, morphology, and thermalization efficiency. The models of Kawaguchi et al. (2016) and Dietrich & Ujevic (2017) assume $\rho \propto {r}^{-2}$ for the density profile, nonspherical geometry, and a constant thermalization efficiency (${\epsilon }_{\mathrm{th}}\approx 0.5$). The model of Barnes et al. (2016) assumes spherical ejecta with $\rho \propto {r}^{-1}$ for the density profile, and the time-dependent (mass-dependent) thermalization efficiency is taken into account. The model of Rosswog et al. (2017) also assumes spherical ejecta with a homogeneously expanding density profile and time-dependent (mass-dependent) thermalization efficiency with the FRDM model.

To check how these differences affect the bolometric light curves, we perform a simple radiation transfer simulation varying the density profile, ejecta morphology, and thermalization efficiency. In this calculation, we assume the flux-limited diffusion approximation of the radiative transfer (Levermore & Pomraning 1981), a constant gray opacity with $10\,{\mathrm{cm}}^{2}\ {{\rm{g}}}^{-1}$, and the heating rate that is employed in Kawaguchi et al. (2016) and Dietrich & Ujevic (2017).

Figure 3 compares the bolometric luminosity for various setups. The figure clearly shows that different ejecta morphologies and thermalization efficiencies change the bolometric luminosity by a factor of $\approx 2$. This explains qualitatively the difference in the bolometric luminosity and light curves in Figure 1. The difference in the model of Rosswog et al. (2017) is also explained by the difference in the ejecta mass, ejecta velocity, and thermalization efficiencies. On the other hand, a different density profile has only a minor effect.

Zoom In Zoom Out Reset image size

These results indicate that for future development of analytical kilonova approximants the focus should be put on modeling the ejecta morphology and the time-dependent thermalization efficiency. We also find that considering a constant thermalization efficiency of ${\epsilon }_{\mathrm{th}}=1$ and then multiplying by ${\epsilon }_{\mathrm{th}}(t)$ (given in Barnes et al. 2016) or directly employing a time-dependent thermal efficiency leads only to differences of $\approx 40 \%$. This suggests that, at least for the bolometric luminosity, the time dependency of the thermalization efficiency can be approximately taken into account just by multiplying its function by the luminosity obtained by the constant efficiency. This is of particular importance for further improvement of the parameterized models, which, at the current stage, are based on simulations employing a constant thermalization efficiency.

In addition to the discussed effects, further uncertainties exist, which make a modeling and prediction of kilonova luminosities difficult. Rosswog et al. (2017) point out that the electron fraction and heating rate are main uncertainties in the current modeling of kilonova light curves. They find that by using two different mass models (DZ31, Duflo & Zuker 1995; Finite Range Droplet Model, Moller et al. 1995) the bolometric luminosity can be different up to $\approx 600 \%$. This is caused by the fact that the nuclear heating rate enters linearly into the bolometric luminosity.

As a first test of the numerical method, we use light curves produced by the parameterized models for BNS and BHNS and also recover the ejecta properties with the same models. The top row of Figure 4 shows light curves of such a comparison, where we assume uncertainties of the models of 1 mag as stated in Dietrich & Ujevic (2017) and Kawaguchi et al. (2016). The top row shows that the injected light curves are recovered properly. We quantify the level of overlap between parameters with "corner" plots (Foreman-Mackey 2016), shown in the bottom row of Figure 4. Shown are 1D and 2D posteriors marginalized over the rest of the parameters. In general, there are a few key features. First, with the $\approx 1$ mag uncertainty associated with these models, a large number of light curves computed with the parameterized models are consistent with the injected/baseline light curve we took. This means that no strict parameter constraints can be obtained. However, although the models have stated $\approx 1$ mag uncertainty, we can study a possible scenario with models having smaller uncertainties, e.g., $\approx 0.2$ mag or even 0.04 mag, which approximate the characteristic uncertainty for observations. In Figure 5, we show histograms for ${M}_{\mathrm{ej}}$ for the case where the uncertainties are varied. The figure demonstrates that ${M}_{\mathrm{ej}}$ constraints are significantly improved when the assigned error to the model is small. In particular, for an uncertainty of 0.2 mag the ejecta mass can be determined up to ${\mathrm{log}}_{10}{M}_{\mathrm{ej}}\approx \pm 0.5$, and in cases where the uncertainty would be limited by the observation (uncertainty of 0.04 mag) the ejecta mass could be determined to ${\mathrm{log}}_{10}{M}_{\mathrm{ej}}\,\approx \pm 0.1$. This motivates the need for further improved parameterized models of kilonova light curves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In contrast to the ejecta mass, the ejecta velocity is poorly constrained in our analysis. This is because the analytic models do not include times $t\lesssim 1\,\mathrm{day}$, where the ejecta are optically thick. However, the dependence on the ejecta velocity is only significant during this stage. Afterward, the light curves are primarily determined by the ejecta mass. Therefore, to improve the estimation of the ejecta velocity, extension of the light-curve models to earlier times is required.

We now perform a comparison between the parameterized models and results from Tanaka et al. (2014). In this and the following subsections, as we are using models that differ from the light curves they are being compared to, there will be a bias above and beyond the statistical uncertainty from sampling over the model parameter space. We will report the statistical uncertainty in the following sections and compare the statistical uncertainty to the true values to estimate bias. As discussed above, the models differ in the thermalization efficiency, opacity calculations, density profiles, and in other ways, and these differences lead to light-curve differences and therefore biases in the parameter estimation. Therefore, the models effectively have larger error bars than the 1 mag error stated in Kawaguchi et al. (2016) and Dietrich & Ujevic (2017).

For this analysis, we distinguish between BNS and BHNS. The BNS setups of Tanaka et al. (2014) are compared to the Dietrich & Ujevic (2017) model, and the BHNS light curves are compared to the Kawaguchi et al. (2016) model. In Figure 6, we show histograms for ${M}_{\mathrm{ej}}$ for uncertainties of 1 mag (dot-dashed lines). The ejecta mass corresponding to the light curves of Tanaka et al. (2014) (vertical dashed lines) is always within the posteriors of the models for the 1 mag posteriors (dot-dashed lines). We find that for 0.2 mag (solid lines) uncertainties some of the true values for BNS systems lie outside the estimated posteriors, which is to be expected because the uncertainties in Dietrich & Ujevic (2017) and Kawaguchi et al. (2016) are 1 mag. But, even for an assigned uncertainty of 0.2 mag, the posteriors of the BHNS setups are consistent with the injected values, which suggests that recovering smaller ejecta masses are in general less accurate (or equivalently have more bias). This might be caused by inaccuracies in the employed bolometric corrections and is already visible in Figure 9 of Dietrich & Ujevic (2017).

Zoom In Zoom Out Reset image size

We now perform a comparison between the parametric models and those of Barnes et al. (2016) and Rosswog et al. (2017). In Figure 7, we take the Barnes et al. (2016) (top panel) model rpft_m005_v2 and the NS12NS12 FRDM model of Rosswog et al. (2017) (bottom panel) and use the Dietrich & Ujevic (2017) model for recovery. One finds that the relative magnitudes between the bands are mostly consistent across the models. However, the models are not able to reproduce the light curves as accurately as for Tanaka et al. (2014). We find that multiple parameters cannot be constrained, while the ejecta mass, shown in Figure 8, is overestimated. Furthermore, for Rosswog et al. (2017) the parameter estimation pipeline leads to a T₀ estimate of the order of a few days, which suggests that follow-up searches using the current parameterized models would not correctly detect transients with light curves similar to those given in Rosswog et al. (2017).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The origin of the difference between the parameterized models and those of Barnes et al. (2016) and Rosswog et al. (2017) is that Dietrich & Ujevic (2017) was built using the light curves of Tanaka et al. (2014). It can be expected that parameterized models approximating the results of Barnes et al. (2016) and Rosswog et al. (2017) can be obtained as well. This shows that for future development, it is urgently required to provide light curves using full radiative transfer simulations that are as realistic as possible, i.e., including different ejecta components, time-dependent efficiency, and complex ejecta morphologies.

We also compare to a few non-kilonova models in Figure 9. Considering the different origin of the EM signal, we expect that the kilonova models cannot capture the injected light curves. We use the Metzger et al. (2015) fiducial model (top panel of Figure 9) describing the blue kilonovae precursor and an SN Ia model from Guy et al. (2007) (bottom panel). Metzger et al. (2015) do have an initially higher blue component. The best-fit curve from Dietrich & Ujevic (2017) is capable of producing time-dependent light-curve approximants. For the SN Ia it was not possible to compute time-dependent light curves with the parameterized models that approximate the SN Ia light curve. This shows that the parameterized models can also help to distinguish transients with different origins.

Zoom In Zoom Out Reset image size

In addition to the large uncertainty of the NR data, the models also contain degeneracies that do not allow the simultaneous extraction of all binary parameters. The ejecta mass and velocity as functions of the binary parameters for BHNS can be approximated by

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

$$\begin{eqnarray}&&{v}_{\mathrm{ej}}={b}_{1}\,q+{b}_{2},\end{eqnarray} \tag{ 2 }$$

with ${\chi }_{\mathrm{eff}}=\chi \,\cos \,{i}_{\mathrm{tilt}}$, where ${i}_{\mathrm{tilt}}$ is the angle between the dimensionless spin of the black hole χ and the orbital angular momentum and ${\tilde{r}}_{\mathrm{ISCO}}$ is the radius of the innermost stable circular orbit normalized by the black hole mass. ${a}_{1},{a}_{2},{a}_{3},{a}_{4},{n}_{1},{n}_{2},{b}_{1},{b}_{2}$ are fitting parameters that are determined by comparison to a large set of NR data (see Kawaguchi et al. 2016). For BNS setups, the ejecta properties are approximated by

$$\begin{eqnarray}{M}_{\mathrm{ej}}^{\mathrm{fit}} & = & {10}^{-3}\left\{\left[a{\left(\displaystyle \frac{{M}_{2}}{{M}_{1}}\right)}^{1/3}\left(\displaystyle \frac{1-2{C}_{1}}{{C}_{1}}\right)+b{\left(\displaystyle \frac{{M}_{2}}{{M}_{1}}\right)}^{n}\right.\right.\\ & & +\left.\left.c\left(1-\displaystyle \frac{{M}_{1}}{{M}_{1}^{* }}\right)\right]{M}_{1}^{* }+(1\leftrightarrow 2)+d\right\},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{v}_{\mathrm{ej}}=\sqrt{{v}_{\rho }^{2}+{v}_{z}^{2}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{v}_{\rho ,z}=\left[{a}_{\rho ,z}\left(\displaystyle \frac{{M}_{1}}{{M}_{2}}\right)(1+{c}_{\rho ,z}\ {C}_{1})\right]+(1\leftrightarrow 2)+{b}_{\rho ,z},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\theta }_{\mathrm{ej}}=\displaystyle \frac{{2}^{4/3}{v}_{\rho }^{2}-{2}^{2/3}{({v}_{\rho }^{2}(3{v}_{z}+\sqrt{9{v}_{z}^{2}+4{v}_{\rho }^{2}}))}^{2/3}}{{({v}_{\rho }^{5}(3{v}_{z}+\sqrt{9{v}_{z}^{2}+4{v}_{\rho }^{2}}))}^{1/3}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\phi }_{\mathrm{ej}}=4{\theta }_{\mathrm{ej}}+\displaystyle \frac{\pi }{2},\end{eqnarray} \tag{ 7 }$$

with the fitting parameters $a,b,c,d,{a}_{\rho },{a}_{z},{b}_{\rho },{b}_{z},{c}_{\rho },{c}_{z},n$ given in Dietrich & Ujevic (2017).

As can be concluded from Equations (1)–(7), the BHNS model depends on the mass ratio q, the "effective" spin of the black hole ${\chi }_{\mathrm{eff}}$, the baryonic mass of the neutron star ${M}_{\mathrm{NS}}^{* }$, the quotient of the neutron star's gravitational mass ${M}_{\mathrm{NS}}$ and baryonic mass ${M}_{\mathrm{NS}}^{* }$, and its compactness C, i.e., five parameters. For the case of BNS systems, the number increases to six: the gravitational masses ${M}_{1},{M}_{2}$, the baryonic masses ${M}_{1}^{* },{M}_{2}^{* }$, and the compactnesses ${C}_{1},{C}_{2}$ of the neutron stars.

As an example to visualize possible degeneracies in Equations (1)–(7), let us suppose that the ejecta mass and the ejecta velocity were measured for a BHNS setup. In Figure 10, we show as red surfaces the allowed binary parameters for which ${M}_{\mathrm{ej}}={10}^{-2}$ under the assumptions of $C=0.13,0.15,0.17$. In addition, we make use of the quasi-universal relation given by Equation (8) (see discussion in the next subsection) to connect the gravitational and baryonic mass to the compactness. As a blue surface, we mark the binary parameters for which ${v}_{\mathrm{ej}}=0.28$. According to Equation (2), the measurement of ${v}_{\mathrm{ej}}$ would determine the mass ratio of the system q but leave the other parameters unconstrained.

Zoom In Zoom Out Reset image size

Figure 10 shows that even if ${M}_{\mathrm{ej}}$ and ${v}_{\mathrm{ej}}$ are accurately known, the binary parameters cannot be determined. The intersections between the red and blue surfaces mark all the allowed regions for which the ejecta properties are consistent with the estimated ${M}_{\mathrm{ej}},{v}_{\mathrm{ej}}$ under the assumption of a given compactness C. Consequently, an accurate measurement of the binary properties is only possible for cases for which more parameters than ${M}_{\mathrm{ej}},{v}_{\mathrm{ej}}$ are determined, e.g., ${\theta }_{\mathrm{ej}}$ and ${\phi }_{\mathrm{ej}}$, or for cases where, due to a simultaneous detection of GWs, some binary parameters are known.

Due to the large number of unknown binary parameters in Equations (1)–(7), several degeneracies exist and binary parameters cannot be constrained uniquely. To reduce this effect, we substitute some parameters with the help of quasi-universal relations. Quasi-universal properties for single neutron stars have been first found by Yagi & Yunes (2013) and were consequently studied for a variety of parameters (see, e.g., Maselli et al. 2013; Pappas & Apostolatos 2014; Yagi et al. 2014); even in BNS systems quasi-universal relations are present (e.g., Bernuzzi et al. 2014). We propose a relation between the quotient of baryonic and gravitational mass ${M}^{* }/M$ and the compactness C of a single neutron star. To construct this relation, we use the EOSs employed for the data set studied in Dietrich & Ujevic (2017) but only consider EOSs that allow nonrotating NS masses above 1.9, which lies even below the highest measured NS mass of $\approx 2.01$. Figure 11 shows ${M}^{* }/M$ as a function of the compactness C for all EOSs. We find that only a small spread is caused by the EOSs. Except for GlendNH3, all curves stay close together.

Zoom In Zoom Out Reset image size

We fit all data with an approximant of the form

$$\begin{eqnarray}&&\displaystyle \frac{{M}^{* }}{M}=1+a\,{C}^{n},\end{eqnarray} \tag{ 8 }$$

where the free fitting parameters are a = 0.8858 and n = 1.2082. The fit is included as a black dashed line in the top panel of Figure 11. By construction, we obtain for $C\to 0$ the correct limit of ${M}^{* }/M\to 1$. The residuals of the fit are shown in the bottom panel of Figure 11. Absolute errors within the compactness interval of $C\in [0.05,0.24]$ are within ±0.01, except for GlendNH3. This leads to fractional errors of $\lesssim 10 \%$ for the term $1-{M}^{* }/M$, which enters directly in the ejecta mass computation for BHNS and BNS systems. On average, fractional errors are $\lesssim 3 \%$. Considering the large uncertainty of Equations (1)–(7), we expect that the error caused by Equation (8) is negligible. But by introducing this relation, the number of free parameters for the BHNS model is reduced by one and for the BNS model is reduced by two. This allows for significantly better extraction of the binary parameters from the ejecta properties.

In the following, we use a similar scheme to that in Section 4 to explore how binary parameters can be recovered from a kilonova detection. We explore the situation where we have made a measurement of ${M}_{\mathrm{ej}}$ and ${v}_{\mathrm{ej}}$. We calculate the likelihood using a kernel density estimator commonly used in GW data analysis (Singer et al. 2014). This technique is useful for cases where the measurements of those distributions arise from parameter estimation with potentially highly correlated estimates among the variables, as is common in GW data analysis. The priors used in the analyses are as follows: for Kawaguchi et al. (2016), $3\leqslant q\leqslant 9$, $0\leqslant {\chi }_{\mathrm{eff}}\leqslant 0.75$, $1\leqslant {M}_{\mathrm{NS}}\leqslant 3$, and $0.1\leqslant C\leqslant 0.2$, while for Dietrich & Ujevic (2017), $1\leqslant {M}_{1}\leqslant 3$, $1\leqslant {M}_{2}\leqslant 3$, $0.08\leqslant {C}_{1}\leqslant 0.24$, and $0.08\leqslant {C}_{2}\leqslant 0.24$. The differences in compactness prior ranges are due to the differences in compactness used in the simulations the models used. The priors are flat over the stated ranges. For this reason, significant structure in the 1D and 2D contours arises from the posterior.

To begin, we explore the correlation between the variables by employing the very optimistic assumption of 1% Gaussian error bars on the measurement, which essentially inverts the equations in the previous section. We show in Figure 12 the parameters consistent with two different choices of ${M}_{\mathrm{ej}}$ and ${v}_{\mathrm{ej}}$. For the BNS case (left panel), we choose ${M}_{\mathrm{ej}}=5\times {10}^{-3}$ and ${v}_{\mathrm{ej}}=0.25$; for the BHNS case (right panel), we choose ${M}_{\mathrm{ej}}=5\times {10}^{-2}$ and ${v}_{\mathrm{ej}}=0.25$. In general, for the BNS systems, the constraints are not strong given the relatively wide variety of parameters that support nonzero ejecta masses and velocities. We choose to plot mass ratio ($q={M}_{1}/{M}_{2}$) and chirp mass (${M}_{{\rm{c}}}={({M}_{1}{M}_{2})}^{3/5}{({M}_{1}+{M}_{2})}^{-1/5}$) instead of M₁ and M₂, due to the clearer peaks in this parameterization. We clearly see in the 2D corner plots degeneracies between ${M}_{{\rm{c}}}$ and q, as well as between C₁ and C₂, which are similar to those described in the previous subsection. These indicate the fundamental limitations of EM-only observations in the measurements of these quantities. For the BHNS systems, the main constraint is on q, which has some correlation with compactness. Due to the significant correlations between q, ${\chi }_{\mathrm{eff}}$, and C, it will be difficult to constrain those individual parameters without measurements from other quantities. The degeneracies can be broken by using assumptions based on priors based on the measured masses of BNS systems, the spins of binary black holes, or the compactness of neutron stars, or indeed, these quantities as measured by coincident GW observations.

Zoom In Zoom Out Reset image size

Figure 13 shows more realistic levels of parameter estimates using the ${M}_{\mathrm{ej}}$ and ${v}_{\mathrm{ej}}$ contours sampled from a light curve with ${M}_{\mathrm{ej}}\approx 5\times {10}^{-3}$ and ${v}_{\mathrm{ej}}\approx 0.2$ with model uncertainties of 0.2 mag. The main difference between these results and the optimistic assumptions above are the relatively poor constraints on ${v}_{\mathrm{ej}}$. For the BNS system (left panel), because the constraints on mass ratio are tied to ${v}_{\mathrm{ej}}$, most values of mass ratio are allowed in this particular case. There are only minimal constraints on ${M}_{{\rm{c}}}$, C₁, and C₂. For the BHNS system (right panel), the only structure visible is the correlation between q, ${\chi }_{\mathrm{eff}}$, and C. In the case of precise measurement of the mass ratio and effective spin by GW parameter estimation, constraints on the neutron star compactness of $C\pm 0.2$ are possible.

Zoom In Zoom Out Reset image size

As a final comparison, we perform parameter estimation for the Dietrich & Ujevic (2017) model with 0.2 mag uncertainty, but instead of sampling in ${M}_{\mathrm{ej}}$ and ${v}_{\mathrm{ej}}$, we sample directly in the system parameters making use of Equations (3)–(7). Figure 14 shows the corner plots for this scenario. We find that the individual binary parameters are almost undetermined; only in the 2D M₁–M₂ or, as shown in the figure, the ${M}_{{\rm{c}}}$–q plane is a clear contour visible. According to the 1D posteriors of q, it seems that high mass ratios are ruled out. Additionally, ${C}_{1},{C}_{2}$ are almost unconstrained, but there seems to be a small preference for larger compactnesses for the shown example.

Zoom In Zoom Out Reset image size

Although we have only discussed the extraction of binary parameters for BNS configurations, similar results are obtained for BHNS systems.