Polarized kilonovae from black hole–neutron star mergers Open Access

ABSTRACT

We predict linear polarization for a radioactively powered kilonova following the merger of a black hole and a neutron star. Specifically, we perform 3D Monte Carlo radiative transfer simulations for two different models, both featuring a lanthanide-rich dynamical ejecta component from numerical-relativity simulations while only one including an additional lanthanide-free disc-wind component. We calculate polarization spectra for nine different orientations at 1.5, 2.5, and 3.5 d after the merger and in the |$0.1\!-\!2\, \mu$|m wavelength range. We find that both models are polarized at a detectable level 1.5 d after the merger while show negligible levels thereafter. The polarization spectra of the two models are significantly different. The model lacking a disc wind shows no polarization in the optical, while a signal increasing at longer wavelengths and reaching |$\sim 1\!-\!6{{\ \rm per\ cent}}$| at |$2\, \mu$|m depending on the orientation. The model with a disc-wind component, instead, features a characteristic ‘double-peak’ polarization spectrum with one peak in the optical and the other in the infrared. Polarimetric observations of future events will shed light on the debated neutron richness of the disc-wind component. The detection of optical polarization would unambiguously reveal the presence of a lanthanide-free disc-wind component, while polarization increasing from zero in the optical to a peak in the infrared would suggest a lanthanide-rich composition for the whole ejecta. Future polarimetric campaigns should prioritize observations in the first ∼48 h and in the |$0.5\!-\!2\, \mu$|m range, where polarization is strongest, but also explore shorter wavelengths/later times where no signal is expected from the kilonova and the interstellar polarization can be safely estimated.

1 INTRODUCTION

Compact object mergers involving at least one neutron star (NS) were long regarded as the most promising scenario for a simultaneous detection of gravitational waves (GWs) and light. This expectation was spectacularly confirmed on 2017 August 17, when the short gamma-ray burst GRB 170817A (Goldstein et al. 2017; Savchenko et al. 2017) and the radioactively powered ‘macronova’ or ‘kilonova’ (hereafter KN) AT2017gfo (Coulter et al. 2017) were detected in coincidence with the GW event GW170817 (Abbott et al. 2017). The worldwide campaign that followed mapped the entire electromagnetic spectrum (e.g. Alexander et al. 2017; Andreoni et al. 2017; Arcavi et al. 2017; Covino et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Evans et al. 2017; Haggard et al. 2017; Hallinan et al. 2017; Kasliwal et al. 2017; Margutti et al. 2017; Pian 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) and led to the consensus that a binary neutron star (BNS) merger was at the origin of these signals. The detection of GRB 170817A provided the smoking gun for the long-thought association between short GRBs and BNS mergers, while the discovery of AT 2017gfo, a KN powered by the radioactive decay of r-process elements synthesized during the coalescence (Li & Paczyński 1998), confirmed that BNS mergers are the prime sites for the production of heavy elements in the Universe (e.g. Kasliwal et al. 2019; Watson et al. 2019).

Following a 19 months break, the Advanced LIGO (LIGO Scientific Collaboration 2015) and Advanced Virgo (Acernese et al. 2015) interferometers began their third observing run in 2019 April. During this run, the LIGO/Virgo Collaboration released 14 real-time public alerts thought to involve at least one NS: six BNS mergers and eight black hole (BH)–NS mergers. These alerts were followed up by several teams (e.g. Gomez et al. 2019; Lundquist et al. 2019; Ackley et al. 2020; Antier et al. 2020; Coughlin et al. 2020a, b; Gompertz et al. 2020; Vieira et al. 2020; Watson et al. 2020) with the hope to detect an electromagnetic counterpart. The larger distances and poorer localizations of these GW alerts, compared to the golden case of GW170817, required commendable efforts by the various teams to optimize each follow up. However, no convincing electromagnetic counterpart has been reported at the time of writing.

