Kilonova Emission from Black Hole–Neutron Star Mergers. I. Viewing-angle-dependent Lightcurves

It has long been proposed that mergers of binary neutron stars (BNS) and black holes and neutron stars (BH–NS) are progenitors of short-duration gamma-ray bursts (sGRBs; Paczynski 1986, 1991; Eichler et al. 1989; Narayan et al. 1992). The rapid accretion of a centrifugally supported disk/torus by the merger remnant (likely a BH or a rapidly rotating, highly magnetized NS) would drive a pair of collimated relativistic jets, which can be observed as an sGRB if the jet points toward Earth (Rezzolla et al. 2011; Paschalidis et al. 2015; Ruiz et al. 2016).⁷ The interaction of the relativistic jets with the surrounding interstellar medium would generate afterglows with emission ranging from radio to X-rays (Rees & Meszaros 1992; Meszaros & Rees 1993; Paczynski & Rhoads 1993; Mészáros & Rees 1997; Sari et al. 1998), which have been confirmed observationally (Berger et al. 2005; Fox et al. 2005; Hjorth et al. 2005).

Besides collimated relativistic jets, BNS and BH–NS mergers are expected to release an amount of neutron-rich matter (Lattimer & Schramm 1974, 1976; Symbalisty & Schramm 1982) that allows elements heavier than iron via the rapid neutron-capture process (r-process) to be synthesized. Li & Paczyński (1998) first predicted a type of transient event powered by the radioactive decay of r-process nuclei following a BNS or BH–NS merger. The more detailed research by Metzger et al. (2010) showed that the peak luminosity of such an optical–infrared transient is around several times 10⁴¹ erg s⁻¹ and hence, named it a "kilonova" (i.e., the luminosity is a factor of ∼10³ higher than a typical nova). Later, the characteristics of kilonova emission were widely investigated theoretically (Kulkarni 2005; Roberts et al. 2011; Barnes & Kasen 2013; Kasen et al. 2013; Tanaka & Hotokezaka 2013; Yu et al. 2013; Grossman et al. 2014; Metzger & Fernández 2014; Metzger & Piro 2014; Perego et al. 2014; Wanajo et al. 2014; Just et al. 2015; Kasen et al. 2015, 2017; Martin et al. 2015; Fernandez & Metzger 2016; Metzger 2017). In the past few years, several kilonova candidates following sGRBs have been claimed from optical–infrared emission in excess of the afterglow emission (Berger et al. 2013; Tanvir et al. 2013; Gao et al. 2015, 2017; Jin et al. 2015, 2016, 2020; Yang et al. 2015; Gompertz et al. 2018; Ascenzi et al. 2019; Rossi et al. 2020). However, due to scarce and ambiguous observational data and the lack of smoking-gun evidence that sGRBs are indeed from BNS or BH–NS mergers, these cases cannot be firmly confirmed. In general, for sGRB-related kilonovae (presumably on-axis mergers), the putative kilonova emission could be outshone by the luminous afterglow. Also, the usual emission timescale of kilonovae is only about a few days to a few weeks. Therefore, for blind optical transient surveys, it is easy to miss the optimal observational time window for searching kilonovae. On the other hand, BNS and BH–NS mergers are gravitational wave (GW) events that could be detected by the Advanced LIGO out to distances of ∼300 and ∼650 Mpc, respectively (Cutler & Thorne 2002). An optimal searching strategy for kilonovae is to take advantage of GW triggers for electromagnetic (EM) follow-up observations (e.g., Metzger & Berger 2012; Cowperthwaite & Berger 2015; Gehrels et al. 2016; Cowperthwaite et al. 2020).

On 2017 August 17, the first BNS merger GW source GW 170817 was detected by the LIGO/Virgo Collaborations (Abbott et al. 2017a). At Δt ∼ 1.7 s after the merger, an sGRB lasting ∼2 s, GRB 170817A, triggered the Fermi Gamma-ray Burst Monitor (GBM; Abbott et al. 2017c; Goldstein et al. 2017; Zhang et al. 2018). About 11 hr later, an associated ultraviolet–optical–infrared kilonova transient, named AT 2017gfo, was discovered in the galaxy NGC 4993, ∼40 Mpc away (Abbott et al. 2017b; Andreoni et al. 2017; Arcavi et al. 2017; Chornock et al. 2017; Coulter et al. 2017; Covino et al. 2017; Cowperthwaite et al. 2017; Díaz et al. 2017; Drout et al. 2017; Hu et al. 2017; Kasliwal et al. 2017; Kilpatrick et al. 2017; Lipunov et al. 2017; McCully et al. 2017; Nicholl et al. 2017; Pian et al. 2017; Pozanenko et al. 2018; Shappee et al. 2017; Smartt et al. 2017; Soares-Santos et al. 2017; Tanvir et al. 2017; Troja et al. 2017; Utsumi et al. 2017; Valenti et al. 2017). At the location coincident with this kilonova, a broadband afterglow source from radio to X-rays was detected, which is consistent with the standard synchrotron afterglow model with an off-axis viewing angle (Alexander et al. 2017; Haggard et al. 2017; Hallinan et al. 2017; Margutti et al. 2017; Troja et al. 2017; D'Avanzo et al. 2018; Dobie et al. 2018; Lazzati et al. 2018; Lyman et al. 2018; Ghirlanda et al. 2019). The discovery by LIGO/Virgo of this GW event from a BNS merger and the subsequent observations of EM counterparts set a milestone for the era of GW-led multimessenger astronomy. The associations between GW 170817, GRB 170817A, and AT 2017gfo provided the smoking-gun evidence for the BNS merger origin of sGRBs and confirmed the theoretical kilonova prediction.

The comprehensive observations of AT 2017gfo showed that its early- and late-stage lightcurve cannot be explained by one single radioactivity-powered component (however, see Waxman et al. 2018). In the literature, the data are widely interpreted by invoking two or even three different radioactivity-powered components (Cowperthwaite et al. 2017; Perego et al. 2017; Tanaka et al. 2017, 2018; Villar et al. 2017; Kawaguchi et al. 2018; Wanajo 2018; Wu et al. 2019). More specifically, the early-stage emission of AT 2017gfo was explained by a lanthanide-free "blue" component (with low opacity) while the late-stage emission is interpreted as a lanthanide-rich "red" component (with high opacity). Introducing both components could account for the AT 2017gfo emission evolving from a blue-component-dominated stage to a red-component-dominated stage. In addition, some authors (e.g., Perego et al. 2017 and Villar et al. 2017) found that an intermediate-opacity "purple" component may be needed to fully account for the broadband kilonova emission data of AT 2017gfo. Theoretically, the lanthanide-free "blue" component is thought to be produced due to heating by the shocks at the contact interface between two merging NSs (e.g., Oechslin & Janka 2006; Wanajo et al. 2014; Radice et al. 2016; Sekiguchi et al. 2016) or by neutrino irradiation from the remnant hypermassive NS (e.g., Metzger & Fernández 2014; Perego et al. 2014; Yu et al. 2018), while the tidal dynamical ejecta is usually interpreted as the "red" component. For the "purple" component, viscous heating and angular momentum transport of the remnant disk could form such intermediate-opacity ejecta (e.g., Fernández & Metzger 2013; Just et al. 2015; Siegel & Metzger 2017; Fujibayashi et al. 2020). One issue with this mainstream interpretation is that the blue component is too bright and too early unless an unreasonably small opacity is introduced (Li et al. 2018). Li et al. (2018; see also Yu et al. 2018; Ren et al. 2019) argued that this may point toward a long-lived central engine that continuously injects energy into the ejecta (Yu et al. 2013).

Besides BNS mergers, BH–NS mergers could be another type of source producing kilonova transients that astronomers expect to detect. During the third observing run (O3) of the LIGO-Virgo Collaboration, several GW events from BH–NS or MassGap merger candidates with low false-alarm rates have been discovered, e.g., S190814bv, S190924h, S190930s, and S200115j (LIGO Scientific Collaboration & Virgo Collaboration 2019a, 2019b, 2019c, 2020). Among them, the BH–NS merger candidate S190814bv has been widely discussed (e.g., Andreoni et al. 2019; Dobie et al. 2019; Gomez et al. 2019; Ackley et al. 2020; Kawaguchi et al. 2020a). Further, there was a BNS merger event (GW 190425) with total mass ∼3.4M_⊙. The possibility that the system is a BH–NS merger cannot be ruled out from GW data (Han et al. 2020; Kyutoku et al. 2020; The LIGO Scientific Collaboration et al. 2020). Mostly because of the incomplete coverage of the error boxes of the GW events, so far follow-up observations have not detected any EM counterpart, in particular, the associated kilonova emission from any of these events.⁸

There are two types of BH–NS mergers: one is the case where the NS directly plunges into the BH without being tidally disrupted (e.g., Shibata et al. 2009), while the other is where the NS undergoes tidal disruption before the merger so that an amount of matter will remain outside the BH after the merger. No sGRB or kilonova is expected in the former case.⁹ Only the latter case (relatively small mass ratio between BH and the NS) is interesting for sGRB and kilonova follow-up observations. Numerical relativity (NR) simulations revealed that the dynamical ejecta from these systems, caused by tidal forces, are highly anisotropic, showing a crescent shape (Kyutoku et al. 2013, 2015; Kawaguchi et al. 2015; Brege et al. 2018). Moreover, the wind outflows from the disk around the remnant BH are directional (e.g., Just et al. 2015; Wu et al. 2016; Siegel & Metzger 2017). As a result, the outflowing materials that power the kilonova emission are highly anisotropic in the BH–NS cases. Therefore, it is of great interest to explore the viewing angle effect on the kilonova lightcurves for these systems. On the other hand, due to the lack of shock heating and neutrino irradiation during or shortly after the merger, only a small fraction (e.g., a few percent Just et al. 2015) of remnant disk can be transformed into lanthanide-free "blue" component ejecta. Kilonovae from the BH–NS mergers would be obviously different from those from BNS mergers (e.g., AT 2017gfo).

In this paper, we study the lightcurves of BH–NS merger kilonovae in detail, paying special attention to the viewing angle effect. This paper is organized as follows. In Section 2, we collect 66 results from recent NR simulations for BH–NS mergers and extrapolate fitting formulae for tidal dynamical ejecta mass and velocity across a wide range of mass ratios. In Section 3, we model the dynamics and temperature evolution of each radioactivity-powered component for a BH–NS merger. In Section 4, we present our simulation results for the BH–NS merger kilonova emission and compare them with the data of AT 2017gfo. In Section 5, we briefly discuss the relationship between kilonovae and GRB afterglows for the on-axis cases. The discussion and conclusions are presented in Section 6. In Appendix A, we summarize the definitions of frequently used variables. In Appendix B, we provide a numerical method to model photosphere evolution and the predicted lightcurves as seen by observers from different lines of sight. In Appendix C, we briefly introduce the sGRB afterglow model.

In this section, we collect 66 simulation results of BH–NS mergers from the published literature and extrapolate the fitting formulae for the mass of tidal dynamical ejecta across a wide range of mass ratios. We also unify the velocity of the dynamical ejecta from different papers and give a new fitting formula covering the range from the near-equal-mass regime to the large-mass-ratio regime. The fitting formulae for the mass and velocity of dynamical ejecta are used to discuss viewing-angle-dependent kilonova lightcurves later.

In a BH–NS merger system, whether the NS directly plunges into the BH or is tidally disrupted by the BH can be described by the relative positions of the radius of the innermost stable circular orbit (ISCO) R_ISCO and the location of the radius at which the tidal disruption occurs, R_tidal (see Shibata & Taniguchi 2011; Shibata & Hotokezaka 2019 for a review). The normalized ISCO radius ${\widetilde{R}}_{\mathrm{ISCO}}={c}^{2}{R}_{\mathrm{ISCO}}/{{GM}}_{\mathrm{BH}}$, determined by the BH mass and the spin of the BH, is given by Bardeen et al. (1972), i.e.,

$$\begin{eqnarray}&&{\widetilde{R}}_{\mathrm{ISCO}}=3+{Z}_{2}-\mathrm{sign}({\chi }_{\mathrm{BH}})\sqrt{(3-{Z}_{1})(3+{Z}_{1}+2{Z}_{2})},\end{eqnarray} \tag{ 1 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{Z}_{1} & = & 1+{\left(1-{\chi }_{\mathrm{BH}}^{2}\right)}^{1/3}\left[{\left(1+{\chi }_{\mathrm{BH}}\right)}^{1/3}+{\left(1-{\chi }_{\mathrm{BH}}\right)}^{1/3}\right],\\ {Z}_{2} & = & \sqrt{3{\chi }_{\mathrm{BH}}^{2}+{Z}_{1}^{2}},\end{array}\end{eqnarray} \tag{ 2 }$$

where M_BH is the BH mass and χ_BH is the dimensionless spin parameter of the BH. On the other hand, the radius R_tidal at which tidal disruption occurs can be estimated by balancing the tidal force and the self-gravitational force at the surface of the NS. In Newtonian theory, it reads

$$\begin{eqnarray}&&{R}_{\mathrm{tidal}}\sim {R}_{\mathrm{NS}}{\left(\displaystyle \frac{3{M}_{\mathrm{BH}}}{{M}_{\mathrm{NS}}}\right)}^{1/3}.\end{eqnarray} \tag{ 3 }$$

If R_tidal ≲ R_ISCO, the NS would plunge into the BH without material outside the remaining BH; if R_tidal ≳ R_ISCO, the NS would be tidally disrupted by the BH while forming an accretion disk around the BH and an unbound tidal dynamical ejecta.

Foucart et al. (2018) presented a nonlinear model wherein the mass outside the BH is determined by the relative locations of R_ISCO and R_tidal. By considering 75 NR simulation results, the fitting model of the total remnant mass outside the remnant BH given by Foucart et al. (2018) is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{M}_{\mathrm{total},\mathrm{fit}}}{{M}_{\mathrm{NS}}^{{\rm{b}}}} & = & \max \left[\left(0.406\displaystyle \frac{1-2{C}_{\mathrm{NS}}}{{\eta }^{1/3}}\right.\right.\\ & & -\ \left.{\left.0.139{\widetilde{R}}_{\mathrm{ISCO}}\displaystyle \frac{{C}_{\mathrm{NS}}}{\eta }+0.255\right)}^{1.761},0\right],\end{array}\end{eqnarray} \tag{ 4 }$$

where $\eta =Q/{(1+Q)}^{2}$, ${M}_{\mathrm{NS}}^{{\rm{b}}}$ is the baryonic mass of the NS, Q = M_BH/M_NS is the mass ratio between the BH mass and the NS mass, and ${C}_{\mathrm{NS}}={{GM}}_{\mathrm{NS}}/{c}^{2}{R}_{\mathrm{NS}}$ is the compactness of the NS. This formula covers the range $Q\in [1,7]$, ${\chi }_{\mathrm{BH}}\in [-0.5,0.9]$, and ${C}_{\mathrm{NS}}\in [0.13,0.182]$.

As for the dynamical ejecta mass, Kawaguchi et al. (2016) presented a fitting model (by referring to Foucart 2012),

$$\begin{eqnarray}&&\begin{array}{r}\begin{array}{r}\begin{array}{l}\displaystyle \frac{{M}_{{\rm{d}},{\rm{f}}{\rm{i}}{\rm{t}}}}{{M}_{{\rm{N}}{\rm{S}}}^{{\rm{b}}}}=max\left[4.464\times {10}^{-2}{Q}^{0.250}\displaystyle \frac{1-2{C}_{{\rm{N}}{\rm{S}}}}{{C}_{{\rm{N}}{\rm{S}}}}-2.269\right.\\ \,\times \,\left.{10}^{-3}{Q}^{1.352}{\widetilde{R}}_{{\rm{I}}{\rm{S}}{\rm{C}}{\rm{O}}}+2.431\left(1-\displaystyle \frac{{M}_{{\rm{N}}{\rm{S}}}}{{M}_{{\rm{N}}{\rm{S}}}^{{\rm{b}}}}\right)-0.4159,0\right],\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 5 }$$

considering 45 NR simulation results. The ranges of the parameters are $Q\in [3,7]$, ${\chi }_{\mathrm{BH}}\in [0,0.9]$, and ${C}_{\mathrm{NS}}\,\in [0.13,0.18]$. Recently, some new NR simulation results of BH–NS mergers including the case with the near-equal-mass regime (e.g., Foucart et al. 2019) and high-spin BH regime (Lovelace et al. 2013) have been published. Here we attempt to fit the 66 NR simulation results we collected (see Table 1) using a similar form to the formula in Equation (4). The numerical models cover the parameter range with Q ∈ [1, 7], χ_BH ∈ [0, 0.97], and C_NS ∈ [0.108, 0.18]. We assume the resulting error estimate combines a 10% relative error and 0.005M_⊙ absolute error similar to Kawaguchi et al. (2016)¹⁰ :

$$\begin{eqnarray}&&{\sigma }_{{\rm{d}},{\rm{N}}{\rm{R}}}=\sqrt{{\left(0.1{M}_{{\rm{d}},{\rm{N}}{\rm{R}}}\right)}^{2}+{\left(0.005{M}_{\odot }\right)}^{2}}.\end{eqnarray} \tag{ 6 }$$

The best least-squares method fitting result is

