Long-duration Gamma-Ray Burst and Associated Kilonova Emission from Fast-spinning Black Hole–Neutron Star Mergers

We collect the observed data of three present z-known merger-origin lGRBs, including GRB 211221A (e.g., Rastinejad et al. 2022; Xiao et al. 2022; Yang et al. 2022), GRB 060614 (Della Valle et al. 2006; Gal-Yam et al. 2006; Zhang et al. 2007), and GRB 211227A (Lü et al. 2022), all characterized by a spiky long-duration ME phase, a rebrightening EE phase, and a temporal lull between these two phases, to study their similarities. Here, the EE phase is defined as the long-lasting lower-level emission phase after the initial intense ME phase (Norris & Bonnell 2006; Lan et al. 2020). The whole emission (WE) phase includes the ME phase and the EE phase. Table 1 lists the observed properties of the ME and the WE phases.

Table 1. The Observed Properties Of GRB 211211A, GRB 060614, ^a and GRB 211227A ^b

	GRB 211211A	GRB 060614	GRB 211227A
ME
Duration (s)	13	6	4
Peak energy (keV)
Energy fluence (erg cm−2)
Isotropic equivalent energy (erg)
Spectral index α
Spectral index β		⋯	⋯
WE
Duration (s)		102 ± 5	84
Peak energy (keV)		10 − 100
Energy fluence c (erg cm−2)
Isotropic equivalent energy (erg)
Spectral index α		⋯
Spectral index β		⋯
Redshift	0.076	0.125	0.228

Notes.

^a The data of GRB 211221A and GRB 060614 are collected from Table 1 of Yang et al. (2022). ^b The data of GRB 211227A are collected from Lü et al. (2022) and Tsvetkova et al. (2022). ^c The energy fluence of the ME for GRB 211227A is calculated by the HEASoft tool (Nasa High Energy Astrophysics Science Archive Research Center (Heasarc), 2014), in the 15–1500 keV energy band.

Download table as:  ASCIITypeset image

In the left panel of Figure 1, we divide the detected GRBs into two populations, i.e., lGRB and sGRB populations, in the duration T₉₀ versus Earth-frame peak energy E_p diagram, through the Gaussian mixture model. The intermediate population at the boundary region between the lGRB and sGRB populations could originate from collaspars or from compact-object coalescences (e.g., Tunnicliffe & Levan 2012; Zaninoni et al. 2016). Both the MEs and WEs of these three merger-origin GRBs fall into the distribution of lGRBs. Therefore, without the redshift, host galaxy information, and associated kilonova detection, these three bursts would be classified as members of the lGRB population, due to their long durations.

Zoom In Zoom Out Reset image size

The correlations between the total isotropic equivalent energy E_γ,iso of the prompt emission and the rest-frame peak energy E_p(1 + z) (the Amati relation; Amati et al. 2002) for both the lGRB and sGRB populations are shown in the right panel of Figure 1. We find that the WEs of these three bursts still behave as normal lGRBs, although their E_γ,iso are lower than those of most observed lGRBs. Conversely, their MEs lie on the middle location of the sGRB track, rather than on the lGRB track. The differences in E_γ,iso between the WEs and MEs, corresponding to the E_γ,iso of the EEs, for these bursts are similar, which are ∼1–2 × 10⁵¹ erg.

After a BNS/disrupted NSBH merger, a compact remnant accretion disk would be formed around the NS or BH. The lifetime of this disk is typically ≲1 s, which is thought to determine the burst duration of an sGRB (e.g., Shapiro 2017; Zhang 2019; Ruiz et al. 2020). In order to explain the long durations of these three merger-origin GRBs, an additional energy/matter injection is needed. When a group of bound ejecta with an energy distribution of dM/dE ∝ E^α , where α ≈ 0 around E = 0, fall back onto the central NS or BH, the fallback rate would track as ∝ t^−5/3 (Rees 1988). As shown in Figure 2, we find that there are long-duration fallback accretion signals appearing in the X-ray light curves of both the MEs and EEs for these three bursts. The fallback rates in the EEs peak at ∼20–40 s and follow ∝ t^−5/3, with a duration of ∼100–200 s. Thus, the fallback accretion might result in their long durations.

Zoom In Zoom Out Reset image size