KN predictions span a wide spectrum depending on the specific properties of the ejected material such as mass, velocity, and composition. These, in turn, depend on system parameters like the binary mass ratio, BH spin, NS compactness, and equation of state. While the detection of AT 2017gfo placed important constraints on the physics of BNS mergers (see Metzger 2019 and Nakar 2020 for recent reviews), the parameter space of BHNS systems is still largely unexplored observationally (but see Yang et al. 2015 and Jin et al. 2015 for a BHNS–KN interpretation of the gamma-ray burst GRB 060614). In particular, predictions from hydrodynamical merger simulations range from systems where the NS is swallowed whole by the BH to systems where material is ejected both dynamically and from a disc (e.g. Rosswog 2005; Foucart et al. 2013; Kyutoku, Ioka & Shibata 2013; Just et al. 2015; Kiuchi et al. 2015; Kyutoku et al. 2015; Fernández, Foucart & Lippuner 2020), with the former producing no electromagnetic emission and the latter a luminous KN (e.g. Barbieri et al. 2019; Kawaguchi, Shibata & Tanaka 2020; Zhu et al. 2020). While the material ejected dynamically carries the extremely low electron fraction of the original NS and is thus rich in lanthanides (Korobkin et al. 2012), the exact composition of the material ejected from the post-merger disc is still debated (Siegel & Metzger 2018; Christie et al. 2019; Fernández et al. 2019; Miller et al. 2019; Fernández et al. 2020; Fujibayashi et al. 2020).

A potential probe to answer some of these open questions is spectropolarimetry. Thomson scattering by free electrons can be an important source of opacity in KNe at early times (Tanaka et al. 2018; Banerjee et al. 2020) and linearly polarize the escaping radiation. The polarization signal received by a distant observer is sensitive to the system inclination, the geometry and composition of the ejected material, as well as the interplay of the different ejecta components. Assuming an idealized two-component model for the material ejected in BNS mergers, Bulla et al. (2019) showed that the presence of two distinct ejecta components with different geometries and compositions leads to an overall polarization signal. Specifically, they predicted a maximum polarization level of |$\sim 0.8{{\ \rm per\ cent}}$| at ∼7000 Å and 1.5 d after the merger under favourable (edge-on) inclinations of the system. The polarization signal was found to be smaller for a more face-on orientation (consistent with AT 2017gfo; e.g. Covino et al. 2017; Hajela et al. 2019; Hotokezaka et al. 2019; Lamb et al. 2019; Troja et al. 2019), at different wavelengths and later times. Focusing on the electron-scattering dominated phase in the first few hours following the merger, Matsumoto (2018) and Li & Shen (2019) predicted a polarization signal up to |$\sim 3{{\ \rm per\ cent}}$| for BNS and BHNS mergers, respectively.

Here, we follow-up on the work by Bulla et al. (2019) and carry out Monte Carlo radiative transfer calculations to predict polarization signatures of BHNS KNe at 1.5, 2.5, and 3.5 d after the merger. Polarized components coming from other emission mechanisms that were invoked in the case of GW170817 (e.g. cocoon, cocoon breakout) are not included in this analysis. We employ a numerical-relativity simulation of BHNS mergers from Kyutoku et al. (2015) and also explore the impact of a spherical post-merger disc wind on the predicted signal. We introduce our BHNS models in Section 2, while describe the set-up of our radiative transfer calculations in Section 3. We summarize our results in Section 4 and discuss them in Section 5.

2 MODELS