$$\begin{eqnarray}\begin{array}{r}\begin{array}{r}\begin{array}{rcl}\displaystyle \frac{{M}_{{\rm{d}},{\rm{f}}{\rm{i}}{\rm{t}}}}{{M}_{{\rm{N}}{\rm{S}}}^{{\rm{b}}}} & = & max\left[\left(0.273\displaystyle \frac{1-2{C}_{{\rm{N}}{\rm{S}}}}{{\eta }^{1/3}}\right.\right.\\ & & -\,\left.{\left.0.035{\widetilde{R}}_{{\rm{I}}{\rm{S}}{\rm{C}}{\rm{O}}}\displaystyle \frac{{C}_{{\rm{N}}{\rm{S}}}}{\eta }-0.153\right)}^{1.491},0\right].\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 7 }$$

Figure 1 shows the NR simulation data of the dynamical ejecta mass versus our fitting results and the normalized residuals. One can see that three points have relatively large dynamical ejecta mass deviating from the fitting result. The deviant data with the largest dynamical ejecta mass (No. 21 in Table 1) are from Lovelace et al. (2013). Due to its large estimated numerical error, it is still located within the region of the 1 − σ confidence interval. The other two points of deviation (No. 53 and No. 54 in Table 1) are from Foucart et al. (2014). We note that similar initial conditions have also been studied by Kyutoku et al. (2015) and Brege et al. (2018), who obtained results consistent with our fitting results. Removing these two points to refit the NR simulation data of the dynamical ejecta, we replot the NR data of the dynamical ejecta versus fitting results in Figure 2 from 0 to 0.1 M_⊙. The best least-squares method fitting result now reads

$$\begin{eqnarray}\begin{array}{r}\begin{array}{r}\begin{array}{rcl}\displaystyle \frac{{M}_{{\rm{d}},{\rm{f}}{\rm{i}}{\rm{t}}}}{{M}_{{\rm{N}}{\rm{S}}}^{{\rm{b}}}} & = & max\left[\left(0.218\displaystyle \frac{1-2{C}_{{\rm{N}}{\rm{S}}}}{{\eta }^{1/3}}\right.\right.\\ & & -\,\left.{\left.0.028{\widetilde{R}}_{{\rm{I}}{\rm{S}}{\rm{C}}{\rm{O}}}\displaystyle \frac{{C}_{{\rm{N}}{\rm{S}}}}{\eta }-0.122\right)}^{1.358},0\right].\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 8 }$$

This fitting formula applies to dynamical ejecta mass for binaries in the following parameter ranges Q ∈ [1, 7], χ_BH ∈ [0, 0.9],¹¹ and C_NS ∈ [0.108, 0.18]. The relative error of numerical data for 0.05 ≲ M_d,NR ≲ 0.1 M_⊙, M_d,NR ≃ 0.04 M_⊙, and M_d,NR ≃ 0.02 M_⊙ are fitted within ∼25%, ∼30%, and ∼40%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Additionally, we compare our dynamical ejecta fitting formula with the fitting model from Kawaguchi et al. (2016). The reduced χ² for our fit is

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{{N}_{{\rm{d}},{\rm{N}}{\rm{R}}}-{N}_{{\rm{p}}}-1}\sum _{i=1}^{{N}_{{\rm{d}},{\rm{N}}{\rm{R}}}}{\left(\displaystyle \frac{{M}_{{\rm{d}},{\rm{N}}{\rm{R}}}-{M}_{{\rm{d}},{\rm{f}}{\rm{i}}{\rm{t}}}}{{\sigma }_{{\rm{d}},{\rm{N}}{\rm{R}}}}\right)}^{2},\end{eqnarray} \tag{ 9 }$$

where N_d,NR = 64 and N_p = 4 are the numbers of the NR data points and the fitting parameters. The number of parameters changes to N_p = 6 when we use the fitting model of Kawaguchi et al. (2016). The χ² values of our result and those of Kawaguchi et al. (2016) are 2.926 and 3.443, respectively, suggesting comparable fitting results between the two fitting formulae.

The fitting formula of the total remnant mass (M_total,fit) and the dynamical ejecta mass (M_d,fit) are fitted with independent data. Within a certain range of parameters, we find that M_d,fit could be approximately, or even larger than, M_total,fit. However, the mass of the dynamical ejecta cannot exceed a few tens of percent of the total remnant mass outside the remnant BH. Therefore, we assume an upper limit on how much total remnant mass can become unbound dynamical ejecta, i.e.,

$$\begin{eqnarray}&&{M}_{{\rm{d}},max}={f}_{max}{M}_{{\rm{t}}{\rm{o}}{\rm{t}}{\rm{a}}{\rm{l}},{\rm{f}}{\rm{i}}{\rm{t}}},\end{eqnarray} \tag{ 10 }$$

where f_max is the maximum fractional factor. The largest unbound component based on the NR results we collected is ∼50% (the model MS1-Q7a5 in Kyutoku et al. 2015; i.e., No. 61 in Table 1). We therefore roughly define the factor f_max = 0.5. The mass of the dynamical ejecta can then be expressed as

$$\begin{eqnarray}&&{M}_{{\rm{d}}}=min({M}_{{\rm{d}},{\rm{f}}{\rm{i}}{\rm{t}}},{f}_{max}{M}_{{\rm{t}}{\rm{o}}{\rm{t}}{\rm{a}}{\rm{l}},{\rm{f}}{\rm{i}}{\rm{t}}}),\end{eqnarray} \tag{ 11 }$$

while the remnant disk mass around the BH can be estimated as

$$\begin{eqnarray}&&{M}_{{\rm{d}}{\rm{i}}{\rm{s}}{\rm{k}}}={M}_{{\rm{t}}{\rm{o}}{\rm{t}}{\rm{a}}{\rm{l}},{\rm{f}}{\rm{i}}{\rm{t}}}-{M}_{{\rm{d}}}.\end{eqnarray} \tag{ 12 }$$

We also give a fitting model for the dynamical ejecta average velocity with 60 valid NR data points. In the literature, the definitions of the average velocity are inconsistent. The Kyoto group (e.g., Kyutoku et al. 2015; Kawaguchi et al. 2016) defines the average velocity as the rms velocity, whereas the linear average velocity is adopted in Foucart et al. (2014, 2017, 2019) and Brege et al. (2018). We transform the linear average velocity into the rms velocity by simply assuming a ratio of ∼1.11 between the two (Foucart et al. 2017). In addition, as discussed in Foucart et al. (2017), there is a Δ = 7.6 MeV/nuc difference in the internal energy of the ejecta matter between their NR simulation results and the Kyoto group's simulation results due to different assumptions of the equations of state (EOSs). Foucart et al. (2017) also mentioned that the ejecta will have a kinetic energy ∼4.6 MeV/nuc lower than predicted due to the energy lost during r-process nucleosynthesis (Metzger et al. 2010). Following this discussion, we unify our collected data points, which are summarized in Table 1. As a simple function of the mass ratio, the fitting model (see Figure 3) for the rms velocity of the dynamical ejecta reads

$$\begin{eqnarray}&&{v}_{{\rm{r}}{\rm{m}}{\rm{s}},{\rm{d}}}=(-0.441{Q}^{-0.224}+0.549)c.\end{eqnarray} \tag{ 13 }$$

This fitting formula can achieve good fitting in a wide Q range (from 1 to 7) within ∼15% of the relative errors.

Zoom In Zoom Out Reset image size

For BH–NS mergers with NS tidal disruption, there are three possible ejecta components. First, a fraction of neutron-rich matter is tidally ejected. Lacking weak-interaction processes such as neutrino irradiation, such dynamical ejecta should be lanthanide-rich, which has a relatively low electron fraction Y_e ≲ 0.2 (e.g., Foucart 2012; Foucart et al. 2014; Kyutoku et al. 2018).

Second, an accretion disk is formed around the remnant BH. At early times after the disk forms, mass loss is driven by neutrino heating thanks to the high temperature in the disk. Due to neutrino illumination, such a neutrino-driven wind ejecta is less neutron rich, with a high electron fraction Y_e ≳ 0.4. Because the neutrino luminosity decreases rapidly with time (Fernández & Metzger 2013; Just et al. 2015), such a wind, directed in the polar direction, only lasts for a short period of time so that only a few percent of disk material is ejected in this form. Third, in the later stages, viscous heating and angular momentum transport play important roles in the mass loss of the disk. The electron fraction of this equatorial-dominated viscosity-driven wind ejecta is low, lying in the range of Y_e ∼ 0.15−0.25 (e.g., Fernández & Metzger 2013; Just et al. 2015; Siegel & Metzger 2017; Metzger 2019).

Tanaka et al. (2019) found that the average effective gray opacities of the mixture of r-process elements are κ ∼ 20–30 cm² g⁻¹ for Y_e ≤ 0.20, κ ∼ 3–5 cm² g⁻¹ for Y_e ≈ 0.25−0.35, and κ ∼ 1 cm² g⁻¹ for Y_e = 0.40 at the temperature T = 5−10 × 10³ K. At lower temperatures, these opacities can decrease steeply. In order to simplify the model, we set the gray opacity of dynamical ejecta, neutrino-driven ejecta, and viscosity-driven ejecta as κ_d = 20 cm² g⁻¹, κ_n = 1 cm² g⁻¹, and κ_v = 5 cm² g⁻¹, respectively.

Figure 4 summarizes the EM counterparts of a BH–NS merger and its mass ejection distribution. In this section, we model the dynamics and temperature evolution of each ejecta component. Because we assume that the three-component ejecta have a homologous expansion of the mass shells for each velocity, i.e., r = vt, one can substitute (r, t) with (v, t).

Zoom In Zoom Out Reset image size

3.1. Tidal Dynamical Ejecta

According to NR simulations of BH–NS mergers (e.g., Kyutoku et al. 2013, 2015; Kawaguchi et al. 2015), the mass distribution of dynamical ejecta is highly anisotropic, with the mass mainly distributed around the equatorial plane and shaped like a crescent. Kyutoku et al. (2013, 2015) and Kawaguchi et al. (2016) indicated that the distribution angle range of the ejecta in the longitudinal direction is almost φ_d ≈ π, while the half-opening angle in the latitudinal direction is typically in the range of θ_d ≈ 10°–20°. Hereafter, we assume that the half-opening angle in the latitudinal direction is a constant with θ_d = 15°. We assume that the dynamical ejecta of the BH–NS mergers have a relatively flat distribution velocity profile with v_min,d < v < v_max,d (Kyutoku et al. 2015; Kawaguchi et al. 2016), i.e., dM_d/dv ≈ const. The mass normalization is determined by

$$\begin{eqnarray}&&{\int }_{{v}_{min,{\rm{d}}}}^{{v}_{max,{\rm{d}}}}\displaystyle \frac{{dM}_{{\rm{d}}}}{dv}dv={M}_{{\rm{d}}}.\end{eqnarray} \tag{ 14 }$$

Similarly, the normalization can be also determined by the mass distribution of the dynamical ejecta in the latitudinal direction. Following Kawaguchi et al. (2016), we simply assume that the mass distribution of ejecta in the latitudinal direction is homogeneous¹² . The volume element of the dynamical ejecta is ${dV}={r}^{2}\cos \theta {drd}\theta d\varphi$ = ${v}^{2}{t}^{3}\cos \theta {dvd}\theta d\varphi$. Therefore, another normalization of the dynamical ejecta mass is

$$\begin{eqnarray}&&{\int }_{{v}_{min,{\rm{d}}}}^{{v}_{max,{\rm{d}}}}{\int }_{0}^{{\theta }_{{\rm{d}}}}{\int }_{-{\varphi }_{{\rm{d}}}/2}^{{\varphi }_{{\rm{d}}}/2}2{\rho }_{{\rm{d}}}(v,t){v}^{2}{t}^{3}\cos \theta dvd\theta d\varphi ={M}_{{\rm{d}}},\end{eqnarray} \tag{ 15 }$$

where φ_d and θ_d are the longitudinal half-opening angle and the latitudinal opening angle, respectively. Expanding the above integral, we get

$$\begin{eqnarray}&&2{\varphi }_{{\rm{d}}}{\theta }_{{\rm{d}}}{t}^{3}{\int }_{{v}_{min,{\rm{d}}}}^{{v}_{max,{\rm{d}}}}\rho (v,t){v}^{2}dv={M}_{{\rm{d}}}.\end{eqnarray} \tag{ 16 }$$

The dynamical ejecta in the latitudinal direction is approximately geometrically thin, i.e., sin θ_d ≈ θ_d ≪ 1. The ejecta profile is approximately

$$\begin{eqnarray}&&\displaystyle \frac{{dM}_{{\rm{d}}}}{dv}=2{\varphi }_{{\rm{d}}}{\theta }_{{\rm{d}}}{\rho }_{{\rm{d}}}(v,t){v}^{2}{t}^{3}.\end{eqnarray} \tag{ 17 }$$

Combining Equation (14) and Equation (17), one gets the density in the form of

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

In Section 2, we give a fitting model of the rms velocity of the dynamical ejecta. The kinetic energy of the dynamical ejecta is

$$\begin{eqnarray}&&{E}_{{\rm{K}},{\rm{d}}}=\displaystyle \frac{1}{2}{M}_{{\rm{d}}}{v}_{{\rm{r}}{\rm{m}}{\rm{s}},{\rm{d}}}^{2},\end{eqnarray} \tag{ 19 }$$

which can be also expressed as

$$\begin{eqnarray}&&{E}_{{\rm{K}},{\rm{d}}}={\int }_{{v}_{min,{\rm{d}}}}^{{v}_{max,{\rm{d}}}}\displaystyle \frac{1}{2}{v}^{2}\displaystyle \frac{{dM}_{{\rm{d}}}}{dv}dv.\end{eqnarray} \tag{ 20 }$$

Based on the NR results (e.g., Kyutoku et al. 2015; Brege et al. 2018), we find that the lower positions at the half-maximum of dM_d/dv are nearly at 0.1c so that we set the minimum velocity of dynamical ejecta to be v_min,d ≃ 0.1c.¹³ Combining Equations (19) and (20), one can express the maximum velocity as

$$\begin{eqnarray}&&{v}_{max,{\rm{d}}}=\sqrt{3{v}_{{\rm{r}}{\rm{m}}{\rm{s}},{\rm{d}}}^{2}-\displaystyle \frac{3}{4}{v}_{min,{\rm{d}}}^{2}}-\displaystyle \frac{1}{2}{v}_{min,{\rm{d}}}.\end{eqnarray} \tag{ 21 }$$

Figure 5 shows a schematic diagram of the dynamical ejecta. We define a critical surface where τ ≈ c/v as the photon diffusion surface (red dashed lines) below which photons are trapped in the ejecta because photon diffusion is slower than ejecta expansion. Above this critical surface, the photon diffusion velocity is larger than the ejecta expansion velocity so that they can escape (even though photon energy can be changed). We also define the photosphere (red solid lines) from where photons are last scattered. This is at an optical depth of τ ≈ 2/3.

Zoom In Zoom Out Reset image size

For a given point in the dynamical ejecta, photons diffuse in all directions. However, for the four directions denoted A, B, C, and D, the shortest diffusion times should be the latitudinal directions, radial direction, and the antiradial direction, respectively. Considering random walk of the photons, the diffusion time can be expressed as t_diff ∼ (R/λ)²(λ/c), where R is the distance between the position of the given point and the surface, and λ = (κ_dρ_d)⁻¹ is the mean free path of photons. Therefore, the diffusion times in the latitudinal, radial, and antiradial directions can be estimated as

$$\begin{eqnarray}&&{t}_{{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}},{\rm{l}}{\rm{a}}{\rm{t}}}(\theta ,t)\sim \displaystyle \frac{{\left({\theta }_{{\rm{d}}}-\theta \right)}^{2}{\kappa }_{{\rm{d}}}{M}_{{\rm{d}}}}{2c{\varphi }_{{\rm{d}}}{\theta }_{{\rm{d}}}({v}_{max,{\rm{d}}}-{v}_{min,{\rm{d}}})t},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{t}_{{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}},{\rm{r}}{\rm{a}}{\rm{d}}}(v,t)\sim \displaystyle \frac{{\left({v}_{max,{\rm{d}}}-v\right)}^{2}{\kappa }_{{\rm{d}}}{M}_{{\rm{d}}}}{2c{\varphi }_{{\rm{d}}}{\theta }_{{\rm{d}}}({v}_{max,{\rm{d}}}-{v}_{min,{\rm{d}}}){v}^{2}t},\end{eqnarray} \tag{ 23 }$$

and

$$\begin{eqnarray}&&{t}_{{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}},-{\rm{r}}{\rm{a}}{\rm{d}}}(v,t)\sim \displaystyle \frac{{\left(v-{v}_{min,{\rm{d}}}\right)}^{2}{\kappa }_{{\rm{d}}}{M}_{{\rm{d}}}}{2c{\varphi }_{{\rm{d}}}{\theta }_{{\rm{d}}}({v}_{max,{\rm{d}}}-{v}_{min,{\rm{d}}}){v}^{2}t}.\end{eqnarray} \tag{ 24 }$$

By comparing these three diffusion times, we can find the two boundaries shown by the light-gray dashed lines in Figure 5, i.e.,