Gompertz et al. (2020a) collected 39 z-known sGRBs (seven events had an EE), and searched for possible NSBH candidates among these sGRBs. However, they did not find any clear evidence for the existence of NSBH GRBs in this complete sGRB sample. If a millisecond magnetar can survive after a BNS merger, the magnetar can lose its rotational energy via the spin-down process, resulting in a long-duration X-ray plateau in the X-ray light curve (Zhang & Mészáros 2001). Gompertz et al. (2020a) suggested that four of the sGRBs with EEs in their sample could be magnetar-like, so they could not originate from NSBH mergers. In Figure 3, we collect the X-ray light curves of the three other sGRBs with EEs, including GRB 061006 (Berger 2007), GRB 061210 (Berger 2007), and GRB 071227 (D'Avanzo et al. 2009), whose X-ray emissions could not be well interpreted via the magnetar model by Gompertz et al. (2020a). We find that GRB 071227 could still be magnetar-like, being characterized by a long-duration X-ray plateau. The X-ray light curves of GRB 061006 and GRB 061210 do not show any signatures of ∝ t^−5/3 decay, and hence they may provide no evidence of an unambiguous fallback accretion signal, though that could be partly due to their limited data points. Thus, the EEs in the present z-known sGRBs may not be caused by the fallback accretions. The similar observed properties of merger-origin lGRBs, i.e., long durations and unambiguous rebrightening fallback accretion signals, do not present in the existing observations of the merger-origin sGRBs with/without EEs, indicating a unique origin for them.

Zoom In Zoom Out Reset image size

In principle, both BNS and NSBH mergers can generate early-time fallback accretions (Rosswog 2007; Metzger et al. 2010a; Fernández et al. 2017). However, most of the disrupted NSBH mergers would eject a lot more materials and lead to more powerful fallback accretions than BNS mergers (Rosswog 2007). Furthermore, the simulations by Metzger et al. (2010a) and Desai et al. (2019) showed that r-process heating might affect the materials that are marginally gravitationally bound. The gravity of the remnants formed after BNS mergers might be too low to drag these r-process heating materials back and result in fallback accretions. However, as predicted by Desai et al. (2019), the late-time fallback accretion of these materials may occur after tens of seconds of the merger, if the remnant BH has a mass of ≳ 6–8 M_⊙, leading to a rebrightening emission with a t^−5/3 power-law decay appearing in the X-ray light curve. The EEs of these three bursts, whose starting times and observational features are consistent with the predictions of Desai et al. (2019), may originate from the fallback accretion of r-process heating materials. We thus suspect that these three bursts are derived from NSBH mergers.

Assuming the unique merger origins of these three lGRBs, we investigate their local event-rate density by tentatively constructing their luminosity function (LF), following Sun et al. (2015, 2022). The limited number of events could lead to the large uncertainty of the LF. However, it could be plausible to make a fine evaluation by adding up the contributions of all three events. Following Equation (12) in Sun et al. (2022), the LF Φ(L) is the sum of each event j in the range from to (with a total of ΔN_L events):

where T_BAT is the total monitoring time by Swift/Burst Alert Telescope (BAT) and is the effective maximum volume that is monitored (with the field of view Ω_BAT), weighted by the density evolution f(z) and time dilation, i.e.,

The is the maximum redshift at which an event can be detected with bolometric peak luminosity L_bol. The bolometric peak luminosities of these three events, defined in the energy range of 1–10⁴ keV, are derived from k-correction based on the spectrum listed in Table 1, following Equation (29) in Sun et al. (2015).

We take the instrument parameters of Swift/BAT, with Ω_BAT = 1.3, T_BAT = 17 yr, and flux sensitivity f_th,BAT = 3 ×10⁻⁸ erg cm⁻² s⁻¹. The f(z) for GRBs with merger origins is adopted from Zhu et al. (2021b). We take the correction of the redshift measurement into account. For both the lGRB and sGRB samples, the ratio between the total number and the number of z-known events from the Swift observations is approximately 4: 1. In addition, since Swift has a relatively softer energy band than BATSE, it tends to detect more lGRBs than sGRBs. One needs to adopt another correction, by a factor of 3, to compensate for the short (hard)-to-long (soft) ratio of the Swift-detected GRBs in comparison with those from the BATSE observations (Sun et al. 2015).

We divide the events into two bins, with , and fit the LF with a single power law of the slope 0.6 ± 0.8. The local event-rate density, which is derived by integrating the LF above the luminosity threshold 3 × 10⁵⁰ erg s⁻¹, is given as

The errors are given at the 1σ confidence level (Gehrels 1986).

The intrinsic local event-rate densities of both lGRBs and sGRBs are of the order of unity, in the unit of Gpc⁻³ yr⁻¹, above the isotropically bolometric luminosity of ∼10⁵⁰ erg s⁻¹ (Sun et al. 2015). We find that the local event-rate density for lGRBs with merger origins is much lower than those for both lGRBs and sGRBs. The rate densities for BNSs inferred using GWs through GWTC-2 (GWTC-3) are (10–1700 Gpc⁻³ yr⁻¹) (e.g., Abbott et al. 2021b; the LIGO Scientific Collaboration et al. 2021; Mandel & Broekgaarden 2022). Abbott et al. (2021a) inferred an NSBH rate density of , by considering the observations of two NSBH mergers, or , assuming a broad NSBH population, while GWTC-3 (the LIGO Scientific Collaboration et al. 2021) reported the NSBH merger rate density as being between 7.8 and 140 Gpc⁻³ yr⁻¹. By considering a jet beaming factor of f_b = 0.01, the beaming-corrected event-rate density for merger-origin lGRBs is , which is much lower than those for BNS and NSBH mergers. However, the beaming-corrected rate density for merger-origin lGRBs is consistent with that for NS-first-born NSBH mergers (≲20% of the total NSBH populations; Román-Garza et al. 2021; Chattopadhyay et al. 2022). By investigating the parameter space for forming NS-first-born NSBH binaries, Hu et al. (2022) found that most of the NSBH binaries that could merge within Hubble time would have BHs with projected aligned spins χ_BH ≳ 0.8, and, hence, could certainly make the tidal disruptions to produce EM counterparts. Only a small fraction of low-mass BHs with χ_BH ∼ 0.2–0.8 could merge with an NS within Hubble time, and this could still allow tidal disruption to happen, if the NSs are not really massive (i.e., M_NS ≳ 1.6 − 2.0 M_⊙). Since the rest of the BH-first-born NSBH mergers would mostly contribute to plunging events, due to the negligible projected aligned spins of the BH components (Zhu et al. 2022), a single explosive population via an NS-first-born NSBH merger could account for their origin.

3.1.1. Fallback Accretion of r-process Heating Materials

For NSBH mergers with NS tidal disruption, a fraction of the r-process heating materials would fall back onto the remnant BH, tens of seconds after the merger, resulting in the EE (Metzger et al. 2010a; Desai et al. 2019), through the Blandford–Znajek (BZ) mechanism (Blandford & Znajek 1977). Following MacFadyen et al. (2001) and Dai & Liu (2012), the fallback rate initially increases with time, as , before the time of t_p corresponds to the peak fallback rate . Then, the late-time fallback accretion behavior would track as , until the break time t_b (Chevalier 1998). While most of the fallback materials are accreted after t_b, we describe the fallback rate as . An empirical three-segment broken power-law function is adopted to model the fallback accretion rate of the r-process heating materials, i.e.,

where t₀ is the starting time of the fallback accretion.

The BZ power is related to the mass and spin of the central BH (e.g., Li & Paczyński 2000; Lei et al. 2013, 2017; Wu et al. 2013; Liu et al. 2017), i.e.,

where M_BH is the central BH mass, χ_BH is the dimensionless aligned spin of the BH, B_BH,15 = B_BH/10¹⁵ G is the magnetic field strength threading the BH horizon, and with . The magnetic field can be estimated through the balance between the magnetic pressure on the horizon and the ram pressure of the innermost part of the accretion flow:

where c is the speed of light and is the radius of the BH horizon, with the Schwarzschild radius r_g = GM_BH/c² and the gravitational constant G. Thus, the BZ power can be also expressed as

Because the central BH would be spun up by accretion and spun down by the BZ mechanism, the conservation of the energy and angular momentum of a BH can be written as

where Ω_H = c χ_BH/(2r_H) is the angular velocity of the BH horizon and and are the specific energy and specific angular momentum of a particle at the innermost stable circular orbit (ISCO) radius (Novikov & Thorne 1973), respectively. represents the normalized radius of the BH ISCO with and (Bardeen et al. 1972). Since the angular momentum of BH is expressed as , one has

By combing Equations (7), (8), and (9), the time-evolving BZ power can thus be calculated. The observed gamma-ray/X-ray light curve caused by the fallback accretion is connected to the BZ power via the gamma-ray/X-ray radiation efficiency η_(γ,X) and the jet beaming factor f_b, i.e.,

3.1.2. Jet Afterglow Emissions

In order to calculate the afterglow light curves, we adopt the Gaussian structured jet model (e.g., Zhang & Mészáros 2002), which was favored by the observations of the GRB 170817A afterglow (e.g., Lamb & Kobayashi 2018; Lazzati et al. 2018; Mooley et al. 2018; Troja et al. 2018; Xie et al. 2018), i.e.,

where E₀ is the on-axis equivalent isotropic energy and θ_c is the characteristic core angle. The spectra of the standard synchrotron emissions from relativistic electrons are employed following Sari et al. (1998), Kumar & Zhang (2015), and Zhang (2018). For more details of the afterglow modeling that we apply to calculate the sGRB light curves along the line of sight, see Appendix C in Zhu et al. (2021b). We constrain the afterglow parameters, including E₀, θ_c, the viewing angle θ_v with respect to the moving direction of the jet, the circumburst number density n, the power-law index of the electron distribution p, and the fraction of the shock energy carried by the magnetic fields ε_B , to fit the multiband light curves of GRB 211211A. The fraction of the shock energy carried by electrons is set to its typical value of ε_e = 0.1.

3.1.3. Ejecta Mass

After NSBH mergers, a fraction of the neutron-rich matter (i.e., an unbound dynamical ejecta) is tidally ejected, while an accretion disk is formed around the remnant BH. The total remnant mass outside the remnant and the dynamical ejecta mass is dependent on the NSBH system parameters, including the BH mass M_BH, the dimensionless spin parameter projected onto the orientation of the orbital angular momentum χ_BH, the NS mass M_NS, and the NS equation of state (EoS), which can be calculated based on an empirical fitting formula, i.e.,

where is the baryonic mass of the NS, C_NS is the compactness of the NS determined by the NS EoS, η = Q/(1 + Q)², and Q = M_BH/M_NS is the mass ratio between the primary BH mass and the secondary NS mass. For the fitting formula of the total remnant mass M_tot,fit (the dynamical ejecta mass M_d,fit), the parameters in Equation (12) are a₁ = 0.406, a₂ = 0.139, a₃ = 0.255, and a₄ = 1.761 (a₁ = 0.218, a₂ = 0.028, a₃ = − 0.122, and a₄ = 1.358), obtained from Foucart et al. (2018) (Zhu et al. 2020). Since the fitting formulas of the total remnant mass and the dynamical ejecta mass are obtained with independent simulation data, one needs to set an upper limit on the maximum fraction of dynamical ejecta mass to total remnant mass, i.e., . We set , based on the simulation results from Kyutoku et al. (2015). Therefore, the final empirical mass of the dynamical ejecta is .

We consider two ejecta components for NSBH kilonova model, i.e., the wind ejecta from the disk around the remnant BH and the dynamical ejecta caused by tidal forces. The wind ejecta mass can be estimated as a constant fraction of the disk mass, i.e., M_w ≈ ξ_w M_disk, where ξ_w ≈ 0.2 (Fernández et al. 2015; Just et al. 2015; Siegel & Metzger 2017).

Numerical simulations have revealed that the dynamical ejecta from NSBH mergers are highly anisotropic and distributed in the equatorial plane (Kyutoku et al. 2015; Kawaguchi et al. 2016; Darbha et al. 2021). Zhu et al. (2020) constructed a viewing angle–dependent model for NSBH kilonovae, and found that the wind ejecta can be covered by the dynamical ejecta for a large θ_v condition. However, for the case of GRB 211211A, observed from the on-axis or near-on-axis view (i.e., θ_v ∼ 0°), one can simultaneously see two components. Hereafter, in order to reduce the computational complexity, we use a simplified model based on Zhu et al. (2020), to separately consider the emissions from the wind ejecta and the dynamical ejecta.

Assuming that the wind ejecta has an isotropic density profile and a homologous expansion, we adopt the common analytic solution derived by Arnett (1982) and Chatzopoulos et al. (2012) for calculating the bolometric luminosity of the wind ejecta:

where is the photon diffusion timescale of the wind ejecta, κ_w is the gray opacity, and β = 13.8 is the dimensionless constant. represents the injection heating rate from the radioactive decay of the r-process nucleus, where , if Y_e ≤ 0.25 ( otherwise) is an electron fraction–dependent term that takes into account extremely neutron-rich ejecta with a decay half-life of a few hours (Perego et al. 2017), _th ≈ 0.5 is the efficiency of thermalization (Metzger et al. 2010b), and is the specific energy injection rate due to radioactive decay, with and α ≈ 1.3 (Korobkin et al. 2012).

In order to calculate the monochromatic light curves of the wind ejecta, we define a photosphere temperature as

where σ_SB is the Stefan–Boltzmann constant, v_w ≈ 0.067 c is assumed to be the wind ejecta velocity (e.g., Just et al. 2015; Perego et al. 2017; Siegel & Metzger 2017), and T_La ≈ 1000 K is the first ionization temperature of lanthanides (Barnes & Kasen 2013). The photosphere radius can be written as

then the flux density contributed from the emission of the wind ejecta is given by

where h, k_B, and ν represent the Planck constant and the Boltzmann constant and frequency, respectively.

Based on the simulations of NSBH mergers (Kyutoku et al. 2015), the mass distribution of the dynamical ejecta is highly anisotropic, with the mass mainly distributed around the equatorial plane and shaped like a crescent. The dynamical ejecta is typically concentrated around the orbital plane, with a half-opening angle in the latitudinal direction of θ_d ≈ 15°, and often sweeps out only a half of the plane, i.e., an opening angle in the longitudinal direction of φ_d ≈ 180°. Since the dynamical ejecta is geometrically thin in the latitudinal direction, the photons would always be diffused from the latitudinal edge. Due to the between the radial velocity range of , based on the numerical relativity simulations (Kyutoku et al. 2015), the bolometric luminosity of the dynamical ejecta can be obtained from Kawaguchi et al. (2016), i.e.,

where is defined as the critical diffuse timescale on which all ejecta can be seen. We set the minimum velocity , following the simulations of Kyutoku et al. (2015), while the maximum velocity can be estimated as , where v_d,rms = ( − 0.441Q^−0.224 + 0.539)c is the root-mean-square velocity of the dynamical ejecta obtained from Zhu et al. (2020).

Because dL_d/2 is released over an area of φ_d vt² dv, one can derive the photosphere temperature of each velocity at a given time:

The total flux density from the dynamical ejecta can be expressed as

For simplicity, we directly use the final BH mass M_BH,f and the final dimensional aligned spin χ_BH,f after NSBH mergers, based on Equation (3) and Equation (5) from Deng (2020), as functions of the initial NSBH system parameters, i.e., M_BH, M_NS, C_NS, and χ_BH, to determine the BZ power. The observed mass distributions of Galactic BNS systems (Lattimer 2012) were inferred to be Gaussian distributions, i.e., . The observations of GW200105 and GW200115 showed that their NS masses are ∼1.9 M_⊙ and ∼1.5 M_⊙ (Abbott et al. 2021a), plausibly more massive than the mass distributions of Galactic BNS systems. Furthermore, some population synthesis simulations have predicted that NSBH mergers may usually contain more massive NSs compared to those in BNS mergers (e.g., Giacobbo & Mapelli 2018; Broekgaarden et al. 2021). However, due to the currently limited observations of NSBH GWs, which could hardly represent a complete mass distribution of the NSs in NSBH mergers, we use the observed mass distribution of the NSs in Galactic BNS systems as the prior of the NS mass. Here, an EoS of AP4 (Akmal & Pandharipande 1997) is adopted, since it is one of the most likely EoSs constrained by GW170817 (Abbott et al. 2018). In our calculations, we set η_X = 0.01 and f_b = 0.01, while η_γ is assumed to be a fitting parameter.

In theory, NSBH kilonovae were thought to be optically dim, but infrared-bright, compared with BNS kilonovae (e.g., Kasen et al. 2017; Kawaguchi et al. 2020; Zhu et al. 2020), because an NSBH merger probably produces a large number of lanthanide-rich dynamical ejecta with opacity κ_d ∼ 10–100 cm² g⁻¹. The disk wind ejecta is hardly likely to be lanthanide-poor, due to the lack of shock heating and neutrino irradiation during or shortly after the merger (e.g., Fernández et al. 2015; Just et al. 2015). However, Fujibayashi et al. (2020) and Kyutoku et al. (2020) recently found that the wind ejecta can still be lanthanide-poor if the viscous coefficient is not extremely high. Given the uncertainty of the wind ejecta, we adopt a wide prior distribution for the wind ejecta opacity, which is κ_w ∼ 0.5–5 cm² g⁻¹.

There are 17 free parameters, summarized in Table 2. The Markov Chain Monte Carlo method with the emcee package (Foreman-Mackey et al. 2013) is adopted to simultaneously fit the data of the gamma-ray/X-ray EEs, afterglow, and kilonova emissions of GRB 211211A. We summarize the 17 free parameters and their fitting results with 1σ credible intervals in Table 2, while the posteriors of these fitting parameters are shown in Figure 4. Except for κ_w and κ_d, these fitting parameters are convergent. The lower value of κ_w indicates that the wind ejecta of NSBH mergers could be lanthanide-poor. The best-fit light curves of the gamma-ray/X-ray EEs, afterglow, and kilonova emissions are shown in Figure 5.

Zoom In Zoom Out Reset image size

Table 2. Priors and Results for the Fitting Parameters

Parameter	Prior	Min	Max	Result
MBH/M⊙	Flat	2.22	15
χBH	Flat	0	0.997
MNS/M⊙	Gaussian	1.0	2.22
κw/cm2 g−1	Flat	0.5	5
κd/cm2 g−1	Flat	10	100
	Log-flat	10−10	1
t0/s	Log-flat	10	12.7
tp/s	Log-flat	12.7	40
tb/s	Log-flat	100	1000
s	Flat	0	15
ηγ	Log-flat	10−2	1
E0/erg	Log-flat	1050	1053
θc/rad	Flat	0	0.2
θv/rad	Flat	0	0.5
p	Flat	2	3
n/g cm−3	Log-flat	10−6	2
B	Log-flat	10−5	1

Note. The prior of M_NS is adopted as a Gaussian distribution, i.e., , consistent with the observed mass distributions of NSs in Galactic BNS systems (Lattimer 2012).

Download table as:  ASCIITypeset image

Our fitting results reveal that GRB 211211A could be a merger between a NS and a BH, with an aligned spin of . The merger would produce ∼ 0.005–0.03 M_⊙ lanthanide-poor wind ejecta and ∼0.015–0.025 M_⊙ lanthanide-rich dynamical ejecta. In Section 2.2, we suspected that GRB 211211A could have originated from an NS-first-born NSBH merger, based on the estimations of the event-rate density. Hu et al. (2022) found that for an NS-first-born NSBH binary system, the companion helium star would be efficiently tidally spun up by the NS, and would thus finally form a fast-spinning BH whose aligned spin would always be χ_BH ≳ 0.8. Thus, the mass and spin of the BH component for GRB 211211A are essentially consistent with those of NS-first-born NSBH mergers, as predicted by Hu et al. (2022).

Zoom In Zoom Out Reset image size

Our fitting results show that the fallback accretion for the interpretation of the EE starts at t₀ = 11.2 s, peaks at t_p = 30.2 s, and breaks at around t_b = 234 s in the rest frame. The peak fallback accretion rate is . Using Equation (4), one can estimate the fallback mass as . The start time and peak accretion rate are consistent with the predictions of Desai et al. (2019). Furthermore, Desai et al. (2019) predicted that rebrightening EEs, caused by the fallback accretion of r-process heating materials, may only occur if the remnant BH has a mass of ≳ 6–8 M_⊙. Based on the fitting results, GRB 211211A would finally form a ∼10M_⊙ BH after the NSBH merger, so our interpretations for the origin of GRB 211211A are self-consistent.

There are some uncertainties in our fitting results. Our results reveal that the mass ratio between the BH and the NS is ∼7 for GRB 211211A. Due to the presently limited simulations for NS-first-born NSBH mergers, it is not known if NSBH systems with such high mass ratios are common in the universe. More detailed simulations of the populations of NS-first-born NSBH mergers, based on population synthesis and detailed binary evolution, are suggested for the future. We also find that a ∼0.005−0.03 M_⊙ lanthanide-poor wind ejecta is needed in order to explain the observations of the kilonova emission associated with GRB 211211A. A large amount of lanthanide-poor wind ejecta are not typically expected to be produced after NSBH mergers (e.g., Fernández et al. 2015; Just et al. 2015), although Fujibayashi et al. (2020) and Kyutoku et al. (2020) have found that NSBH mergers can still lead to considerable lanthanide-poor wind ejecta, if the viscous coefficient is not extremely high. Future multimessenger observations between NSBH GWs and their associated kilonova emissions will help us to constrain the mass fraction of r-process elements for the wind ejecta from NSBH mergers.