In this work, we study two separate BHNS models with different ejecta geometries and compositions. The first model, referred to as BHNS-Dyn, employs 3D dynamical ejecta from the numerical-relativity simulations of Kyutoku et al. (2015). Specifically, we adopt the H4-Q3a75 model, which assumes a |$4.05\, \mathrm{M}_\odot$| BH with the dimensionless spin parameter equal to 0.75, a |$1.35\, \mathrm{M}_\odot$| NS and the H4 equation of state (Lackey, Nayyar & Owen 2006, see table II of Kyutoku et al. 2015 for more details). For our radiative transfer simulations (see Section 3), we extrapolate the profile of the dynamical ejecta from 10 ms (Kyutoku et al. 2015) to 1.5–3.5 d after the merger by assuming homologous expansion. This procedure overestimates the expansion of the inner part of the ejecta, which is significantly affected by the gravitational potential of the remnant BH at 10 ms after the onset of merger. Although this simplification is crude, we find that it does not affect our results strongly since the high degree of polarization is observed for directions not affected by this assumption (see |${\bf n}_{\#{\rm d}}$| and |${\bf n}_{\#{\rm a}}$| in Section 4). Modelling the long-term evolution of the ejecta (Rosswog et al. 2014) is necessary for more accurate predictions of the KN and is left for future work. After restricting to the amount of unbound material in the model, we obtain an ejecta mass of |$M_{\rm BHNS-Dyn}=M_{\rm dyn}=0.043\, \mathrm{M}_\odot$|⁠. As shown in the left-hand panels of Fig. 1, the dynamical ejecta are concentrated about the x–y (orbital) plane in a crescent-like shape. A lanthanide-rich composition is assumed for this component (see Section 1), with specific opacities selected from Tanaka et al. (2018; see Section 3). The one-component BHNS-Dyn model may approximate the emission from massive BHNS binary mergers, for which the dynamical ejecta could dominate the disc wind (Kyutoku et al. 2015, see their fig. 11). Radiative transfer simulations assuming homologous expansion were carried out by Tanaka et al. (2014) for a similar BHNS model without a disc-wind component, although the focus was on predicting light curves and colours and no polarization predictions were made.

The second model, referred to as BHNS-DynWind and shown in the right-hand panels of Fig. 1, is a variant of the first model where a post-merger disc wind is added to the dynamical ejecta from Kyutoku et al. (2015). The wind is assumed to be spherical and lanthanide-free in composition (see Section 3). The disc mass predicted for the H4-Q3a75 model of Kyutoku et al. (2015) is |$M_{\rm disc}=0.3\, \mathrm{M}_\odot$|⁠, |$20$|–|$40{{\ \rm per\ cent}}$| of which is predicted to be ejected in the form of a disc wind (Just et al. 2015; Siegel & Metzger 2018; Fernández et al. 2019, 2020; Miller et al. 2019; Fujibayashi et al. 2020). Here, we assume that the ejected disc wind is |$30{{\ \rm per\ cent}}$| of the disc mass, i.e. |$M_{\rm wind}=0.3\, M_{\rm disc}=0.09\, \mathrm{M}_\odot$|⁠. This material is distributed assuming a broken power law for the density profile, with ρ ∝ v⁻³ for 0.05 < v < 0.3c and ρ ∝ v⁻⁷ for v > 0.3c, where c is the speed of light. This particular density profile is chosen to mimic the one in Bulla et al. (2019). A density floor is assumed for the inner regions of the model (v < 0.05c), adding an extra |$0.016\, \mathrm{M}_\odot$| of material to the wind ejecta. Combining the wind and the dynamical component, the total ejecta mass of the BHNS-DynWind model is equal to |$M_{\rm BHNS-DynWind}=0.149\, \mathrm{M}_\odot$|⁠. When combining the dynamical ejecta and disc wind, we assume each grid cell to be dominated by the denser component. This approximation reduces the mass of the dynamical ejecta in the BHNS-DynWind model only by 1.35 × 10⁻⁴M_⊙ (⁠|$\sim 0.3{{\ \rm per\ cent}}$|⁠) compared to that in the BHNS-Dyn model.

3 RADIATIVE TRANSFER SIMULATIONS

Linear polarization is calculated for the BHNS-Dyn and BHNS-DynWind models using the 3D Monte Carlo radiative transfer code possis (Bulla 2019).¹ Briefly, possis discretizes the radiation in N_ph Monte Carlo quanta and follows them as they diffuse out the ejecta and interact with matter. The linear polarization is described in terms of the Stokes vector S = (I, Q, U)² or the normalized dimensionless Stokes vector s = S/I = (1, q, u). For each Monte Carlo quantum, the normalized Stokes vector is initialized to s₀ = (1, 0, 0) and updated after each interaction with matter. Specifically, electron scattering is assumed to polarize while bound–bound transitions to depolarize the radiation (see Bulla et al. 2015 and Bulla 2019 for more details). To speed up the calculations, synthetic observables including polarization are extracted using the Event Based Technique described in Bulla et al. (2015).

In this work, we carry out radiative transfer simulations using N_ph = 10⁷ Monte Carlo quanta. Given that the models are symmetric about the x–y equatorial plane (see top panels of Fig. 1), we restrict our investigation to the Northern hemisphere (z ≥ 0) and select nine observer orientations defined by their angles Θ_obs and ϕ_obs (see left-hand panels of Fig. 1 for a definition of these angles). The first viewing angle, n_1a, is chosen along the z-axis (cos Θ_obs = 1 and ϕ_obs = 0°). Four viewing angles are then selected 60° away from the pole (cos Θ_obs = 0.5) and with different azimuthal angles: n_2a (ϕ_obs = 0°), n_2b (ϕ_obs = 90°), n_2c (ϕ_obs = 180°), and n_2d (ϕ_obs = 270°). Similarly, four viewing angles are selected in the equatorial plane (cos Θ_obs = 0): n_3a (ϕ_obs = 0°), n_3b (ϕ_obs = 90°), n_3c (ϕ_obs = 180°), and n_3d (ϕ_obs = 270°). The angles ϕ_obs and Θ_obs are chosen to sample uniformly the 2π solid angle. We note that the dynamical ejecta are distributed in a crescent-like shape centred at ϕ_obs ∼ 270° (see bottom left-hand panel of Fig. 1).

Polarization levels are computed at three individual epochs: 1.5, 2.5, and 3.5 d since merger. These epochs correspond to the first three polarimetric observations of AT 2017gfo (Covino et al. 2017) and are chosen for ease of comparison with the results in Bulla et al. (2019). Coordinates and densities for the BHNS-Dyn and BHNS-DynWind models are remapped to the desired epochs according to homologous expansion (i.e. assuming constant velocity v and density dropping as ρ ∝ t⁻³).

Two main changes are made compared to the calculations presented in Bulla et al. (2019) for BNS mergers. First, the initial location of the Monte Carlo quanta is sampled from the ejecta mass distribution (i.e. more quanta are created at higher compared to lower densities). This approach is a better description of the underlying physics, where a fraction of the γ-ray photons from the radioactive decay of r-process elements are absorbed by matter and re-emitted at longer wavelengths as a thermal component (Li & Paczyński 1998). This assumption improves on the approximation adopted in Bulla et al. (2019), where Monte Carlo quanta were emitted only from a spherical photosphere predefined at some given velocity/radius. Secondly, the full wavelength information of opacities is used to predict polarization spectra between 0.1 and |$2\, \mu$|m rather than polarization levels at individual wavelengths. In particular, realistic opacities based on atomic calculations from Tanaka et al. (2018) are used for both lanthanide-free (lf) and lanthanide-rich (lr) compositions at each epoch. At the wavelengths investigated in this study, electron scattering and bound–bound interactions are the main source of opacity. Electron-scattering and bound–bound opacities used in this work are shown in Fig. 2. Their relative contribution is sensitive to the density/temperature conditions within the ejecta, with electron scattering (bound–bound) opacities typically decreasing (increasing) with time as the ejecta expand, cool down, and the different elements recombine to lower-ionization stages. The steep decrease in bound–bound opacities with wavelength is due to the number of transitions being larger at bluer wavelengths. The values adopted for the electron scattering opacities at the different epochs are |$(\kappa ^{\rm lf}_{\rm es},\kappa ^{\rm lr}_{\rm es})_{1.5\, {\rm d}}=(0.0080,0.0068)\, {\rm cm}^2\, {\rm g^{-1}}$|⁠, |$(\kappa ^{\rm lf}_{\rm es},\kappa ^{\rm lr}_{\rm es})_{2.5\, {\rm d}}=(0.0044,0.0046)\, {\rm cm}^2\, {\rm g^{-1}}$|⁠, and |$(\kappa ^{\rm lf}_{\rm es},\kappa ^{\rm lr}_{\rm es})_{3.5\, {\rm d}}=(0.0040,0.0032)\, {\rm cm}^2\, {\rm g^{-1}}$|⁠. Bound–bound opacities |$(\kappa ^{\rm lf}_{\rm bb},\kappa ^{\rm lr}_{\rm bb})$| at each wavelength and time are computed through tenth-order polynomial fits of the opacities from Tanaka et al. (2018; see thick solid lines in Fig. 2).

4 RESULTS

In the following, we summarize the results of the simulations described in Section 3. We begin with the BHNS-Dyn model in Section 4.1 and then discuss the BHNS-DynWind model in Section 4.2. The implications of our results in terms of observations will be discussed in Section 5.

4.1 Dynamical ejecta

Polarization spectra computed for the BHNS-Dyn model are shown in Fig. 3. Viewing angles are grouped in different rows depending on the polar angle Θ_obs (polar to equator from top to bottom), while different columns show predictions at different epochs (1.5, 2.5, and 3.5 d from left to right). The polarization signal shows a clear wavelength-, angular-, and time-dependence. We address each of these three behaviours below.

We start by focusing on the polar orientation (cos Θ_obs = 1, n_1a) at 1.5 d after the merger, shown in the top-left panel of Fig. 3. The predicted level is strongly wavelength dependent, with polarization consistent with zero below |$\sim 1.2\, \mu$|m and increasing moving to longer wavelengths. The polarization level reaches |$P\sim 1{{\ \rm per\ cent}}$| at |$\sim 2\, \mu$|m, at the boundary of the simulated wavelength domain. The wavelength-dependence of the signal is a direct consequence of the adopted opacities and it is controlled by the ratio κ_es/κ_bb between polarizing electron-scattering contributions and depolarizing bound–bound interactions (see Section 3). In the BHNS-Dyn model, the ejecta are assumed to be lanthanide-rich and thus characterized by relatively high bound–bound opacities. As shown in the left-hand panel of Fig. 2 (see red lines), κ_es/κ_bb ≪ 1 below |$\sim 1.2\, \mu$|m, hence the negligible polarization levels obtained in this wavelength range. Starting from |$\sim 1.2\, \mu$|m, the κ_es/κ_bb ratio reaches values larger than 0.1 and eventually exceed unity, hence the rise in polarization moving farther into the infrared.

The signal is strongly viewing-angle dependent due to the large-scale asymmetries in the BHNS-Dyn model and the resulting projection effects. Polarization levels are found to be the highest for orientations facing the crescent-like ejecta distribution, namely ϕ_obs = 0° (n_2a and n_3a) and ϕ_obs = 270° (n_2d and n_3d), cf. with bottom-left panel of Fig. 3. As mentioned in Section 2, these viewing angles are the least affected by the assumption of homologous expansion in our calculations. At ϕ_obs = 0°, the polarization level at |$\sim 2\, \mu$|m increases from |$\sim 1{{\ \rm per\ cent}}$| at the Pole (n_1a) to |$\sim 1.5{{\ \rm per\ cent}}$| at Θ_obs = 60° (n_2a) and |$\sim 3{{\ \rm per\ cent}}$| at the Equator (n_3a). At ϕ_obs = 270°, the polarization peak at |$\sim 2\, \mu$|m reaches even higher levels: |$\sim 3{{\ \rm per\ cent}}$| at Θ_obs = 60° (n_2d) and |$\sim 6{{\ \rm per\ cent}}$| at the Equator (n_3d). In contrast, the polarization is much lower (⁠|$\lesssim 0.5{{\ \rm per\ cent}}$|⁠) for orientations at the other side of the crescent-like ejecta distribution (ϕ_obs = 90° and ϕ_obs = 180°) as they receive fewer polarizing contributions. However, the results for these directions should be taken with care since the assumption of the homologous expansion with a central cavity may introduce significant artificial effects (see Section 2). Quantifying the impact of this assumption on the polarization predicted along these orientations is left for future studies. Here, we note that lanthanide-rich material filling the cavity would be associated with high densities/opacities, with radiation created in these regions diffusing out on long time-scales and affecting only marginally the early-time polarization predicted in this study.

The time-dependence of the polarization signal is very strong. At 2.5 d after the merger, i.e. only 1 d after the first simulated epoch, the levels are negligible at all wavelengths and for all the different viewing angles (⁠|$\lesssim 0.3{{\ \rm per\ cent}}$| at |$2\, \mu$|m). Polarization signals are even weaker and essentially consistent with zero at 3.5 d after the merger. The extremely rapid drop in polarization is a consequence of the time evolution of the κ_es/κ_bb ratio and of the electron-scattering optical depth τ_es. As shown in Fig. 2, the bound–bound opacities rapidly increase with time as the ejecta cool and atoms recombine, while the variation in electron-scattering opacity is much smaller. Therefore, the κ_es/κ_bb ratio decreases extremely fast, with κ_es/κ_bb ∼ 0.05 (0.003) at |$\sim 2\, \mu$|m and 2.5 (3.5) d after the merger. In addition, the density drop due to the expansion (a factor of ∼5 between the first and the second epoch and of ∼13 between the first and the third epoch, ρ ∝ t⁻³) causes the electron-scattering optical depth τ_es to rapidly decrease with time. These effects lead to a rapid decrease in polarizing contribution and therefore resulting polarization signal.

4.2 Dynamical  + wind ejecta

Polarization spectra computed for the BHNS-DynWind model are shown in Fig. 4, with predictions for different viewing angles and epochs organized in different panels as in Fig. 3 (see Section 4.1). Similarly to what is seen in the BHNS-Dyn model, the polarization signal predicted for the BHNS-DynWind model depends on wavelength, viewing angle, and time.

Again, we begin by focusing on the polar viewing angle (cos Θ_obs = 1, n_1a) at 1.5 d after the merger, shown in the top-left panel of Fig. 4. The overall polarization level is |$\lesssim 0.5{{\ \rm per\ cent}}$| at all wavelengths investigated and shows a characteristic ‘double-peak’ structure: one peak is seen in the optical (⁠|$\sim 0.7\, \mu$|m) and one in the infrared (⁠|$\sim 1.5\, \mu$|m). This behaviour can be understood by inspecting the wavelength-dependence of opacities from Fig. 2 (see left-hand panel). While the BHNS-Dyn model is characterized by lanthanide-rich opacities, the BHNS-DynWind model features both lanthanide-free opacities from the wind and lanthanide-rich opacities from the dynamical ejecta. The interplay between the two components is responsible for shaping the polarization spectrum, which can be broadly summarized by looking at three separate wavelength ranges. In the first range, |$\lambda \lesssim 0.5\, \mu$|m, κ_es/κ_bb ≪ 1 for both components and thus the polarization signal is consistent with zero. In the second range, |$0.5\lesssim \lambda \lesssim 1.2\, \mu$|m, |$\kappa ^{\rm lf}_{\rm es}/\kappa ^{\rm lf}_{\rm bb}\gtrsim 1$| and |$\kappa ^{\rm lr}_{\rm es}/\kappa ^{\rm lr}_{\rm bb}\ll 1$| and thus the observer receives polarizing contributions from the spherical wind and depolarizing contributions from the dynamical ejecta. While all the polarizing contributions would cancel to zero if the spherical wind was the only component in the model, the presence of a depolarizing asymmetric component (i.e. the dynamical ejecta) leads to an overall polarization signal. This effect, which is maximized at |$\sim 0.7\, \mu$|m, is the same as the one first discussed by Bulla et al. (2019) (see also Section 5). In the third range, |$\lambda \gtrsim 1.2\, \mu$|m, κ_es/κ_bb ≳ 1 in both components and thus the ejecta are dominated by electron scattering. While this resulted in large polarization signals for the strongly asymmetric BHNS-Dyn model, the infrared polarization predicted in the BHNS-DynWind model is modest due to the configuration approaching a pure-electron scattering case in a spherical photosphere (Bulla et al. 2019).

The viewing-angle dependence of the polarization signal is shown in the different rows of Fig. 4. In general, the polarization spectra continue to display a double-peak profile as the one described above. In the optical, the strongest polarization is seen for orientations at ϕ_obs = 0° and ϕ_obs = 270° as it was the case in the BHNS-Dyn model. This is because the dynamical ejecta are distributed preferentially towards these orientations, thus blocking and breaking the cancellation of polarizing contribution from the wind more effectively (see above). Quantitatively, a peak of |$P\sim 1.2{{\ \rm per\ cent}}$| is found at |$\sim 0.7\, \mu$|m when the system is viewed from the n_2d and n_3d orientation. In contrast, the infrared polarization is less viewing-angle dependent reflecting the overall symmetry achieved in the polarizing contributions at these wavelengths (see above).

The BHNS-DynWind model is also characterized by a rapid depolarization of the signal, with levels below |$\sim 0.3{{\ \rm per\ cent}}$| at 2.5 and 3.5 d after the merger. As explained in Section 4.1, this behaviour is caused by both the κ_es/κ_bb ratio and the (polarizing) electron-scattering opacity decreasing with time. We note, however, that the presence of a lanthanide-free component in the BHNS-DynWind model leads to relatively higher levels compared to those in the BHNS-Dyn model, especially at wavelengths shorter than |$\sim 1.3\, \mu$|m.

5 DISCUSSION AND CONCLUSIONS

Using the Monte Carlo radiative transfer code possis, we predict polarization spectra of two 3D BHNS models as viewed from nine different orientations. While both models feature the presence of a lanthanide-rich dynamical ejecta component (Kyutoku et al. 2015), the BHNS-DynWind model includes also a lanthanide-free spherical wind component that is not present in the BHNS-Dyn model.

A clear polarization signal is found in the wavelength range investigated (⁠|$0.1\!-\!2\, \mu$|m), reaching |$\sim 6{{\ \rm per\ cent}}$| for the BHNS-Dyn model viewed from an equatorial viewing angle. The signal is restricted to early epochs, |$1.5\,$| d after the merger, and disappears within a day or two. This indicates that a polarization signal could be detected in future KNe resulting from the merger of an NS and a BH. Our quantitative predictions call for spectropolarimetric observations that are rapid and cover a wavelength range as wide as possible. While prioritizing observations in the first ∼48 h and in the |$0.5\!-\!2\, \mu$|m range, where the signal is expected to be the highest, the best strategy would be to obtain additional observations at later times and shorter wavelengths. This would allow to explore epochs and wavelengths where the intrinsic signal is predicted to be negligible and thus the interstellar polarization can be robustly estimated and removed from all the observations (see also Bulla et al. 2019). In this respect, we note that the predicted wavelength dependence in the polarization of the KN is quite different from that seen in (Galactic and host-galaxy) interstellar dust, making the separation of these two components easier even for low signal-to-noise data (Patat et al. 2015).

The polarization spectra are significantly different between the two models. While the BHNS-Dyn model is polarized only in the infrared (⁠|$\lambda \gtrsim 1.2\, \mu$|m), the BHNS-DynWind model features a ‘double-peak’ polarization spectrum with a peak in the optical (⁠|$\sim 0.7\, \mu$|m) and another in the infrared (⁠|$\sim 1.3\!-\!2\, \mu$|m depending on the viewing angle). Even in cases where the two peaks are difficult to resolve in data (e.g. because of low amplitudes), their distinction should be facilitated by the corresponding polarization angles (i.e. q/u) changing differently with time and wavelength (see e.g. bottom panels of Fig. 4). The clear difference between the two models is ascribed to the presence of a lanthanide-free spherical wind component, indicating that polarimetric data of future events might constrain the neutron richness of the disc-wind component in BHNS mergers, a property that is yet not settled in the literature (Siegel & Metzger 2018; Christie et al. 2019; Fernández et al. 2019, 2020; Miller et al. 2019; Fujibayashi et al. 2020). Specifically, the detection of a polarization signal at optical wavelengths would be a smoking gun for the presence of a lanthanide-free disc-wind component, especially if accompanied by an early-blue emission as in AT 2017gfo. Alternatively, a non-detection in the optical together with a strong (⁠|$\gtrsim 2{{\ \rm per\ cent}}$| at |$\sim 2\, \mu$|m) signal in the infrared would point to a lanthanide-rich composition for the entire ejecta. Knowledge about the composition can constrain the disc mass and the launching mechanism/time-scales of the post-merger ejecta (Fujibayashi et al. 2020).

Our predictions call for (spectro)polarimetric observations extending as much as possible into the infrared domain. While infrared polarimeters are rare on large telescopes, Tinyanont et al. (2019) showed that high signal-to-noise polarimetry of bright transients should be feasible at these wavelengths with the WIRC + Pol instrument mounted on the 200 inch (∼5 m) Hale Telescope at Palomar Observatory. As shown by Tinyanont et al. (submitted), the WIRC + Pol instrument can observe sources as faint as 14.5 mag in the J band to a polarimetric accuracy of d|$P\sim 0.1{{\ \rm per\ cent}}$| per spectral channel in less than 2 h. This accuracy would be sufficient to detect a polarization signal at ∼1 d after the merger in a nearby (∼10 Mpc) AT2017gfo-like KN. An even more promising option is offered by the IRCS instrument at the 8.2-m Subaru telescope (Watanabe et al. 2018). A nominal³ accuracy of d|$P\sim 0.5\, {{\ \rm per\ cent}}$| can be achieved with 1-h exposure for J ∼ 20.7 mag, H ∼ 20.0 mag, and K ∼ 19.5 mag, which would enable the detection of high polarization signals as those predicted here in a AT2017gfo-like KN up to ∼150–200 Mpc. Additional near-infrared polarimeters suitable for bright transients include LIRIS (Long-slit Intermediate Resolution Infrared Spectrograph) on the 4.2-m William Herschel Telescope (see e.g. Manchado et al. 2004; Wiersema et al. 2012) and SofI on the 3.6-m New Technology Telescope (see e.g. Higgins et al. 2019).

Within a given model, a clear viewing-angle dependence is found in the polarization signal. This brings the potential to pinpoint the system inclination of future BHNS mergers using spectropolarimetric observations, especially when combined with constraints from different probes (e.g. KN optical/near-infrared spectra and light curves, radio observations). In the future, systematic studies of polarization from various configurations of the ejecta components will be crucial to determine the viewing angle accurately. Estimating the inclination of NS mergers is important to understand the physics of these systems and in particular their ejecta distribution.

Our polarization predictions can be compared to those in the literature. Matsumoto (2018) and Li & Shen (2019) predicted polarization arising from a fast ejecta component (≳0.5c) in BNS and BHNS, respectively, where free neutrons survive and their β-decay electrons dominate the scattering opacity in the first few hours after the merger. This high-velocity component is neglected in this work and thus a one-to-one comparison with Matsumoto (2018) and Li & Shen (2019) is not possible. A more meaningful comparison can be done with Bulla et al. (2019), which focused on BNS mergers and assumed a two-component (wind + dynamical) ejecta geometry similar to that in the BHNS-DynWind model presented here. However, the dynamical ejecta in their model have a torus-like shape while ours a crescent-like shape covering only ∼ half of the azimuthal angle. The assumed axial symmetry in the Bulla et al. (2019) model leads to a polarization signal that is zero at the pole, when the ejecta are symmetric in projection, and increases monotonically towards the equator (see their fig. 2). Observationally, this means that a BNS KN would be unpolarized or very lowly polarized when viewed closed to face-on, as GW170817/AT2017gfo was (e.g. Covino et al. 2017; Hajela et al. 2019; Hotokezaka et al. 2019; Lamb et al. 2019; Troja et al. 2019). In contrast, our BHNS-DynWind model is not axially symmetric and thus it is polarized in any orientation investigated. This suggests that BHNS systems might be easier to detect in polarization compared to BNS, although we note that this statement depends on whether dynamical ejecta deviating substantially from the axial symmetry can be produced in BNS mergers and/or whether a lanthanide-free component is present in BHNS mergers. Finally, we note that predictions from Bulla et al. (2019) for an equatorial viewing angle match remarkably well the values obtained by BHNS-DynWind model when viewed from the n_3a orientation (see yellow stars in Fig. 4). This is not too surprising since the torus-like geometry in the Bulla et al. (2019) model is relatively similar in projection to the crescent-like shape in our BHNS-DynWind model when viewed from n_3a.

ACKNOWLEDGEMENTS

The authors are thankful to the anonymous reviewer for helpful comments that improved this paper. Koutarou Kyutoku was supported by Japanese Society for the Promotion of Science (JSPS) Kakenhi Grant-in-Aid for Scientific Research (No. JP16H06342, No. JP17H01131, No. JP18H05236, No. JP19K14720, No. JP20H00158).

DATA AVAILABILITY

The polarization spectra presented in this work are publicly available at https://github.com/mbulla/kilonova_models.

Footnotes

REFERENCES