$$\begin{eqnarray}\begin{array}{r}\begin{array}{r}\begin{array}{rcl}{v}_{{\rm{b}}{\rm{o}}{\rm{u}}1} & = & {v}_{min,{\rm{d}}}{\left(\theta -{\theta }_{{\rm{d}}}+1\right)}^{-1}\approx {v}_{min,{\rm{d}}}({\theta }_{{\rm{d}}}-\theta +1),\\ {v}_{{\rm{b}}{\rm{o}}{\rm{u}}2} & = & {v}_{max,{\rm{d}}}{\left({\theta }_{{\rm{d}}}-\theta +1\right)}^{-1}\approx {v}_{max,{\rm{d}}}(\theta -{\theta }_{{\rm{d}}}+1).\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 25 }$$

The latter approximations of the above two equations are due to the thin half-opening angle of the dynamical ejecta in the latitudinal directions. These two boundaries divide the dynamical ejecta into four different regions. For Regions A and B in Figure 5, the directions of the temperature gradient is along the latitudinal directions; for Region C, it is along the radial direction; and for Region D, it is along the antiradial direction. The photon diffusion directions are opposite to the directions of the temperature gradients.

One can also calculate the optical depths in each of the three directions. The optical depths in the latitudinal, radial, and antiradial directions are

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}{\tau }_{{\rm{d}},{\rm{l}}{\rm{a}}{\rm{t}}}(v,\theta ,t) & \approx & {\int }_{\theta }^{{\theta }_{{\rm{d}}}}{\kappa }_{{\rm{d}}}{\rho }_{{\rm{d}}}(v,t)vtd\theta \\ & = & \displaystyle \frac{{\kappa }_{{\rm{d}}}{M}_{{\rm{d}}}({\theta }_{{\rm{d}}}-\theta )}{2{\varphi }_{{\rm{d}}}{\theta }_{{\rm{d}}}({v}_{max,{\rm{d}}}-{v}_{min,{\rm{d}}}){vt}^{2}},\end{array}\end{array}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}{\tau }_{{\rm{d}},{\rm{r}}{\rm{a}}{\rm{d}}}(v,t) & \approx & {\int }_{v}^{{v}_{max,{\rm{d}}}}{\kappa }_{{\rm{d}}}{\rho }_{{\rm{d}}}(v,t)tdv\\ & = & \displaystyle \frac{{\kappa }_{{\rm{d}}}{M}_{{\rm{d}}}({v}^{-1}-{v}_{max,{\rm{d}}}^{-1})}{2{\varphi }_{{\rm{d}}}{\theta }_{{\rm{d}}}({v}_{max,{\rm{d}}}-{v}_{min,{\rm{d}}}){t}^{2}},\end{array}\end{array}\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}{\tau }_{{\rm{d}},-{\rm{r}}{\rm{a}}{\rm{d}}}(v,t) & \approx & {\int }_{{v}_{min,{\rm{d}}}}^{v}{\kappa }_{{\rm{d}}}{\rho }_{{\rm{d}}}(v,t)tdv\\ & = & \displaystyle \frac{{\kappa }_{{\rm{d}}}{M}_{{\rm{d}}}({v}_{min,{\rm{d}}}^{-1}-{v}^{-1})}{2{\varphi }_{{\rm{d}}}{\theta }_{{\rm{d}}}({v}_{max,{\rm{d}}}-{v}_{min,{\rm{d}}}){t}^{2}},\end{array}\end{array}\end{eqnarray} \tag{ 28 }$$

respectively. By comparing these optical depths, one can derive two boundaries, v_bou1 = v_min,d(θ − θ_d + 1) and v_bou2 = v_max,d(θ − θ_d + 1), which are equal to those derived in Equation (25). The fact that the boundaries defined by the photon diffusion surface and by the photosphere are consistent with each other makes sense, as the optical depths at the two surfaces are connected through the same mathematical conditions. It also suggests that the physical conditions in adjacent regions (e.g., A versus C, A versus D, B versus C, and B versus D) are continuous.

In the following, we will use an approximation method to model the temperature gradient of the dynamical ejecta for different regions. First, we consider the temperature gradient in the latitudinal direction. Due to the thin half-opening angle of the dynamical ejecta in the latitudinal direction, we assume that all photons escape the ejecta along the v_z-axis. The skin-depth angle θ_diff,d (below the ejecta surface) where the photons can diffuse vertically out of the ejecta surface within the dynamical time is given by (θ_d − θ_diff,d)vt ≈ (c/τ_d,lat)t and τ_d,lat ≈ κ_dρ_d (θ_d − θ_diff,d)vt, i.e.,

$$\begin{eqnarray}&&{\theta }_{{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}},{\rm{d}}}(t)={\theta }_{{\rm{d}}}\left(1-\displaystyle \frac{t}{{t}_{{\rm{c}}}}\right),\end{eqnarray} \tag{ 29 }$$

where we set the critical diffuse timescale that all ejecta can be seen for each velocity v as

$$\begin{eqnarray}&&{t}_{{\rm{c}}}=\sqrt{\displaystyle \frac{{\kappa }_{{\rm{d}}}{M}_{{\rm{d}}}{\theta }_{{\rm{d}}}}{2c{\varphi }_{{\rm{d}}}({v}_{max,{\rm{d}}}-{v}_{min,{\rm{d}}})}}.\end{eqnarray} \tag{ 30 }$$

Consider the luminosity dL_d from the part of ejecta in a velocity bin of v to v + dv. We assume that the specific energy injection rate due to radioactive decay is $\dot{\epsilon }{(t)={\epsilon }_{{Y}_{e}}{\dot{\epsilon }}_{0}(t/\mathrm{day})}^{-s}$ with ${\dot{\epsilon }}_{0}\approx 1.58\times {10}^{10}\,\mathrm{erg}\,{{\rm{g}}}^{-1}\,{{\rm{s}}}^{-1}$ and s ≈ 1.3, where ${\epsilon }_{{Y}_{e}}$ is an electron-fraction-dependent term which takes into account extremely neutron-rich ejecta with a decay half-life of a few hours (Perego et al. 2017), i.e.,

$$\begin{eqnarray}{\epsilon }_{{Y}_{e}}=\left\{\begin{array}{ll}0.5+2.5{\left[1+{e}^{4(t/\mathrm{day}-1)}\right]}^{-1}, & \mathrm{if}\ {Y}_{e}\geqslant 0.25,\\ 1, & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 31 }$$

At t < t_c, only the photons at a depth smaller than θ_d can escape within a dynamical time and hence contribute to the emission; at t ≥ t_c, all the photons can escape. Therefore, the luminosity per unit velocity is

$$\begin{eqnarray}\begin{array}{r}\displaystyle \frac{{dL}_{{\rm{d}}}}{dv}=\displaystyle \frac{1}{2}{\epsilon }_{{\rm{t}}{\rm{h}}}{\epsilon }_{{Y}_{e}}{\dot{\epsilon }}_{0}\displaystyle \frac{{dM}_{{\rm{d}}}}{dv}\left\{\begin{array}{ll}\displaystyle \frac{t}{{t}_{{\rm{c}}}}{\left(\displaystyle \frac{t}{{\rm{d}}{\rm{a}}{\rm{y}}}\right)}^{-s}, & t\lt {t}_{{\rm{c}}},\\ {\left(\displaystyle \frac{t}{{\rm{d}}{\rm{a}}{\rm{y}}}\right)}^{-s}, & t\geqslant {t}_{{\rm{c}}}\end{array}\right.,\end{array}\end{eqnarray} \tag{ 32 }$$

where the efficiency of thermalization is taken to be _th = 0.5 (Metzger et al. 2010) and the factor of 1/2 accounts for the two surfaces of the ejecta in the latitudinal directions.

The expression for dL_d/dv enables us to calculate the local emissivity per unit area, D_d(v, t), which is approximately a perfect blackbody. Because dL_d is released over an area of φ_drdr, we get

$$\begin{eqnarray}\begin{array}{r}{D}_{{\rm{d}}}(v,t)=\displaystyle \frac{{dL}_{{\rm{d}}}}{{\varphi }_{{\rm{d}}}rdr}\propto \left\{\begin{array}{ll}{v}^{-1}{t}^{-(s+1)}, & t\lt {t}_{{\rm{c}}},\\ {v}^{-1}{t}^{-(s+2)}, & t\geqslant {t}_{{\rm{c}}}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 33 }$$

Therefore, one can derive the effective temperature of each velocity at a given time,

$$\begin{eqnarray}&&{T}_{{\rm{e}}{\rm{f}}{\rm{f}},{\rm{d}}}(v,t)={\left(\displaystyle \frac{{D}_{{\rm{d}}}(v,t)}{{\sigma }_{{\rm{S}}{\rm{B}}}}\right)}^{1/4},\end{eqnarray} \tag{ 34 }$$

where σ_SB is the Steffan–Boltzmann constant. We then use the Eddington approximation (Mihalas 1970; Rybicki & Lightman 1979) to describe thermal temperature of each point in velocity space. The internal thermal temperature can be written as

$$\begin{eqnarray}&&{T}_{{\rm{d}}}(v,\theta ,t)={T}_{{\rm{e}}{\rm{f}}{\rm{f}},{\rm{d}}}(v,t){\left[\displaystyle \frac{3}{4}\left({\tau }_{{\rm{d}},{\rm{l}}{\rm{a}}{\rm{t}}}(v,\theta ,t)+\displaystyle \frac{2}{3}\right)\right]}^{1/4}.\end{eqnarray} \tag{ 35 }$$

We have considered the temperature gradient if all photons escape the ejecta along the v_z direction, which has equal gradient directions in Regions A and B. In the following, we model the direction of temperature gradient in Regions C and D. If all photons escape from the radial direction, as for the dynamical ejecta with a relatively flat velocity distribution profile in the radial direction, the temperature of each velocity v would be equal. Because the physical conditions at the boundaries between adjacent regions (A versus C, A versus D, B versus C, and B versus D) are continuous, one can simply map the velocity gradient (i.e., temperature gradient) in Regions C and D by that in Regions A and B through the adjacent boundaries. The temperature distributions of the dynamical ejecta in all four regions are qualitatively depicted with the color scheme in Figure 5.

3.2. Black Hole Disk Outflows

3.2.1. Neutrino-driven Wind Ejecta

Different from BNS mergers that may have significant neutrino-driven ejecta due to neutrino heating by a remnant supramassive or hypermassive NS that can release copious thermal neutrinos (e.g., Dessart et al. 2009; Perego et al. 2014; Martin et al. 2015), BH–NS mergers only have a limited amount of neutrino-driven ejecta. Just et al. (2015) noted that only ∼1% of the disk mass around the remnant BH can contribute to the neutrino-driven ejecta. In our calculation, we assume that the neutrino-driven ejecta is a constant fraction of the remnant disk mass, i.e., M_n = ξ_nM_disk with ξ_n ∼ 0.01. Martin et al. (2015) indicated that the neutrino-driven ejecta is a polar emission, which has a rather uniform mass distribution, i.e., F_n(θ) ≈ const, for θ ≤ θ_n, where θ_n is the half-opening angle of the neutrino-driven ejecta in the latitudinal direction. The normalization of M_n is

$$\begin{eqnarray}&&{M}_{{\rm{n}}}=2{\int }_{0}^{{\theta }_{{\rm{n}}}}\displaystyle \frac{{M}_{{\rm{n}}}}{2(1-\cos {\theta }_{{\rm{n}}})}\sin \theta d\theta .\end{eqnarray} \tag{ 36 }$$

We note that the latitudinal angle θ of the dynamical ejecta and the latitudinal angle θ for the neutrino-driven ejecta and viscosity-driven ejecta are defined differently. The former is defined as the angle with respect to the v_x–v_y plane, while the latter is defined as the angle from the ${{\boldsymbol{v}}}_{z}$ direction.

Wollaeger et al. (2018) indicated that for a spherically symmetric radiation-dominated outflow with adiabatic index Γ = 4/3, its expansion profile is ${dm}/{dv}\propto {\left(1-{\left(v/{v}_{\max }\right)}^{2}\right)}^{3}$. We then simply define the neutrino-driven ejecta mass distribution in the radial direction according to a similar profile, i.e., ${{dm}}_{{\rm{n}}}/{dv}\propto {\left(1-{\left(v/{v}_{\max ,{\rm{n}}}\right)}^{2}\right)}^{3}$, where dm_n is the neutrino-driven ejecta mass of the material with a certain θ and φ, and ${v}_{\max ,{\rm{n}}}$ is its maximum velocity. Also, dm_n/dv can be expressed as

$$\begin{eqnarray}&&2\pi {\int }_{0}^{{v}_{\max ,{\rm{n}}}}\displaystyle \frac{{{dm}}_{{\rm{n}}}}{{dv}}{dv}=\displaystyle \frac{{{dM}}_{{\rm{n}}}}{\sin \theta d\theta }.\end{eqnarray} \tag{ 37 }$$

Combining Equations (36) and (37), one gets the mass per velocity with a certain θ and φ:

$$\begin{eqnarray}&&\displaystyle \frac{{{dm}}_{{\rm{n}}}}{{dv}}=\displaystyle \frac{35{M}_{{\rm{n}}}}{64\pi (1-\cos {\theta }_{{\rm{n}}}){v}_{\max ,{\rm{n}}}}{\left[1-{\left(\displaystyle \frac{v}{{v}_{\max ,{\rm{n}}}}\right)}^{2}\right]}^{3}.\end{eqnarray} \tag{ 38 }$$

Another normalization of M_n is

$$\begin{eqnarray}&&2\times 2\pi {\int }_{0}^{{\theta }_{{\rm{n}}}}{\int }_{0}^{{v}_{\max ,{\rm{n}}}}{\rho }_{{\rm{n}}}(v,t){v}^{2}{t}^{3}\sin \theta {dvd}\theta ={M}_{{\rm{n}}},\end{eqnarray} \tag{ 39 }$$

where ρ_n is the density of the neutrino-driven ejecta and ${dV}={r}^{2}\sin \theta {drd}\theta d\varphi$ is the volume element. The density can be calculated by combining Equations (37), (38), and (39):

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{{\rm{n}}}(v,\theta ,t) & = & \displaystyle \frac{{{dm}}_{{\rm{n}}}}{{v}^{2}{t}^{3}{dv}}\\ & = & \displaystyle \frac{35{M}_{{\rm{n}}}}{64\pi (1-\cos {\theta }_{{\rm{n}}}){v}_{\max ,{\rm{n}}}{v}^{2}{t}^{3}}{\left[1-{\left(\displaystyle \frac{v}{{v}_{\max ,{\rm{n}}}}\right)}^{2}\right]}^{3}.\end{array}\end{eqnarray} \tag{ 40 }$$

It is easy to calculate the rms velocity:

$$\begin{eqnarray}&&{v}_{\mathrm{rms},{\rm{n}}}=\sqrt{\tfrac{2}{{M}_{{\rm{n}}}}\displaystyle \int \tfrac{1}{2}{\rho }_{{\rm{n}}}{v}^{2}{dV}}=\displaystyle \frac{1}{3}{v}_{\max ,{\rm{n}}}.\end{eqnarray} \tag{ 41 }$$

Hereafter we define ${v}_{\max ,{\rm{n}}}=3{v}_{\mathrm{rms},{\rm{n}}}$. The simulation results of Martin et al. (2015) showed the rms velocity of neutrino-driven ejecta to be mainly in the range of ${v}_{\mathrm{rms},{\rm{n}}}\simeq 0.055\mbox{--}0.075\,c$, while the best-fitting results for the lightcurve of AT 2017gfo from Perego et al. (2017) showed the rms velocity to be v_rms,n ≈ 0.0667c. Based on these results, we directly set v_rms,n = 0.0667c in our calculations.

3.2.2. Viscosity-driven Wind Ejecta

Simulation results indicated that the mass fraction of the viscosity-driven ejecta in the remnant disk mass is determined by the spin of the remnant BH. It can range from ∼5% for a low-spin BH to ∼30% for a high-spin BH with the dimensionless spin parameter χ_BH ≃ 0.95 (Fernández et al. 2015; Just et al. 2015; Fujibayashi et al. 2020). As for the BH–NS mergers, the spin of the remnant BH is mainly dependent on the spin of the initial BH and the mass ratio between the BH and the NS (Kyutoku et al. 2011; Zappa et al. 2019). An NS can easily plunge into an antialigned BH which can hardly form a remnant disk and tidal dynamical ejecta. Pannarale (2013) and Zappa et al. (2019) found that the spin distribution of the remnant BH is in the range of χ_BH ≃ 0.6–0.9 if the spin of the initial BH is aligned with (i.e., χ_BH ≥ 0) the orbital angular momentum.¹⁴ We therefore simply set the viscosity-driven ejecta mass as a constant fraction of the disk mass M_v = ξ_vM_disk, where ξ_v = 0.2 (Fernández et al. 2015; Just et al. 2015; Siegel & Metzger 2017; Fujibayashi et al. 2020). Perego et al. (2017) assumed an equatorial-dominated outflow of the viscosity-driven ejecta F_v(θ) = sin² θ based on simulations (Wu et al. 2016; Lippuner et al. 2017; Siegel & Metzger 2017). We still use this mass distribution to model the viscosity-driven ejecta. For the total mass, the normalization with θ is

$$\begin{eqnarray}&&{M}_{{\rm{v}}}\propto {\int }_{0}^{\pi /2}{F}_{{\rm{v}}}(\theta )\sin \theta d\theta .\end{eqnarray} \tag{ 42 }$$

Therefore, one has

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{{\rm{v}}}}{d\theta }=\displaystyle \frac{3}{4}{M}_{{\rm{v}}}{\sin }^{3}\theta .\end{eqnarray} \tag{ 43 }$$

Similar to the neutrino-driven ejecta, we define dm_v as the mass of the material with a certain θ and φ. The mass per velocity at certain θ and φ is ${{dm}}_{{\rm{v}}}/{dv}\propto {(1-{(v/{v}_{\max ,{\rm{v}}})}^{2})}^{3}$, where ${v}_{\max ,{\rm{v}}}=3{v}_{\mathrm{rms},{\rm{v}}}$. The simulation results from Just et al. (2015) showed that the rms velocity of the viscosity-driven ejecta ejected from the BH remnant disk mainly lies in the range of v_rms,v ∼ 0.03−0.04c. Hereafter, we directly assume v_rms,v = 0.03c. The normalization can also be expressed as

$$\begin{eqnarray}&&2\pi {\int }_{0}^{{v}_{\max ,{\rm{v}}}}\displaystyle \frac{{{dm}}_{{\rm{v}}}}{{dv}}{dv}=\displaystyle \frac{{{dM}}_{{\rm{v}}}}{\sin \theta d\theta }.\end{eqnarray} \tag{ 44 }$$

Therefore, we have

$$\begin{eqnarray}&&\displaystyle \frac{{{dm}}_{{\rm{v}}}}{{dv}}=\displaystyle \frac{105{M}_{{\rm{v}}}{\sin }^{2}\theta }{128\pi {v}_{\max ,{\rm{v}}}}{\left[1-{\left(\displaystyle \frac{v}{{v}_{\max ,{\rm{v}}}}\right)}^{2}\right]}^{3}.\end{eqnarray} \tag{ 45 }$$

The density can be easily calculated as

$$\begin{eqnarray}&&{\rho }_{{\rm{v}}}(v,\theta ,t)=\displaystyle \frac{{{dm}}_{{\rm{v}}}}{{v}^{2}{t}^{3}{dv}}=\displaystyle \frac{105{M}_{{\rm{v}}}{\sin }^{2}\theta }{128\pi {v}_{\max ,{\rm{v}}}{v}^{2}{t}^{3}}{\left[1-{\left(\displaystyle \frac{v}{{v}_{\max ,{\rm{v}}}}\right)}^{2}\right]}^{3}.\end{eqnarray} \tag{ 46 }$$

3.2.3. Temperature Evolution Model

In the following, we present a method similar to the dynamical ejecta model to delineate the temperature distribution of the two-component wind ejecta.

First, we assume that the radiation direction and temperature gradient are along the radial direction for the two-component wind ejecta. For the neutrino-driven ejecta, the optical depth along the radial direction is ${\tau }_{{\rm{n}}}(v)\approx {\int }_{v}^{{v}_{\max ,{\rm{n}}}}{\kappa }_{{\rm{n}}}{\rho }_{{\rm{n}}}{tdv}$. Photons will escape if $({v}_{\max ,{\rm{n}}}-{v}_{\mathrm{diff},{\rm{n}}})t\approx {ct}/{\tau }_{{\rm{n}}}({v}_{\mathrm{diff},{\rm{n}}})$, where ${\tau }_{{\rm{n}}}({v}_{\mathrm{diff},{\rm{n}}})\,\approx {\int }_{{v}_{\mathrm{diff},{\rm{n}}}}^{{v}_{\max ,{\rm{n}}}}{\kappa }_{{\rm{n}}}{\rho }_{{\rm{n}}}{tdv}$ and v_diff,n is the diffusion velocity for photons to escape. One can then obtain v_diff,n(t) at a given time. Therefore, the radial comoving luminosity is

$$\begin{eqnarray}&&{{dL}}_{{\rm{n}}}^{{\prime} }=\displaystyle \frac{{{dL}}_{{\rm{n}}}}{{{ \mathcal D }}^{2}}=\dot{\epsilon }(t){{dM}}_{{\rm{n}}},\end{eqnarray} \tag{ 47 }$$

where ${{dM}}_{{\rm{n}}}={\rho }_{{\rm{n}}}{v}^{2}{t}^{3}\sin \theta {dvd}\theta d\varphi$, ${ \mathcal D }=1/[{\rm{\Gamma }}(1-\beta )]$ is the Doppler factor, and ${\rm{\Gamma }}=1/\sqrt{1-{\beta }^{2}}$. Similar to dynamical ejecta, one can calculate the total emissivity per unit area D_n. Because dL_n is released over an area of ${v}_{\mathrm{phot},{\rm{n}}}^{2}{t}^{2}\sin \theta d\theta d\varphi$, where the photosphere position ${\tau }_{{\rm{n}}}({v}_{\mathrm{phot},{\rm{n}}})=2/3$, one gets

$$\begin{eqnarray}&&{D}_{{\rm{n}}}(t)=\displaystyle \frac{\dot{\epsilon }(t)t{\int }_{{v}_{\mathrm{diff},{\rm{n}}}}^{{v}_{\max ,{\rm{n}}}}{{ \mathcal D }}^{2}{\rho }_{{\rm{n}}}{v}^{2}{dv}}{{v}_{\mathrm{phot},{\rm{n}}}^{2}}.\end{eqnarray} \tag{ 48 }$$

The effective temperature of the neutrino-driven ejecta is

$$\begin{eqnarray}&&{T}_{\mathrm{eff},{\rm{n}}}(t)={\left(\displaystyle \frac{{D}_{{\rm{n}}}(t)}{{\sigma }_{\mathrm{SB}}}\right)}^{1/4}.\end{eqnarray} \tag{ 49 }$$

With the Eddington approximation (Mihalas 1970; Rybicki & Lightman 1979), the thermal temperature of each point in the velocity space of the neutrino-driven ejecta is given by

$$\begin{eqnarray}&&{T}_{{\rm{n}}}(v,t)={T}_{\mathrm{eff},{\rm{n}}}(t){\left[\displaystyle \frac{3}{4}\left({\tau }_{{\rm{n}}}(v,t)+\displaystyle \frac{2}{3}\right)\right]}^{1/4}.\end{eqnarray} \tag{ 50 }$$

The temperature gradient of the viscosity-driven ejecta is along the radial direction, similar to the neutrino-driven ejecta. The optical depth of the viscosity-driven ejecta is ${\tau }_{{\rm{v}}}(v)\,\approx {\int }_{v}^{{v}_{\max ,{\rm{v}}}}{\kappa }_{{\rm{v}}}{\rho }_{{\rm{v}}}{tdv}$. Photons will escape if $({v}_{\max ,{\rm{v}}}-{v}_{\mathrm{diff},{\rm{v}}})t\,\approx {ct}/{\tau }_{{\rm{v}}}({v}_{\mathrm{diff},{\rm{v}}})$, where ${\tau }_{{\rm{v}}}({v}_{\mathrm{diff},{\rm{v}}})\approx {\int }_{{v}_{\mathrm{diff},{\rm{v}}}}^{{v}_{\max ,{\rm{v}}}}{\kappa }_{{\rm{v}}}{\rho }_{{\rm{v}}}{tdv}$. One can obtain v_diff,v(θ, t) as a function of time. We define the total emissivity per unit area D_v as

$$\begin{eqnarray}&&{D}_{{\rm{v}}}(t)=\displaystyle \frac{\dot{\epsilon }(t)t{\int }_{{v}_{\mathrm{diff},{\rm{v}}}}^{{v}_{\max ,{\rm{v}}}}{{ \mathcal D }}^{2}{\rho }_{{\rm{v}}}{v}^{2}{dv}}{{v}_{\mathrm{phot},{\rm{v}}}^{2}},\end{eqnarray} \tag{ 51 }$$

where v_phot,v is for the photosphere position where τ_v(v_phot,v) = 2/3. The effective temperature of the viscosity-driven ejecta is

$$\begin{eqnarray}&&{T}_{\mathrm{eff},{\rm{v}}}(\theta ,t)={\left(\displaystyle \frac{{D}_{{\rm{v}}}(\theta ,t)}{{\sigma }_{\mathrm{SB}}}\right)}^{1/4}.\end{eqnarray} \tag{ 52 }$$

The thermal temperature of each point in the velocity space of the viscosity-driven ejecta can also be derived by

$$\begin{eqnarray}&&{T}_{{\rm{v}}}(v,\theta ,t)={T}_{\mathrm{eff},{\rm{v}}}(\theta ,t){\left[\displaystyle \frac{3}{4}\left({\tau }_{{\rm{v}}}(v,\theta ,t)+\displaystyle \frac{2}{3}\right)\right]}^{1/4}.\end{eqnarray} \tag{ 53 }$$

Based on the above assumptions, there is an overlapping region between the neutrino-driven ejecta and the viscosity-driven ejecta in the parameter regimes of v ≤ v_max,v and 0 ≤ θ ≤ θ_n or (π/2 − θ_n) ≤ θ ≤ π/2. For simplicity, we directly overlay the density of the two-component wind ejecta in this region, but consider that each point of this region has a single temperature. The viscosity-driven ejecta is equatorial-dominated so that the mass of the overlapping portion accounts for only a small part of the viscosity-driven ejecta. Therefore, this assumption has a slight effect on our final results. A schematic diagram of the density and temperature is present in Figure 6.

Zoom In Zoom Out Reset image size

So far, we have assumed that the temperature gradient direction is along the radial direction for both neutrino-driven ejecta and viscosity-driven ejecta. However, photons can also escape from the side of the neutrino-driven ejecta. As shown in Figure 6, we divide the wind ejecta into three regions. We define the temperature gradient of Region A and Region C to be along the radial direction, while the temperature gradient of Region B along the direction perpendicular to the edge of the neutrino-driven wind. The boundaries shown by the two light-gray dashed lines are, respectively,

$$\begin{eqnarray}\begin{array}{rcl}{v}_{\mathrm{bou}1} & \approx & {v}_{\max ,{\rm{n}}}(\theta -{\theta }_{{\rm{n}}}+1),\\ {v}_{\mathrm{bou}2} & = & {v}_{\max ,{\rm{v}}}\left(\displaystyle \frac{\theta }{{\theta }_{{\rm{n}}}}\right)+{v}_{\max ,{\rm{n}}}(1-{\theta }_{{\rm{n}}})\left(1-\displaystyle \frac{\theta }{{\theta }_{{\rm{n}}}}\right).\end{array}\end{eqnarray} \tag{ 54 }$$

The upper boundary can be derived by comparing the diffusion time similar to the dynamical ejecta. The lower boundary is our hypothetical boundary. Both the photon diffusion surface and the photosphere are nearly continuous at the boundaries.

Using a similar method to deal with the temperature evolution of the dynamical ejecta, we can calculate the thermal temperatures, i.e., T_n,bou1(v_bou1, θ, t) and T_n,bou2(v_bou2, θ, t), located at the two boundaries using Equation (50). The optical depth in the direction perpendicular to the edge of the wind is τ_n,side(v, θ, t) = κ_nρ_n(θ_n − θ)vt. Using the Eddington approximation, one can obtain the thermal temperature of Region B of the wind ejecta.

With the above preparation, in this section we investigate the evolution of the temperature profile of the ejecta and the observed photosphere emission. We will calculate the evolution of the emergent spectra and lightcurves as a function of the viewing angle using a spatial discretization model as described in detail in Appendix B. In order to discuss the viewing angle effect, we define a spherical coordinate system with the v_z-axis aligned with the jet axis and the v_x–v_z plane dividing the dynamical ejecta equally and symmetrically. As shown in Figure 7, the line of sight forms a latitudinal viewing angle θ_view with respect to the v_z-axis and a longitudinal angle φ_view with respect to the v_x–v_z plane. The ranges of these two viewing angles are θ_view ∈ [0, π/2] and φ_view ∈ [0, π], respectively.

Zoom In Zoom Out Reset image size

As discussed in Section 2, the remnant disk mass, the dynamical ejecta mass, and the rms velocity of the dynamical ejecta depend on the following parameters: the mass ratio (or the BH mass), the dimensionless spin of the BH, the NS (gravitational) mass, the NS baryonic mass, and the compactness of the NS. The NS baryonic and gravitational masses are related through the compactness parameter (e.g., Lattimer & Prakash 2001; Coughlin et al. 2017; Gao et al. 2020),

$$\begin{eqnarray}&&{M}_{\mathrm{NS}}^{{\rm{b}}}={M}_{\mathrm{NS}}\left(1+0.8858{C}_{\mathrm{NS}}^{1.2082}\right),\end{eqnarray} \tag{ 55 }$$

so that there are only four independent parameters. The masses of the neutrino-driven ejecta and the viscosity-driven ejecta can be described as constant fractions of the total mass of the remnant disk. Therefore, by setting four parameters, one can calculate all the input parameters (the neutrino-driven ejecta mass, the viscosity-driven ejecta mass, the dynamical ejecta mass, and the rms velocity of the dynamical ejecta) needed to calculate kilonova emission. In the following, we give two example cases, including a small mass remnant and a large mass remnant outside of the BH. The four relevant parameters in the two cases are summarized in Table 2. Hereafter, we mark Case I and Case II as corresponding to the small remnant mass case and the large remnant mass case, respectively. At the end of this section, we also extend our results to different mass-ratio regimes.

4.1. Temperature Profile Evolution

As discussed in Section 3, we assume a homologous expansion of the mass shells for each velocity and specific mass distribution, and model the dynamical evolution of each component ejecta. We also assume that the kilonovae of BH–NS mergers are powered only by the radioactive decay of r-process nuclei. In Figure 8, we show the results of the evolution of the ejecta temperature profiles, Case I and Case II. The evolution of the photosphere (red solid lines) and photon diffusion surface (red dashed lines) in the rest frame is shown in both cases.

Zoom In Zoom Out Reset image size

Consistent with intuition, Case II has a larger mass outside of the remnant BH and therefore has a higher temperature profile compared with Case I. In both cases, at the beginning of the merger, the ejecta cool rapidly due to the rapid expansion of the ejecta (the size of the ejecta increases by orders of magnitude in a short duration). At later times, the cooling slows down as the relative expansion is smaller (the size only increases linearly with time). This can be clearly seen in Figure 8.

For the two-component BH disk wind ejecta, the matter is mainly concentrated in the region where the velocity is smaller than the rms velocity. The photosphere and the photon diffusion surface evolve quickly at the beginning of the emission due to the density distribution and rapid density change caused by the initial expansion. Compared with Case II, the photosphere and photon diffusion surface of Case I can penetrate deeper into the ejecta at the same time after the merger. As these two surfaces pass through the low-density part of the ejecta, their evolution rates gradually slow down, and they are close to no evolution at ∼10 days. The photosphere cannot penetrate into the central matter of the BH disk wind ejecta, indicating that it cannot become completely optical thin. In addition to the photosphere, the photon diffusion surface cannot penetrate the entire ejecta even for the small mass remnant case. The emission of a large portion of matter below the photon diffusion surface cannot contribute to the luminosity of the kilonova.

Different from the BH disk wind ejecta, the density distribution of the dynamical ejecta is relatively homogeneous. As a result, the evolution rate of the photosphere and photon diffusion surface does not show the tendency of slowing down. The photon diffusion surface spends ∼4 days and ∼7 days passing through the entire dynamical ejecta for Cases I and II, respectively. After this, the energy deposited in the entire ejecta can contribute to the luminosity of the kilonova. On the other hand, the dynamical ejecta can hardly become completely optically thin within 10.5 days after the merger. This is especially true for Case II, in which the evolution rate of the photosphere is so slow that it stays at the front edge of the dynamical ejecta at early times after the merger. For Case I, the photosphere can quickly penetrate into the ejecta but the photosphere never completely penetrates the ejecta within 10.5 days. For both cases, the low-velocity region of the dynamical ejecta has higher density and temperature, indicating that this region contributes more to the observed luminosity.

Some caveats are worth mentioning. For simplicity, we model the temperature profile evolution with a constant gray opacity approximation. The opacities of a mixture of r-process elements are strongly dependent on temperature and wavelength (e.g., Kasen et al. 2013; Tanaka & Hotokezaka 2013; Tanaka et al. 2017). Tanaka et al. (2019) found that the gray opacities are nearly constant for temperatures T = 5–10 × 10³ K and decease steeply at lower temperatures. However, a steep decrease of the gray opacities significantly occurs at the temperature T ≲ 1500 K. The fitting temperature of the emergent spectrum, shown in the following sections (see Figure 10), cools to ∼1500 K at t ∼ 7–10 days after the merger. This means that the photosphere can penetrate deeper into the ejecta at late times, and the ejecta may become completely optically thin at late times. We also use a simple blackbody approximation for our model. For more sophisticated simulations, one should also consider temperature "floor" behavior. For example, Barnes & Kasen (2013) found that the the effective temperature remains unchanged as the photosphere recedes when lanthanide-rich ejecta cool to the first ionization temperature of the lanthanides (T_La ≈ 2500 K).

4.2. Photosphere Evolution in the Observer Frame

The calculations presented in Section 4.1 are presented in the rest frame of the source, also called the laboratory frame. The photosphere defined there does not reflect the photosphere of a certain observer. The photosphere in the frame of an observer is defined by the optical depth in a particular direction, so that it is dependent on viewing angle. One should also consider the light-propagation effect, such that the photosphere is defined at "retarded" times for the same observational time, i.e., emission comes from the so-called equal-arrival-time surface. These effects are fully considered in our calculations (see Appendix B).

Figure 9 shows the 2D sectional drawings of the photosphere evolution in the observer frame in velocity space for the two cases (Case I, left, and Case II, right) with different viewing angles. We take φ_view = 0 and adopt five values of the viewing angle θ_view, i.e., θ_view = −90°, −45°, 0°, 45°, and 90°.

Zoom In Zoom Out Reset image size

The photosphere in the observer frame also evolves rapidly in the beginning and slows down later, becoming very slow after ∼10 days of merger. For all directions, the photosphere as seen by the observer cannot pass through the central matter of the BH disk wind ejecta due to its high density. The matter behind the core is also blocked. Therefore, an observer cannot see emission from the ejecta during the entire evolution period of kilonova emission.

Within 14.5 days (which is the timescale of our calculation), the photosphere in the frame of any observer also cannot pass through the entire dynamical ejecta. However, compared with the photosphere in the source frame, the photosphere in the observer frame continuously evolves, with the evolution rate not slowing down with time. For long-enough time, one may finally see the entire dynamical ejecta matter. A special case is for θ_view = 90°. The wind ejecta will be hidden from view by the dynamical ejecta. The contribution of the wind ejecta emission can be essentially ignored after the neutrino-driven ejecta become optically thin. The observer can only see the emission from the dynamical ejecta.

Because the half-opening angle in the latitudinal direction of the dynamical ejecta is really thin, there is a significant change of the projected photosphere area for observers with different lines of sight. Obviously, the face-on projected photosphere area would be larger than those at large latitudinal viewing angles θ_view. The projected photosphere area variation with latitudinal viewing angle is roughly a factor of ∼(1.5−3). The change of projected photosphere area would affect the kilonova luminosity. Also, the projected photosphere shape is highly asymmetric. Both the change of the area and the change of the projected photosphere shape would have effects on the polarization of the kilonova (e.g., Shapiro & Sutherland 1982; Li & Shen 2019), which we do not study in this work.

Our calculation results for photosphere evolution in the observer frame are again based on the constant gray opacity approximation. As we mention before, the photosphere properties can also be affected by the steeply decreasing gray opacity when the temperature drops below ∼5000 K (Tanaka et al. 2019). The actual evolution of the observer-frame photosphere emission at late times could be more complicated.

4.3. Evolution of the Emergent Spectra

With the projected surface and thermal temperature of the photosphere in the observer's frame solved, the total observed flux density can be expressed as (see Appendix B for details)

$$\begin{eqnarray}&&{F}_{\lambda }(\lambda ,{t}_{\mathrm{obs}})\approx \displaystyle \frac{2}{{h}^{4}{c}^{3}{D}_{{\rm{L}}}^{2}}{\displaystyle \iint }_{S}\displaystyle \frac{{{ \mathcal D }}^{3}{\left({hc}/{ \mathcal D }\lambda \right)}^{5}}{\exp ({hc}/{ \mathcal D }{k}_{{\rm{B}}}\lambda {T}_{\mathrm{mesh}}^{{ij}})-1}d{\sigma }^{{\prime} {ij}},\end{eqnarray} \tag{ 56 }$$

where D_L is the luminosity distance, h is the Planck constant, k_B is the Boltzmann constant, λ is the wavelength, ${T}_{\mathrm{mesh}}^{{ij}}$ is the temperature at the mesh grid of our spatial discretization model, and $d{\sigma }^{{\prime} {ij}}$ is the infinitesimal projected photosphere area, which reads

$$\begin{eqnarray}&&d{\sigma }^{{\prime} {ij}}={\left(\displaystyle \frac{{t}_{\mathrm{obs}}}{1-{p}_{\mathrm{phot}}^{{ij}}/c}\right)}^{2}{{dv}}_{\mathrm{mesh},x}^{{\prime} {ij}}{{dv}}_{\mathrm{mesh},y}^{{\prime} {ij}}.\end{eqnarray} \tag{ 57 }$$

Here, t_obs is the observational time, ${p}_{{\rm{phot}}}^{{ij}}$ is the velocity space distance between the photosphere points and the mesh grid plane, and ${v}_{\mathrm{mesh},x}^{{ij}}$, ${v}_{\mathrm{mesh},y}^{{ij}}$ are the velocity components in the mesh grid. The factor of ${t}_{\mathrm{obs}}/(1-{p}_{\mathrm{phot}}^{{ij}}/c)$ is introduced to account for the light-propagation effect. We set the luminosity distance to be D_L = 10 pc hereafter, so that magnitude stands for the absolute magnitude.

We show examples of face-on (θ_view = 0°, φ_view = 0°) emergent spectra and different component contributions for both cases in Figure 10. One can see that throughout the evolution, the emission is mainly contributed by the radiation of the dynamical ejecta. For both cases, the peak flux density of the neutrino-driven ejecta is only about half of that of the dynamical ejecta, even though the neutrino-driven ejecta has a higher temperature. The emission from the viscosity-driven ejecta is overshone by the neutrino-driven ejecta at the early time of emission, so that one can only see a little radiation from the viscosity-driven ejecta for the face-on geometry. Because the neutrino-driven ejecta has a lower opacity and a smaller mass than other components, its radiation fades out rapidly with time. After ∼1.5 days post-merger, the wavelengths of the peak emission move from optical to infrared. As the neutrino-driven ejecta becomes optically thin, the contribution of the viscosity-driven ejecta increases significantly. At late times of emission, for Case I, the viscosity-driven ejecta contribute to relatively short wavelengths due to its low opacity compared with the dynamical ejecta. The dynamical ejecta, on the other hand, contribute more to the long-wavelength band. As for Case II, the late-time emission is always contributed by the dynamical ejecta. We fit the total emergent spectra with a single-blackbody model. For both cases, the spectra can be approximately fitted by the single-temperature blackbody model. Consistent with intuition, the temperatures of Case II are higher than those of Case I, because more mass is outside the BH remnant. Case II has a total ejecta mass that is only about four times of that of Case I. The blackbody-fit temperature is only higher by a factor of ∼(1.1−1.2).

Zoom In Zoom Out Reset image size

In Figures 11 and 12, we show the φ_view- and θ_view-dependent emergent spectra at the same observational epochs after the merger. We first discuss the φ_view-dependent emergent spectra. One can find that for increasing φ_view, the early-stage observed luminosity would decrease while the late-stage observed luminosity would increase. For θ_view = 30° and θ_view = 60°, the observer can simultaneously observe all three components. Changing φ_view has little effect on the observed wavelength at peak flux density. As the observed temperature is not significantly modified by the change of viewing angle, the difference in the peak flux density can approximately represent that in the observed luminosity. Changing φ_view results in a difference of observed luminosity by a factor of ∼(1.05−1.1) for both cases at early times (t ≲ 1.5 days) and a factor of ∼(0.75–0.85) at late times (t ≳ 7.5 days). For θ_view = 90°, the dynamical ejecta can outshine the other two components completely when φ_view ≈ 0°–60°. Therefore, under this condition, the variation of φ_view can have a relatively larger effect on the fitting temperature of the emergent spectrum and the observed luminosity. The variation can be particularly significant for Case I at late times because its late-stage luminosity is mainly contributed by the viscosity-driven ejecta. The observed luminosity of φ_view-dependent emergent spectra varies by a factor from ∼1.1 for Case I and ∼1.3 for Case II at early times (t ≲ 1.5 days) to ∼0.2 for Case I and ∼0.5 for Case II at late times (t ≳ 7.5 days).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We next discuss the θ_view dependence. Figure 12 shows the θ_view-dependent emergent spectra, where we set φ_view = 0°, which always has the largest observed luminosity for each θ_view at early stages. The observed luminosity is the largest at the same epoch if θ_view = 0°. This is because the projected photosphere area is the largest along the line of sight. The observed luminosity would decrease with increasing θ_view. At each epoch for Case I, the maximum observed luminosity is ∼1.5–7 times of the minimum observed luminosity. For Case II, the factor varies from ∼1.5 to 3 times the minimum observed luminosity.

One can conclude that the dynamical ejecta is the main contributor to the kilonova emission of BH–NS mergers. This is mainly because the projected photosphere area of the dynamical ejecta is far larger than those of the other two components along any line of sight. Another reason is that BH–NS mergers lack lanthanide-free ejecta. Due to the large projected photosphere area of the dynamical ejecta, there is only a small part of the observed photosphere that can be covered by the other two ejecta components if we only change φ_view. This is why changing φ_view has little effect on the variation of the observed spectra. The light-propagation effect would mainly affect the late-stage observed luminosity. For a large θ_view condition, the observer would see the photons emitted from the dynamical ejecta ∼(v_max,d/c)t_obs days earlier if the dynamical ejecta move away from the observer (i.e., φ_view ∼ 180°). The earlier high-temperature dynamical ejecta can enhance the observed luminosity. This is the reason why the late-stage observed luminosity would be enhanced with increasing φ_view. However, the light-propagation effect has little effect on the early-stage observed luminosity due to the rapid change of the projected area. In contrast, the projected photosphere area would significantly decrease with increasing θ_view. However, this decreasing trend cannot cause a significant change in luminosity. With the change in viewing angle, the variations of the observed luminosity are mainly derived from the change of the projected photosphere area of the dynamical ejecta, which causes a difference of only a factor of ∼(2−3) of the observed luminosity for different epochs.

4.4. Viewing-angle-dependent Lightcurve

With the time-dependent kilonova emergent spectra calculated in Section 4.3, one can calculate viewing-angle-dependent bolometric and color lightcurves. The total bolometric luminosity is given by ${L}_{\mathrm{bol}}(t)=4\pi {D}_{{\rm{L}}}^{2}{\int }_{0}^{\infty }{F}_{\lambda }d\lambda$, where the flux density at the photon wavelength λ is given by Equation (56). Based on the conversion between frequency and wavelength, i.e., F_νdν = F_λdλ, the flux density at photon frequency ν can be expressed as F_ν = λ²F_λ/c. One can also obtain the multifrequency monochromatic AB magnitude defined by ${M}_{\nu }=-2.5{\mathrm{log}}_{10}({F}_{\nu }/3631\,\mathrm{Jy})$. Because we set D_L = 10 pc, M_ν actually is the AB absolute magnitude. In Figure 13, we show the examples of face-on bolometric lightcurves and ugrizJHK-band lightcurves for both cases. We also show the contributions of different components to the face-on bolometric lightcurves. Similar to the discussion on the emergent spectra in Section 4.3, one can clearly see that most of the radiation energy is contributed by the dynamical ejecta. The face-on bolometric lightcurves can be approximately described by a power law with two breaks. In the broken-power-law description, the features of the face-on bolometric lightcurves are as follows. (i) Before t ∼ 1 day, the luminosity of a BH–NS merger kilonova is about a few times 10⁴¹ erg s⁻¹. During this time range, the radiation of the neutrino-driven ejecta is significant, even though only ∼1/5 of the radiation energy is contributed by the neutrino-driven ejecta at t ∼ 0.5 days. The bolometric lightcurves are well represented by a power-law decay with a power-law index of −0.37 for both cases. (ii) With the neutrino-driven ejecta becoming optically thin, almost all of the kilonova emission is contributed by the dynamical ejecta. Between t ≳ 1 day and t = t_c (t_c ∼ 4 days for Case I and t_c ∼ 7 days for Case II), the temporal index for both cases becomes ∼−0.26. (iii) There is another break at t = t_c for the dynamical ejecta emission, which is followed by a steeper decay with index ∼−0.93 for Case I and ∼−1 for Case II. The contribution from the viscosity-driven ejecta becomes progressively important at late times.

Zoom In Zoom Out Reset image size

The resulting lightcurves on a timescale of ∼0.5 days have a bolometric luminosity ∼1.7 × 10⁴¹ erg s⁻¹ for Case I and ∼3.7 × 10⁴¹ erg s⁻¹ for Case II. Case II also has a longer evolution time compared with Case I. One can see that kilonovae from BH–NS mergers are typically fainter than ∼−14.5 mag in the optical and ∼−15 mag in the infrared for Case I, and fainter than ∼−15 mag in the optical and ∼−16 mag in the infrared for Case II. As for multiband lightcurves, different bands peak at different times but for the same band ,there is no significant difference in the peak time for the two cases except the K band, which shows a coincidence of the peak time with t_c. Because t_c carries the information of ejecta mass, one may use the K-band peak time to estimate the mass in the ejecta. Compared with the peak AB absolute magnitudes of each filter for the two cases, the differences are in the range of 0.6–0.8 mag.

Figures 14 and 15 show the predicted φ_view- and θ_view-dependent lightcurves. Both φ_view and θ_view can change the peak time in each band, which is caused by the light-propagation effect. The variation of the peak time in each band can be approximately estimated by $\sim -({v}_{max,{\rm{d}}}/c)\sin {\theta }_{{\rm{v}}{\rm{i}}{\rm{e}}{\rm{w}}}\cos {\varphi }_{{\rm{v}}{\rm{i}}{\rm{e}}{\rm{w}}}{t}_{{\rm{p}}{\rm{e}}{\rm{a}}{\rm{k}}}$, where t_peak is roughly the peak timescale of each band. Therefore, the variation of the peak time depends on the relative motion direction of the dynamical ejecta: the peak time would increase if the dynamical ejecta moves away from the observer (i.e., φ_view = 90°–180°) and decrease if it moves toward the observer (i.e., φ_view = 0°–90°). As φ_view varies when we set θ_view = 30° and θ_view = 60°, the change in the shape and the magnitude of multiband lightcurves is tiny. For both cases, the differences in magnitude in each band between the maximum to the minimum are ∼0.1 mag and ∼0.2 mag at θ_view = 30° and θ_view = 60°, respectively. At θ_view = 90°, varying φ_view can give a relatively obvious effect on the shape of late-time multiband lightcurves. The differences in the maximum magnitude of each band can be larger, i.e., in the range of 0.4–0.6 mag. At φ_view ∼ 0°–45°, the decay index of each band can be steeper. This is because along the line of sight in this range, the dynamical ejecta would block the emission of the other two components. However, at late times, the viscosity-driven ejecta would mainly contribute to the near-infrared band while the dynamical ejecta contribute to longer wavelengths. The lack of contribution from the viscosity-driven ejecta causes this steeper decay of each multiband lightcurve.

Zoom In Zoom Out Reset image size

We show θ_view-dependent lightcurves in Figure 15, where we set φ_view = 0°. One can see that the multiband lightcurves at θ_view = 0° are almost the same as those at θ_view = 30°. The differences in the maximum magnitude of each filter are in the range of 0.5–0.6 mag. Altogether, the magnitude differences caused by the viewing angle changes are approximately ∼1 mag.

Zoom In Zoom Out Reset image size

Kyutoku et al. (2015) indicated that the dynamical ejecta is always smaller than 0.1 M_⊙. This means that BH–NS merger kilonovae are always less luminous than ∼4.5 × 10⁴¹ erg s⁻¹. Corresponding to AB absolute magnitudes, they are fainter than ∼−15 mag in optical and ∼−16 mag in infrared. Both φ_view and θ_view can change the peak time of each filter due to the light-propagation effect. They have little effect on the shape of the multiband lightcurves, expect for θ_view = 90° and φ_view ∼ 0°–45°, in which case the decay index of the late-time multiband lightcurves can be steeper. The observed luminosity differences caused by viewing angle changes are in the range of a factor of ∼2−3, which corresponds to approximately ∼1 mag.

4.5. Viewing Angle Dependence for Different Mass Ratios

In the above subsections, we discussed the viewing angle effect on the emergent spectra and lightcurves by setting two cases with certain parameters. The masses of the dynamical ejecta and the viscosity-driven ejecta for both cases have a similar order of magnitude. A BH–NS binary system with a high mass ratio, e.g., Q = 5, can produce relativistic tidal dynamical ejecta whose projected surface area is much larger than those of the other components. Under these premises, we can thus draw the conclusions that the dynamical ejecta is the main contributor to the kilonova emission of BH–NS mergers, and that the variation of the observed luminosity is only a factor of ∼(2−3) as the viewing angle varies. However, for the near-equal-mass BH–NS mergers, simulation results from Foucart et al. (2019) showed that only a small amount of matter can become unbound dynamical ejecta. Moreover, Kyutoku et al. (2015) indicated that the lower-energy material remaining outside the apparent horizon for a smaller value of Q tends to form slower dynamical ejecta. The typical velocity for a near-equal-mass BH–NS merger lies in the range of ∼(0.1−0.15)c (Foucart et al. 2019), which can reduce the area of the dynamical ejecta. Therefore, it is necessary to discuss the properties of kilonovae produced from the near-equal-mass BH–NS mergers and the applicability of our conclusions with different mass ratios Q.

In the left and middle panels of Figure 16, by using our new fitting formula presented in Section 2, we plot the dynamical ejecta mass and the mass ratio between the dynamical ejecta and the viscosity-driven ejecta associated with Q and χ_BH, where we set a certain NS mass (M_NS = 1.35 M_⊙) and a stiff compactness of the NS (C_NS = 0.130). Here, as the mass of the neutrino-driven ejecta is only ∼1/20 of the mass of the viscosity-driven ejecta, we mainly compare the difference between the dynamical ejecta and the viscosity-driven ejecta. The right panel in Figure 16 also describes the relation of Q and the face-on projected surface area ratio between the dynamical ejecta and the viscosity-driven ejecta. As shown in Figure 16, for a mass range with Q ≳ 3, the projected area of the dynamical ejecta would be significantly larger than that of the viscosity-driven ejecta. Also, the lanthanide-rich dynamical ejecta would occupy a considerable portion of the material outside the remnant BH. Therefore, the emission from the dynamical ejecta would mainly contribute to BH–NS merger kilonovae, while the observed luminosity only varies by a factor of ∼(2−3), due to the variation of the projected photosphere area of the dynamical ejecta with respect to the viewing angle.

Zoom In Zoom Out Reset image size

For the near-equal-mass regime, as shown in Figure 16, most of the materials outside the remnant would form a bound disk to produce the viscosity-driven ejecta. Besides, the velocities of the unbound dynamical ejecta decrease, which can significantly reduce the area ratio between the dynamical ejecta and the viscosity-driven ejecta. It can be predicted that the kilonovae from near-equal-mass BH–NS mergers are much dimmer due to the lack of fast-moving dynamical ejecta. The dynamical ejecta also cannot block the emission of the other two components at any viewing angle, so that observers would simultaneously observe the emission from all components along the line of sight. As a result, the viewing angle variation of the observed luminosity and the variation of the lightcurve peak time due to the photon-propagation effect are both minor.

4.6. Comparison With GW 170817/AT 2017gfo

Another EM counterpart for BH–NS mergers is an sGRB and its broadband afterglows. The sGRB afterglow lightcurve sensitively depends on the viewing. It is therefore interesting to simultaneously model the predicted lightcurves for both kilonova and sGRB afterglow from the BH–NS mergers. In this section, we discuss the viewing-angle-dependent lightcurves for kilonovae and sGRB afterglows, and how the different parameter values affect the detectability of the kilonova for the on-axis configuration.

In BH–NS mergers with tidal disruption, the remnant BH would accrete from the remnant disk and launch a relativistic jet via the Blandford–Znajek mechanism (Blandford & Znajek 1977). The kinetic energy of the jet (Barbieri et al. 2019) may be estimated as

$$\begin{eqnarray}&&{E}_{{\rm{K}},{\rm{j}}{\rm{e}}{\rm{t}}}=\epsilon (1-{\xi }_{{\rm{n}}}-{\xi }_{{\rm{v}}}){M}_{{\rm{d}}{\rm{i}}{\rm{s}}{\rm{k}}}{c}^{2}{{\rm{\Omega }}}_{{\rm{H}}}^{2}f({{\rm{\Omega }}}_{{\rm{H}}}),\end{eqnarray} \tag{ 58 }$$

where  = 0.015, Ω_H is the dimensionless angular frequency at the horizon, which is determined by the final spin of the BH,

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{H}}}=\displaystyle \frac{{\chi }_{\mathrm{BH},{\rm{f}}}}{2(1+\sqrt{1-{\chi }_{\mathrm{BH},{\rm{f}}}^{2}})},\end{eqnarray} \tag{ 59 }$$

and f(Ω_H) = 1 + 1.38 ${{\rm{\Omega }}}_{{\rm{H}}}^{2}$ − $9.2{{\rm{\Omega }}}_{{\rm{H}}}^{4}$ is a correction factor for high-spin values (Tchekhovskoy et al. 2010). We use Equation (11) from Pannarale (2013) to calculate the final spin of the BH.

We apply a power-law structured jet model (Rossi et al. 2002; Zhang & Mészáros 2002). The angular distributions of the kinetic energy and Lorentz factor Γ are adopted to be (Ghirlanda et al. 2019; Salafia et al. 2019)

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{dE}}{d{\rm{\Omega }}}(\theta ) & = & \displaystyle \frac{{E}_{{\rm{c}}}/4\pi }{1+{\left(\theta /{\theta }_{{\rm{c}}}\right)}^{{s}_{1}}},\\ {\rm{\Gamma }}(0,\theta ) & = & 1+\displaystyle \frac{{{\rm{\Gamma }}}_{c}-1}{1+{\left(\theta /{\theta }_{{\rm{c}}}\right)}^{{s}_{2}}},\end{array}\end{eqnarray} \tag{ 60 }$$

where we set ${{\rm{\Gamma }}}_{{\rm{c}}}=250$, θ_c = 5°, s₁ = 5.5 and s₂ = 3.5, and ${E}_{{\rm{c}}}={E}_{{\rm{K}},\mathrm{jet}}/\pi {\theta }_{{\rm{c}}}^{2}$. These parameters have been used to model GW 170187/GRB 170817A. The standard GRB afterglow model is briefly introduced in Appendix C.

5.1. Viewing Angle Dependence

In our calculations to discuss the viewing-angle-dependent lightcurves of sGRB afterglows, the following typical afterglow model parameters are adopted: the fraction of shock energy carried by the magnetic field _B = 0.001, the fraction of shock energy carried by electrons _e = 0.1, and the ISM number density n = 5 × 10⁻³ cm⁻³.

In Figure 17, we present the lightcurves in three representative filters, i.e., the r, J, and K bands. Both the kilonova lightcurves and afterglow lightcurves are shown for Case I (left panel) and Case II (right panel). We adopt several possible viewing angles from θ_view = 0° to 30°. One can see that the change in the kilonova lightcurves is essentially negligible.

Zoom In Zoom Out Reset image size

For an on-axis view, i.e., θ_view = 0°–5° (θ_view = 0 − θ_c), the r-band flux density from the kilonovae is much less than the r-band flux density from the afterglow. For the J band, the kilonovae are less luminous than the afterglow most of the time, but may show up in ∼(2–7) days when the J-band flux densities between the two are comparable. The K-band emission after ∼2 days becomes dominated by the kilonova emission. Therefore, along the on-axis line of sight, the best filter to observe kilonovae from BH–NS mergers is the K band, even though the emission shows up at a relatively late epoch.

For an off-axis view, the early-stage kilonova emission is always much more luminous than the sGRB afterglow emission. The effect becomes more prominent as the viewing angle increases. However, as the sGRB blastwave decelerates (i.e., the Lorentz factor Γ of the external shock decreases), the more energetic jet core becomes visible, so that the afterglow lightcurves continuously rise and eventually outshine the kilonova emission in the r and J bands. The K band is likely dominated by the kilonova emission for an extended period of time.

5.2. Parameter Dependence

In order to discuss the parameter space where the kilonova can be observed in the on-axis view, we only show the comparison between the kilonova lightcurves and the afterglow lightcurves for Case I, because similar results can be obtained for Case II.

The left and middle panels of Figure 18, respectively, show the effects of the ISM number density n and the fraction of the shock energy carried by the magnetic field _B in the afterglow lightcurve. Here, we set the fraction of shock energy carried by electrons to be _e = 0.1, as the variation of this parameter constrained by observations has a relatively small change (e.g., Santana et al. 2014). The right panels of Figure 18 compare the difference in the flux density between the afterglow component and the kilonova component at the peak time of the kilonova lightcurve for each band. One can still conclude that the best filter to observe the kilonova from a BH–NS merger is the K band. Compared with the J band and K band in which the kilonova emission is observable for a large parameter space, the optical emission of the kilonova can only be observed in a very low-density environment with low _B.

Zoom In Zoom Out Reset image size

Observationally, Fong et al. (2015; see also O'Connor et al. 2020) found that most sGRBs occur in low-density environments with median density 3–15 × 10⁻³ cm⁻³. The systematic studies on magnetic fields in external forward shocks (Santana et al. 2014; Wang et al. 2015) showed that _B has a wide range of ∼10⁻⁸–10⁻³ (10⁻⁶–10⁻¹) and is mainly centered at ∼few × 10⁻⁵ (∼few × 10⁻³) by taking n = 1 cm⁻³ (n = 10⁻² cm⁻³). Therefore, as shown in the right panels of Figure 18, the optical emission of kilonovae is often outshone by the afterglow and not detectable, while the emission in the J band and K band can be easily detected. Only a small amount of BH–NS merger kilonovae can have their optical emission observable along the jet axis.

By considering three radioactivity-powered components, i.e., lanthanide-rich tidal dynamical ejecta, intermediate-opacity viscosity-driven wind ejecta, and lanthanide-free neutrino-driven wind ejecta, we modeled the dynamics and temperature profile evolution of the BH–NS merger kilonovae in great detail. Because numerical simulations show that these three components are highly anisotropic, we have paid special attention to the viewing angle effect on BH–NS merger kilonovae. We presented a numerical method to model the evolution of the photosphere in the observer frame and studied the emergent spectra and lightcurves for different lines of sight in two degrees of freedom for the viewing angles, i.e., the latitudinal viewing angle θ_view and the longitudinal viewing angle φ_view. The ejecta models presented here are simplified with the assumption of homologous expansion, constant gray opacity, and simple radiative transfer. More realistic models should consider complex dynamical evolution, thermodynamical evolution, and temperature- and wavelength-dependent opacities. Nevertheless, this simplified model provides valuable information for us to comprehend the characteristics of BH–NS merger kilonovae and the viewing angle effect on multiband kilonova lightcurves.

We find that the dynamical ejecta would contribute to the majority of kilonova emission from BH–NS mergers due to its projected photosphere area being the largest. Because NR simulations of BH–NS mergers revealed that the mass of the dynamical ejecta is always ≲0.1 M_⊙ (Kyutoku et al. 2015),¹⁶ the peak luminosity of BH–NS merger kilonovae would always be less luminous than ≲4.5 × 10⁴¹ erg s⁻¹. Corresponding to the AB absolute magnitudes commonly used by observers, its maximum absolute magnitude is ∼−15 mag in the optical and ∼−16 mag in the infrared.

Long-lived energy injection from the remnant BH may produce an additional source of ejecta heating in excess of the contribution from r-process radioactivity. A fraction of gravitationally bound material would fall back onto the BH and enter the disk over a range of timescales from seconds to days (Rosswog 2007; Kyutoku et al. 2015). The energy release from fallback accretion may enhance the peak brightness of the kilonova. For example, Ma et al. (2018) proposed an energy injection mechanism invoking a wind driven by the Blandford–Payne mechanism from an accretion disk (Blandford & Payne 1982). According to our prediction, a radioactivity-powered BH–NS merger kilonova would be always fainter than the critical luminosity, i.e., ∼4.5 × 10⁴¹ erg s⁻¹, at the peak time. Future observations of a BH–NS merger kilonova brighter than this critical luminosity would suggest additional energy injection from the central engine.

We compare our theoretical results on BH–NS merger kilonovae with the observational properties of AT 2017gfo. At each epoch after the merger, we find that the blackbody-fit temperature of BH–NS mergers kilonovae is lower than that of AT 2017gfo. Due to the lack of abundant lanthanide-free ejecta as in AT 2017gfo, BH–NS merger kilonovae are optically dim, but possibly infrared bright.

We showed that the observed luminosity of BH–NS merger kilonovae varies with the projected photosphere area determined by the viewing angles. The variation of the longitudinal viewing angle φ_view has little effect on the variation of the observed luminosity, while the variation of the latitudinal viewing angle θ_view can significantly change the projected photosphere area, and hence, affect the observed luminosity. In total, the difference of the observed luminosity caused by the variation in viewing angle is only a factor of ∼(2−3), corresponding to the change of the multiband magnitude by ∼1 mag. This is similar to the result of Roberts et al. (2011), who adopted three-dimensional radiation simulations to study the viewing angle effect on the emission of tidal tails. Furthermore, the blackbody-fitting temperature and the shape of the observed multiband lightcurves are not significantly dependent on the line of sight. However, for the case where the dynamical ejecta blocks most of the emission from the disk wind outflows along the line of sight, the decay index of multiband lightcurves would become steeper at late times. In addition, both φ_view and θ_view can affect the peak time of the multiband lightcurves. The variation of the peak time, caused by the light-propagation effect, depends on the relative motion direction of the dynamical ejecta. More specifically, the peak time would increase if the dynamical ejecta moves away from the observer and decrease if it moves toward the observer. Recently, Darbha & Kasen (2020) provided a simple analytic estimate of the viewing-angle-dependent lightcurves as a function of the projected surface area along the line of sight by considering three specific geometries, i.e., an ellipsoid, a ring torus, and a conical section embedded in a sphere. Their conclusions are qualitatively similar to ours. However, they derived that the viewing angle effect can cause a factor of ∼5−10 difference in the BH–NS merger kilonova luminosity, which is larger than ours. This may be introduced by their simplified treatment of the geometry. In particular, they modeled the BH–NS merger tidal tail as an oblate ellipsoid with an axial ratio R = 5 or a torus with a radius ratio K = 5. However, according to simulation results (e.g., Roberts et al. 2011; Kyutoku et al. 2015; Brege et al. 2018), the dynamical ejecta is concentrated near the equator with the opening angle in the longitudinal direction filling an arc of about π, which is shaped like a crescent. Therefore, one may simply assume that the dynamical ejecta is like a moving ellipsoid with a relatively broad axial ratio R ∼ 2−3. This would reduce the peak luminosity variation due to viewing angle variation to be consistent with our results.

The sGRB afterglows are very sensitive to the viewing angle and can significantly affect the detectability of BH–NS merger kilonovae. For an on-axis observer, an optical filter is not recommended to be used to observe the emission from the kilonova because it is completely outshone by the afterglow emission. The optical emission can only be observed if the sGRB occurs in a very low-density environment with a low _B. Redder filters are preferred to detect kilonova emission. In the J band, kilonova emission may be barely detected several days after the merger. The K band is the most ideal band as the kilonova emission becomes dominant after ∼2 days and lasts for an extended period of time.

For an off-axis geometry, the early-stage kilonova emission is always much more luminous compared with the afterglow emission. The effect is more prominent for larger viewing angles. In relatively blue bands, the afterglow emission will outshine the kilonova emission at late times as the bright jet core becomes visible.

Fujibayashi et al. (2020) recently found that the viscosity-driven ejecta, whose electron fraction mainly lies in the range of Y_e ∼ 0.3–0.4, if the viscous coefficient is not extremely high, can be lanthanide-free. A substantial mass of the blue-component ejecta may be formed after BH–NS merger, which may be similar to or even more than the mass of the blue component invoked to explain AT 2017gfo. As a result, an AT 2017gfo–like kilonova may be possible to generate for a BH–NS merger as well. We may simply estimate the viewing-angle-dependent observed luminosity using the fitting results of AT 2017gfo. For the face-on configuration, the peak observed luminosity may reach ∼1 × 10⁴² erg s⁻¹. The observed luminosity drops in viewing angles where the dynamical ejecta blocks the emission from the blue-component ejecta, with the minimum at e.g., θ_view ∼ 90° and φ_view ∼ 0°–45°. Because the mass of the dynamical ejecta mainly lies in range of (0.01–0.05) M_⊙ (e.g., Cowperthwaite et al. 2017; Murguia-Berthier et al. 2017; Kasen et al. 2017; Kasliwal et al. 2017; Perego et al. 2017; Tanaka et al. 2017; Villar et al. 2017), the peak observed luminosity would be in the range of ∼(1−2) × 10⁴¹ erg s⁻¹ for the case when the blue component is blocked. Overall, if BH–NS mergers indeed have a blue component, the peak observed luminosity can vary by a factor of ∼5–10 for different viewing angles.

We point out that our results on the viewing angle effect only apply to the mass ratio range of Q ≳ 3. The population-synthesis simulation results of the final mass distribution of BH–NS mergers (Giacobbo & Mapelli 2018; Mapelli & Giacobbo 2018) showed that the mass ratio of almost all merging systems is larger than Q > 2.5, and the most likely value is Q = 5. We can conclude that our results are always relevant for BH–NS merger kilonovae. Moreover, the existence of low-mass BHs cannot be theoretically ruled out. We also discuss the properties of BH–NS merger kilonovae in the near-equal-mass regime. We extend the fitting formulae for the mass and velocity of the dynamical ejecta across a wider range of mass ratios from Q = 1 to Q = 7 validated with 66 simulations. For near-equal-mass regime, our results show that the dynamical ejecta is slow for ejecta produced during the tidal disruption of a neutron star (see also Foucart et al. 2019). The projected photosphere of the dynamical ejecta would decrease remarkably so that the viewing angle effects on the peak luminosity of the lightcurves would change. Besides, a less massive dynamical ejecta can be ejected after the merger. One can predict that kilonovae of near-equal-mass BH–NS mergers are much dimmer than those of large mass ratios. More detailed studies on the parameter dependence and viewing angle dependence of the properties of BH–NS merger kilonovae, detection rate, and polarization are subject to further studies in the subsequent articles in this series. Apparently, the characteristics of BH–NS merger kilonovae are complex. As more BH–NS mergers are detected by LIGO/Virgo, one expects that their associated kilonovae will be eventually detected (Bhattacharya et al. 2019); see Zappa et al. (2019) and Tsujimoto et al. (2020). Future GW-led multimessenger observations of BH–NS mergers can help us to better understand the characteristics of these kilonovae.

We thank an anonymous referee for constructive suggestions, Koutarou Kyutoku and Masaomi Tanaka for valuable discussion and comments, and Kyohei Kawaguchi for providing helpful data and information. The work of J.P.Z. is partially supported by the National Science Foundation of China under grant No. 11721303 and the National Basic Research Program of China under grant No. 2014CB845800. L.D.L. is supported by the National Postdoctoral Program for Innovative Talents (grant No. BX20190044), China Postdoctoral Science Foundation (grant No. 2019M660515), and the "Li Yun" postdoctoral fellow of Beijing Normal University. Z.L. is supported by the National Natural Science Foundation of China under grant Nos. 11773003 and U1931201. Y.W.Y. is supported by the National Natural Science Foundation of China under grant Nos. 11822302 and 1183303. H.G. is supported by the National Natural Science Foundation of China under grant Nos. 11722324, 11690024, 11603003, and 11633001; the Strategic Priority Research Program of the Chinese Academy of Sciences, grant No. XDB23040100; and the Fundamental Research Funds for the Central Universities.

Software: Matlab, https://www.mathworks.com.

In Table A1, we summarize the definitions of frequently used variables.

Table A1.  Definitions of Frequently Used Variables

Case	Variable	Definition
BH–NS merger	MBH	Mass of the BH
	MNS	Mass of the NS
	MNSb	Baryon mass of the NS
	χBH	Dimensionless spin of the BH
	CNS	Compactness of the NS
	Q	Mass ratio between the BH mass and NS mass
Ejecta	Ye	Electron fraction
	${M}_{\mathrm{disk}}$	Mass of the remnant disk
	${M}_{{\rm{d}}}$	Mass of the tidal dynamical ejecta
	${\kappa }_{{\rm{d}}}$	Opacity of the dynamical ejecta; the value is ${\kappa }_{{\rm{d}}}=20\,{{\rm{c}}{\rm{m}}}^{2}\,{{\rm{g}}}^{-1}$
	${\theta }_{{\rm{d}}}$	Half-opening angle in the latitudinal direction of the dynamical ejecta; the value we set is ${\theta }_{{\rm{d}}}\approx {15}^{\circ }$
	${\varphi }_{{\rm{d}}}$	Opening angle in the longitudinal direction of the dynamical ejecta; the value is ${\varphi }_{{\rm{d}}}\approx \pi$
	${v}_{\min ,{\rm{d}}}$	Minimum velocity of the dynamical ejecta; the value we set is ${v}_{min,{\rm{d}}}=0.1\,c$
	${v}_{\mathrm{rms},{\rm{d}}}$	rms velocity of the dynamical ejecta
	${M}_{{\rm{n}}}$	Mass of the neutrino-driven wind ejecta
	${\kappa }_{{\rm{n}}}$	Opacity of the neutrino-driven ejecta; the value we set is ${\kappa }_{{\rm{n}}}=1\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$
	${\theta }_{{\rm{n}}}$	Opening angle in the latitudinal direction of the neutrino-driven ejecta; the value is ${\theta }_{{\rm{n}}}\approx 30^\circ$
	${v}_{\max ,{\rm{n}}}$	Maximum velocity of the neutrino-driven ejecta; the value we set is ${v}_{\max ,{\rm{n}}}=0.2\,c$
	${v}_{\mathrm{rms},{\rm{n}}}$	rms velocity of the neutrino-driven ejecta
	${M}_{{\rm{v}}}$	Mass of the viscosity-driven wind ejecta
	${\kappa }_{{\rm{v}}}$	Opacity of the viscosity-driven ejecta; the value we set is ${\kappa }_{{\rm{v}}}=5\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$
	${v}_{\max ,{\rm{v}}}$	Maximum velocity of the viscosity-driven ejecta; the value we set is ${v}_{\max ,{\rm{v}}}=0.09\,c$
	${v}_{\mathrm{rms},{\rm{v}}}$	rms velocity of the viscosity-driven ejecta
Viewing angle	${\theta }_{\mathrm{view}}$	Latitudinal viewing angle
	${\varphi }_{\mathrm{view}}$	Longitudinal viewing angle
Numerical method	${v}_{\mathrm{mesh},{\rm{x}}}^{{\prime} {ij}}$	${v}_{x}^{{\prime} }$ component of mesh grid points
	${v}_{\mathrm{mesh},{\rm{y}}}^{{\prime} {ij}}$	${v}_{y}^{{\prime} }$ component of mesh grid points
	${p}_{\mathrm{phot}}^{{ij}}$	Velocity space distance between the photosphere and the mesh grid plane
	${T}_{\mathrm{mesh}}^{{ij}}$	Thermal temperature of mesh grid points
Observational parameter	${{\rm{D}}}_{{\rm{L}}}$	Luminosity distance; the value we set is ${D}_{{\rm{L}}}=10\,\mathrm{pc}$
	Fλ	Flux density at photon wavelength λ
	Fν	Flux density at photon frequency ν
	Mν	AB absolute magnitude
	${t}_{\mathrm{obs}}$	Observational time
Afterglow	${\nu }_{{\rm{m}}}^{{\prime} }$	Synchrotron frequency of the accelerated electrons with the minimum Lorentz factor
	${\nu }_{{\rm{c}}}^{{\prime} }$	Cooling frequency
	${\nu }_{{\rm{a}}}^{{\prime} }$	Synchrotron self-absorption frequency
	n	Interstellar medium number density
	${\epsilon }_{B}$	Fraction of shock energy carried by the magnetic field
	${\epsilon }_{e}$	Fraction of shock energy carried by electrons

Download table as:  ASCIITypeset images: 1 2

In Section 2, we have established the dynamics and temperature evolution of different components in velocity space. In order to reconstruct viewing-angle-dependent lightcurves, one needs to solve the observed photosphere along the line of sight of an observer and integrate the projected photosphere. Rosswog et al. (2014), Grossman et al. (2014), Martin et al. (2015), Perego et al. (2017), and Barbieri et al. (2019, 2020) used a semianalytical method that divided ejecta into finite slices and summed up individual contributions of each slice to compute lightcurves with a one-dimensional viewing angle (see e.g., Kawaguchi et al. 2018, 2020b; Darbha & Kasen 2020; Korobkin et al. 2020 used radiative transfer simulations to calculate the viewing-angle-dependent lightcurve). This appendix is devoted to a numerical method of calculating the photosphere evolution and lightcurve as seen by observers with different viewing angles. We define the two-dimensional viewing angles (θ_view, φ_view) within the ranges of ${\theta }_{\mathrm{view}}\in [0,\pi /2]$ and φ_view ∈ [0, π], which we have shown in Figure 7.

The steps of establishing the evolution of the photosphere in the observer frame include the following, in brief:

• 1.  
  Comparing with the local coordinate system of all the ejecta components whose origin is O, we set up a three-dimensional Cartesian coordinate system with origin O' and the axes v'_x, ${v}_{y}^{{\prime} }$, and ${v}_{z}^{{\prime} }$. We get a mesh grid on the plane of ${v}_{x}^{{\prime} }{O}^{{\prime} }{v}_{y}^{{\prime} }$ and set ${v}_{z}^{{\prime} }{O}^{{\prime} }$ as the line of sight. The points of the mesh grid are marked as $({v}_{\mathrm{mesh},x}^{{\prime} {ij}},{v}_{\mathrm{mesh},y}^{{\prime} {ij}},0)$ where the superscripts (i and j) represent the IDs of the points at the mesh grid.
• 2.  
  Rotate the density profile of all the ejecta components from the coordinate system O to the coordinate system O' simultaneously. In order to rotate velocity points located at the profile of the ejecta, we follow Goldstein et al. (2002) and define three Euler rotation matrices which are rotations about the v_x-, v_y-, and v_z-axes using the right-hand rule:

  $$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{R}}}_{x}(\theta ) & = & \left(\begin{array}{ccc}1 & 0 & 0\\ 0 & \cos \theta & -\sin \theta \\ 0 & \sin \theta & \cos \theta \end{array}\right),\\ {{\boldsymbol{R}}}_{y}(\theta ) & = & \left(\begin{array}{ccc}\cos \theta & 0 & -\sin \theta \\ 0 & 1 & 0\\ \sin \theta & 0 & \cos \theta \end{array}\right),\\ {{\boldsymbol{R}}}_{z}(\theta ) & = & \left(\begin{array}{ccc}\cos \theta & \sin \theta & 0\\ -\sin \theta & \cos \theta & 0\\ 0 & 0 & 1\end{array}\right).\end{array}\end{eqnarray} \tag{ B1 }$$

  Take one point ${\boldsymbol{v}}={\left[{v}_{x},{v}_{y},{v}_{y}\right]}^{{\rm{T}}}$ located at the profile of the ejecta; for example, the point's location at the coordinate system O' after rotating is

  $$\begin{eqnarray}&&{{\boldsymbol{v}}}^{{\prime} }={{\boldsymbol{R}}}_{y}({\theta }_{\mathrm{view}}){{\boldsymbol{R}}}_{z}({\varphi }_{\mathrm{view}}){\boldsymbol{v}}.\end{eqnarray} \tag{ B2 }$$
• 3.  
  At one point of the mesh grid $({v}_{\mathrm{mesh},x}^{{\prime} {ij}},{v}_{\mathrm{mesh},y}^{{\prime} {ij}},0)$, one can find the intersection points of the ejecta profile along the line of sight. We set the v_z components of each intersection point to be ${p}_{1}^{{ij}},{p}_{2}^{{ij}},{p}_{3}^{{ij}},\cdots$ from large to small, which are the velocity space distances between the intersection points and the mesh grid plane.
• 4.  
  Transform the mesh grid into the coordinate system O. The rotation rule is

  $$\begin{eqnarray}&&{\boldsymbol{v}}={{\boldsymbol{R}}}_{z}(-{\varphi }_{\mathrm{view}}){{\boldsymbol{R}}}_{y}(-{\theta }_{\mathrm{view}}){{\boldsymbol{v}}}^{{\prime} }.\end{eqnarray} \tag{ B3 }$$

  We mark the point of the mesh grid as $({v}_{\mathrm{mesh},x}^{{ij}},{v}_{\mathrm{mesh},y}^{{ij}},{v}_{\mathrm{mesh},z}^{{ij}})$.
• 5.  
  In order to solve the photosphere in the observer frame, one needs to do the line integral to find the positions where the optical depth is τ_view = 2/3 along the line of sight through each mesh grid point at the coordinate system O. The line parameter functions through one mesh grid point $({v}_{\mathrm{mesh},x}^{{ij}},{v}_{\mathrm{mesh},y}^{{ij}},{v}_{\mathrm{mesh},z}^{{ij}})$ are

  $$\begin{eqnarray}\begin{array}{rcl}{v}_{x}-{v}_{\mathrm{mesh},x}^{{ij}} & = & p\sin {\theta }_{\mathrm{view}}\cos {\varphi }_{\mathrm{view}},\\ {v}_{y}-{v}_{\mathrm{mesh},y}^{{ij}} & = & p\sin {\theta }_{\mathrm{view}}\sin {\varphi }_{\mathrm{view}},\\ {v}_{z}-{v}_{\mathrm{mesh},z}^{{ij}} & = & p\cos {\theta }_{\mathrm{view}},\end{array}\end{eqnarray} \tag{ B4 }$$

  where p is the parameter of the line functions. As for a given point (v_x, v_y, v_z) which is located on the line, p will be a velocity space distance between the given point and the mesh grid plane. The optical depth's line integral of the dynamical ejecta, neutrino-driven ejecta, and viscosity-driven ejecta are, respectively, as follows:

  $$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}{\tau }_{{\rm{v}}{\rm{i}}{\rm{e}}{\rm{w}},{\rm{d}}}^{ij} & = & {\int }_{{\rm{L}}}{\kappa }_{{\rm{d}}}{\rho }_{{\rm{d}}}dl\\ & = & \displaystyle \frac{{\kappa }_{{\rm{d}}}{M}_{{\rm{d}}}{t}^{-2}}{2{\varphi }_{{\rm{d}}}\sin {\theta }_{{\rm{d}}}({v}_{max,{\rm{d}}}-{v}_{min,{\rm{d}}})}{\int }_{p}\displaystyle \frac{1}{{v}^{2}}dp,\\ {\tau }_{{\rm{v}}{\rm{i}}{\rm{e}}{\rm{w}},{\rm{n}}}^{ij} & = & {\int }_{{\rm{L}}}{\kappa }_{{\rm{n}}}{\rho }_{{\rm{n}}}dl\\ & = & \displaystyle \frac{35{\kappa }_{{\rm{n}}}{M}_{{\rm{n}}}{t}^{-2}}{64\pi (1-\cos {\theta }_{{\rm{n}}}){v}_{max,{\rm{n}}}}{\int }_{p}\displaystyle \frac{1}{{v}^{2}}{\left[1-{\left(\displaystyle \frac{v}{{v}_{max,{\rm{n}}}}\right)}^{2}\right]}^{3}dp,\\ {\tau }_{{\rm{v}}{\rm{i}}{\rm{e}}{\rm{w}},{\rm{v}}}^{ij} & = & {\int }_{{\rm{L}}}{\kappa }_{{\rm{v}}}{\rho }_{{\rm{v}}}dl\\ & = & \displaystyle \frac{105{\kappa }_{{\rm{v}}}{M}_{{\rm{v}}}{t}^{-2}}{128\pi {v}_{max,{\rm{v}}}}{\int }_{p}\displaystyle \frac{{v}^{2}-{v}_{z}^{2}}{{v}^{4}}{\left[1-{\left(\displaystyle \frac{v}{{v}_{max,{\rm{v}}}}\right)}^{2}\right]}^{3}dp,\end{array}\end{array}\end{eqnarray} \tag{ B5 }$$

  where ${v}_{z}={v}_{\mathrm{mesh},z}^{{ij}}+p\cos {\theta }_{\mathrm{view}}$, ${v}^{2}={p}^{2}+2p({v}_{\mathrm{mesh},x}^{{ij}}\sin {\theta }_{\mathrm{view}}\cos {\varphi }_{\mathrm{view}}$ $+\,{v}_{\mathrm{mesh},y}^{{ij}}\sin {\theta }_{\mathrm{view}}\sin {\varphi }_{\mathrm{view}}+{v}_{\mathrm{mesh},z}^{{ij}}\cos {\theta }_{\mathrm{view}})+{{v}_{\mathrm{mesh}}^{{ij}}}^{2}$, and ${{v}_{\mathrm{mesh}}^{{ij}}}^{2}={{v}_{\mathrm{mesh},x}^{{ij}}}^{2}\,+{{v}_{\mathrm{mesh},y}^{{ij}}}^{2}\,+{{v}_{\mathrm{mesh},z}^{{ij}}}^{2}$.We have obtained the intersection points of the ejecta profile and their v_z components ${p}_{1}^{{ij}},{p}_{2}^{{ij}},{p}_{3}^{{ij}},\cdots$ in step 3. At a given time, one can do the line integration along the light of sight. As for the first ejecta component that the observer will see, we calculate the entire optical depth ${\tau }_{\mathrm{view}}^{{ij}}$ by assuming that the range of integration is from ${p}_{2}^{{ij}}$ to ${p}_{1}^{{ij}}$. If ${\tau }_{\mathrm{view}}^{{ij}}\gt 2/3$, the photosphere cannot penetrate this ejecta component. One can solve the photosphere's position of ${\tau }_{\mathrm{view}}^{{ij}}=2/3$ and calculate its velocity space distance ${p}_{\mathrm{phot}}^{{ij}}$ between ${p}_{1}^{{ij}}$ and ${p}_{2}^{{ij}}$. On the other hand, if ${\tau }_{\mathrm{view}}^{{ij}}\lt 2/3$, the photosphere can penetrate this ejecta component and enter the next ejecta component that the observer will see. ${p}_{\mathrm{phot}}^{{ij}}$ will be located between ${p}_{2}^{{ij}}$ and ${p}_{3}^{{ij}}$. Repeat above steps until ${p}_{\mathrm{phot}}^{{ij}}$ is solved.
• 6.  
  However, due to the light-propagation effect, photons emitted at a given time t reach the observer at different arrival times t_obs. We set the arrival time t_obs = t when a photon is emitted at the mesh grid plane. Therefore, the observer would see a photon emitted from the photosphere at

  $$\begin{eqnarray}&&{t}_{\mathrm{obs}}=t\left(1-\displaystyle \frac{{p}_{\mathrm{phot}}^{{ij}}}{c}\right).\end{eqnarray} \tag{ B6 }$$
• 7.  
  At a given observational time t_obs, the photosphere location assumed to be $({v}_{\mathrm{phot},x}^{{ij}},{v}_{\mathrm{phot},y}^{{ij}},{v}_{\mathrm{phot},z}^{{ij}})$ can be solved using Equation (B4) if ${p}_{\mathrm{phot}}^{{ij}}$ have been solved in the coordinate system O'. We have obtained the thermal temperature ${T}_{\mathrm{mesh}}^{{ij}}$ in Section 2 when the observational time t_obs and the observed photosphere location are known. Throughout all the points of the mesh grid, one can obtain the entire photosphere in the observer frame and the thermal temperature at the photosphere. A schematic diagram is presented in Figure B1. In order to calculate the observed flux, it is convenient to directly do the interpolation and integration by projecting the parameter information of the entire photosphere onto the mesh grid if the source angular size can be ignored. The observed flux for a given frequency ν is

  $$\begin{eqnarray}\begin{array}{rcl}{F}_{\nu }(\nu ,{t}_{\mathrm{obs}}) & \approx & \displaystyle \frac{2}{{h}^{2}{c}^{2}{D}_{{\rm{L}}}^{2}}\\ & & \times \ {\displaystyle \iint }_{S}\displaystyle \frac{{ \mathcal D }{\left(h\nu /{ \mathcal D }\right)}^{3}}{\exp (h\nu /{ \mathcal D }{k}_{{\rm{B}}}{T}_{\mathrm{mesh}}^{{ij}})-1}d{\sigma }^{{\prime} {ij}},\end{array}\end{eqnarray} \tag{ B7 }$$

  where D_L is the luminosity distance, h is the Planck constant, k_B is the Boltzmann constant, and $d{\sigma }^{{\prime} {ij}}$ is the infinitesimal projected photosphere area, i.e.,

  $$\begin{eqnarray}&&d{\sigma }^{{\prime} {ij}}={\left(\displaystyle \frac{{t}_{\mathrm{obs}}}{1-{p}_{\mathrm{phot}}^{{ij}}/c}\right)}^{2}{{dv}}_{\mathrm{mesh},x}^{{\prime} {ij}}{{dv}}_{\mathrm{mesh},y}^{{\prime} {ij}}.\end{eqnarray} \tag{ B8 }$$

  The Doppler factor is ${ \mathcal D }=1/[{\rm{\Gamma }}(1-\beta \cos {\rm{\Delta }}\theta )]$, where ${\rm{\Gamma }}=1/\sqrt{1-{\beta }^{2}}$, β = v/c, and Δθ is the angle between the moving direction of the point located at the photosphere and the line of sight. In the coordinate system O', the moving direction of the point located at the photosphere is $({v}_{\mathrm{mesh},x}^{{\prime} {ij}},{v}_{\mathrm{mesh},y}^{{\prime} {ij}},{p}_{\mathrm{phot}}^{{ij}})$ while the line of sight is toward to the ${v}_{z}^{{\prime} }$ direction. Following the included angle formula, one can calculate $\cos {\rm{\Delta }}\theta \,={p}_{\mathrm{phot}}^{{ij}}/\sqrt{{{v}_{\mathrm{mesh},x}^{{\prime} {ij}}}^{2}+{{v}_{\mathrm{mesh},y}^{{\prime} {ij}}}^{2}+{{p}_{\mathrm{phot}}^{{ij}}}^{2}}$. The conversion between frequency and wavelength is obtained from F_νdν = F_λdλ. Therefore, the observed flux for a given wavelength λ is

  $$\begin{eqnarray}&&{F}_{\lambda }(\lambda ,{t}_{\mathrm{obs}})\approx \displaystyle \frac{2}{{h}^{4}{c}^{3}{D}_{{\rm{L}}}^{2}}{\displaystyle \iint }_{S}\displaystyle \frac{{{ \mathcal D }}^{3}{\left({hc}/{ \mathcal D }\lambda \right)}^{5}}{\exp ({hc}/{ \mathcal D }{k}_{{\rm{B}}}\lambda {T}_{\mathrm{mesh}}^{{ij}})-1}d{\sigma }^{{\prime} {ij}}.\end{eqnarray} \tag{ B9 }$$

Zoom In Zoom Out Reset image size

In this section, we briefly describe the afterglow model we used to calculate the sGRB lightcurves along the line of sight.

After producing the sGRB, the relativistic jet sweeps into the ISM, driving a forward shock. The mass swept per unit solid angle at a given radius R is

$$\begin{eqnarray}&&\mu (R)=\displaystyle \frac{{R}^{3}}{3}{{nm}}_{p},\end{eqnarray} \tag{ C1 }$$

where n is the ISM number density and m_p is the proton mass. By assuming energy conservation, the Lorentz factor of the shocked material (Panaitescu & Kumar 2000; Granot & Kumar 2003) is

$$\begin{eqnarray}&&{\rm{\Gamma }}(R,\theta )=\displaystyle \frac{{\mu }_{0}}{2\mu }\left[\sqrt{1+\tfrac{4\mu ({c}^{-2}{dE}/d{\rm{\Omega }}+\mu +{\mu }_{0})}{{\mu }_{0}^{2}}}-1\right],\end{eqnarray} \tag{ C2 }$$

where ${\mu }_{0}(\theta )=({dE}/d{\rm{\Omega }})/[{\rm{\Gamma }}(0,\theta )-1]{c}^{2}$, dE/dΩ is the kinetic energy angular distribution of the jet, and c is the speed of light.

For a relativistic shock propagating into a cold ISM, the physical condition of the shocked plasma is obtained from the shock-jump conditions (Blandford & McKee 1976), i.e., the conservations of baryon number and energy and momentum fluxes. Based on these conservation conditions, the electron number density n_s of the shocked material is

$$\begin{eqnarray}&&{n}_{{\rm{s}}}=\displaystyle \frac{{\gamma }_{\mathrm{ad}}{\rm{\Gamma }}+1}{{\gamma }_{\mathrm{ad}}-1}n,\end{eqnarray} \tag{ C3 }$$

while its Lorentz factor is

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{s}}}=[{\gamma }_{\mathrm{ad}}({\rm{\Gamma }}-1)+1]\sqrt{\displaystyle \frac{{\rm{\Gamma }}+1}{{\gamma }_{\mathrm{ad}}(2-{\gamma }_{\mathrm{ad}})({\rm{\Gamma }}-1)+2}},\end{eqnarray} \tag{ C4 }$$

where γ_ad is the post-shock adiabatic index. We adopt the γ_ad, which is a function of Γ, by Pe'er (2012): ${\gamma }_{\mathrm{ad}}=(5-1.21937z+0.18203{z}^{2}-0.96583{z}^{3}+2.32513{z}^{4}\,-2.39332{z}^{5}+1.07136{z}^{6})/3$ where z = Θ/(0.24 + Θ) and

$$\begin{eqnarray}&&{\rm{\Theta }}=\left(\displaystyle \frac{{\rm{\Gamma }}\beta }{3}\right)\left(\displaystyle \frac{{\rm{\Gamma }}\beta +1.07{({\rm{\Gamma }}\beta )}^{2}}{1+{\rm{\Gamma }}\beta +1.07{({\rm{\Gamma }})\beta )}^{2}}\right).\end{eqnarray} \tag{ C5 }$$

By assuming that the shock material is concentrated in a thin layer behind the shock and has a uniform radial density distribution, the thickness of the shocked layer is given by Salafia et al. (2019),

$$\begin{eqnarray}&&{\rm{\Delta }}R=\displaystyle \frac{R({\gamma }_{\mathrm{ad}}-1)}{3({\gamma }_{\mathrm{ad}}+1){\rm{\Gamma }}}.\end{eqnarray} \tag{ C6 }$$

The shock surface brightness is

$$\begin{eqnarray}&&{I}_{\nu }(\nu ,R,\theta ,\phi )={{ \mathcal D }}^{3}{\rm{\Delta }}{R}^{{\prime} }{j}_{{\nu }^{{\prime} }}^{{\prime} }(\nu /{ \mathcal D }),\end{eqnarray} \tag{ C7 }$$

where the Doppler factor is ${ \mathcal D }(R,\theta ,\varphi ,{\theta }_{\mathrm{view}})\,={{\rm{\Gamma }}(R,\theta )}^{-1}{[1-\beta (R,\theta )\cos \alpha ]}^{-1}$, β = (1 − Γ⁻²)^1/2 and ${\rm{\Delta }}{R}^{{\prime} }\,={\rm{\Gamma }}(R,\theta ){\rm{\Delta }}R$, and ${j}_{{\nu }^{{\prime} }}^{{\prime} }$ is the comoving emissivity due to synchrotron emission. The synchrotron spectrum for a distribution of electrons depends on the ordering of characteristic break frequencies, i.e., the typical synchrotron frequency of the accelerated electrons with the minimum Lorentz factor ${\nu }_{{\rm{m}}}^{{\prime} }$, the cooling frequency ${\nu }_{{\rm{c}}}^{{\prime} }$, and the synchrotron self-absorption frequency ${\nu }_{{\rm{a}}}^{{\prime} }$. The spectrum for ${\nu }_{{\rm{m}}}^{{\prime} }\lt {\nu }_{{\rm{c}}}^{{\prime} }$ is classified as the "slow cooling" case, i.e., (Sari et al. 1998)

$$\begin{eqnarray}{j}_{{\nu }^{{\prime} }}^{{\prime} }={j}_{{\nu }^{{\prime} },\max }^{{\prime} }\left\{\begin{array}{ll}{\left(\displaystyle \frac{{\nu }_{{\rm{a}}}^{{\prime} }}{{\nu }_{{\rm{m}}}^{{\prime} }}\right)}^{1/3}{\left(\displaystyle \frac{{\nu }^{{\prime} }}{{\nu }_{{\rm{a}}}^{{\prime} }}\right)}^{2}, & {\nu }^{{\prime} }\lt {\nu }_{{\rm{a}}}^{{\prime} }\\ {\left(\displaystyle \frac{{\nu }^{{\prime} }}{{\nu }_{{\rm{m}}}^{{\prime} }}\right)}^{1/3}, & {\nu }_{{\rm{a}}}^{{\prime} }\lt {\nu }^{{\prime} }\lt {\nu }_{{\rm{m}}}^{{\prime} }\\ {\left(\displaystyle \frac{{\nu }^{{\prime} }}{{\nu }_{{\rm{m}}}^{{\prime} }}\right)}^{-(p-1)/2}, & {\nu }_{{\rm{m}}}^{{\prime} }\lt {\nu }^{{\prime} }\lt {\nu }_{{\rm{c}}}^{{\prime} }\\ {\left(\displaystyle \frac{{\nu }_{{\rm{c}}}^{{\prime} }}{{\nu }_{{\rm{m}}}^{{\prime} }}\right)}^{-(p-1)/2}{\left(\displaystyle \frac{{\nu }^{{\prime} }}{{\nu }_{{\rm{c}}}^{{\prime} }}\right)}^{-p/2}, & {\nu }^{{\prime} }\gt {\nu }_{{\rm{c}}}^{{\prime} }\end{array}\right.\end{eqnarray} \tag{ C8 }$$

and that for ${\nu }_{{\rm{c}}}^{{\prime} }\lt {\nu }_{{\rm{m}}}^{{\prime} }$, the "fast cooling" case, is

$$\begin{eqnarray}{j}_{{\nu }^{{\prime} }}^{{\prime} }={j}_{{\nu }^{{\prime} },\max }^{{\prime} }\left\{\begin{array}{ll}{\left(\displaystyle \frac{{\nu }_{{\rm{a}}}^{{\prime} }}{{\nu }_{{\rm{c}}}^{{\prime} }}\right)}^{1/3}{\left(\displaystyle \frac{{\nu }^{{\prime} }}{{\nu }_{{\rm{a}}}^{{\prime} }}\right)}^{2}, & {\nu }^{{\prime} }\lt {\nu }_{{\rm{a}}}^{{\prime} }\\ {\left(\displaystyle \frac{{\nu }^{{\prime} }}{{\nu }_{{\rm{c}}}^{{\prime} }}\right)}^{1/3}, & {\nu }_{{\rm{a}}}^{{\prime} }\lt {\nu }^{{\prime} }\lt {\nu }_{{\rm{c}}}^{{\prime} }\\ {\left(\displaystyle \frac{{\nu }^{{\prime} }}{{\nu }_{{\rm{c}}}^{{\prime} }}\right)}^{-1/2}, & {\nu }_{{\rm{c}}}^{{\prime} }\lt {\nu }^{{\prime} }\lt {\nu }_{{\rm{m}}}^{{\prime} }\\ {\left(\displaystyle \frac{{\nu }_{{\rm{m}}}^{{\prime} }}{{\nu }_{{\rm{c}}}^{{\prime} }}\right)}^{-1/2}{\left(\displaystyle \frac{{\nu }^{{\prime} }}{{\nu }_{{\rm{m}}}^{{\prime} }}\right)}^{-p/2}, & {\nu }^{{\prime} }\gt {\nu }_{{\rm{m}}}^{{\prime} }.\end{array}\right.\end{eqnarray} \tag{ C9 }$$

In our calculations, we set p = 2.2. The comoving synchrotron emissivity of electrons behind the shock at the peak of the spectrum is

$$\begin{eqnarray}&&{j}_{\nu ^{\prime} ,\max }^{{\prime} }\approx 0.66\displaystyle \frac{{q}_{e}^{3}}{{m}_{e}^{2}{c}^{4}}\displaystyle \frac{p-2}{3p-1}\displaystyle \frac{B^{\prime} {\epsilon }_{e}e}{{\gamma }_{{\rm{m}}}},\end{eqnarray} \tag{ C10 }$$

where $e=({\rm{\Gamma }}-1){n}_{{\rm{s}}}{m}_{p}{c}^{2}$, m_e is the mass of electrons, q_e is the electron charge, and ${\epsilon }_{e}$ is the fraction of internal energy that is given to electrons. The comoving magnetic field strength $B^{\prime}$ is

$$\begin{eqnarray}&&B^{\prime} ={\left(32\pi {m}_{p}{\epsilon }_{B}n\right)}^{1/2}{\rm{\Gamma }}c,\end{eqnarray} \tag{ C11 }$$

where ${\epsilon }_{B}$ is the internal energy that goes to magnetic fields.

More specifically, for each characteristic frequency, ${\nu }_{{\rm{m}}}^{{\prime} }$ is the synchrotron frequency, which is ${\nu }_{{\rm{m}}}^{{\prime} }={\gamma }_{{\rm{m}}}^{2}{q}_{e}B^{\prime} /2\pi {m}_{e}c$, where m_e is the electron mass,

$$\begin{eqnarray}&&{\gamma }_{{\rm{m}}}=\max \left[1,\displaystyle \frac{p-2}{p-1}({\rm{\Gamma }}-1)\displaystyle \frac{{m}_{p}}{{m}_{e}}\right],\end{eqnarray} \tag{ C12 }$$

${\nu }_{{\rm{c}}}^{{\prime} }$ is the synchrotron frequency corresponding to the Lorentz factor ${\gamma }_{{\rm{c}}}$, which is

$$\begin{eqnarray}&&{\gamma }_{{\rm{c}}}=\displaystyle \frac{6\pi {m}_{e}{c}^{2}{\rm{\Gamma }}\beta }{{\sigma }_{{\rm{T}}}{B}^{{\prime} 2}R},\end{eqnarray} \tag{ C13 }$$

where σ_T is the Thomson cross section. The frequency ${\nu }_{{\rm{a}}}^{{\prime} }$ is the synchrotron self-absorption, which we calculate based on Shen & Zhang (2009):

$$\begin{eqnarray}&&2{k}_{{\rm{B}}}T^{\prime} \displaystyle \frac{{\nu }_{{\rm{a}}}^{{\prime} 2}}{{c}^{2}}=C(p){I}_{{\nu }_{{\rm{a}}}^{{\prime} }}^{{\prime} },\end{eqnarray} \tag{ C14 }$$

where k_B is the Boltzmann constant and ${k}_{{\rm{B}}}{T}^{{\prime} }=\max ({\gamma }_{{\rm{a}}},\min ({\gamma }_{{\rm{m}}},{\gamma }_{{\rm{c}}})){m}_{e}{c}^{2}$. The correction factor is

$$\begin{eqnarray}C(p)=\left\{\begin{array}{ll}{C}_{1}(p)=\tfrac{(p+2)(p-1/3)}{p+2/3}, & \mathrm{if}\ {\nu }_{{\rm{a}}}^{{\prime} }\lt \min ({\nu }_{{\rm{m}}}^{{\prime} },{\nu }_{{\rm{c}}}^{{\prime} })\\ {C}_{2}(p)=\sqrt{2}(p+1)\tfrac{{\rm{\Gamma }}\left(\tfrac{3p+22}{12}\right){\rm{\Gamma }}\left(\tfrac{3p+2}{12}\right)}{{\rm{\Gamma }}\left(\tfrac{3p+19}{12}\right){\rm{\Gamma }}\left(\tfrac{3p-1}{12}\right)}, & \mathrm{if}\ \min ({\nu }_{{\rm{m}}}^{{\prime} },{\nu }_{{\rm{c}}}^{{\prime} })\lt {\nu }_{{\rm{a}}}^{{\prime} }\lt \max ({\nu }_{{\rm{m}}}^{{\prime} },{\nu }_{{\rm{c}}}^{{\prime} }),\end{array}\right.\end{eqnarray} \tag{ C15 }$$

where Γ denotes the gamma function in this function.

Due to the aberration effect, photons emitted at (R, θ, φ) arrive at the observer with a time delay with respect to a photon emitted at R = 0 by

$$\begin{eqnarray}&&{t}_{\mathrm{obs}}={\int }_{0}^{R}\displaystyle \frac{{dR}(1-{\beta }_{{\rm{s}}}\cos \alpha )}{{\beta }_{{\rm{s}}}c},\end{eqnarray} \tag{ C16 }$$

where ${\beta }_{{\rm{s}}}={(1-{{\rm{\Gamma }}}_{{\rm{s}}}^{-2})}^{1/2}$ and $\cos \alpha =\cos \theta \cos {\theta }_{\mathrm{view}}\,+\sin \theta \sin \varphi \sin {\theta }_{\mathrm{view}}$. By integrating the equal-arrival time surfaces, the flux density is given by

$$\begin{eqnarray}&&{F}_{\nu }(\nu ,{t}_{\mathrm{obs}})=\displaystyle \frac{2}{{D}_{{\rm{L}}}^{2}}{\int }_{0}^{1}d\,\cos \theta {\int }_{-\pi /2}^{\pi /2}d\varphi {R}^{2}{I}_{\nu }(\nu ,R).\end{eqnarray} \tag{ C17 }$$

In order to calculate the flux density from the counter-jet, one can set θ_view to θ_view + π.