Black Hole - Neutron Star mergers: using kilonovae to constrain the equation of state

The merging of a binary system involving two neutron stars (NSs), or a black hole (BH) and a NS, often results in the emission of an electromagnetic (EM) transient. One component of this EM transient is the epic explosion known as a kilonova (KN). The characteristics of the KN emission can be used to probe the equation of state (EoS) of NS matter responsible for its formation. We predict KN light curves from computationally simulated BH-NS mergers, by using the 3D radiative transfer code \texttt{POSSIS}. We investigate two EoSs spanning most of the allowed range of the mass-radius diagram. We also consider a soft EoS compatible with the observational data within the so-called 2-families scenario in which hadronic stars coexist with strange stars. Computed results show that the 2-families scenario, characterized by a soft EoS, should not produce a KN unless the mass of the binary components are small ($M_{\rm BH} \leq 6M_{\odot}$, $M_{\rm NS} \leq 1.4M_{\odot}$) and the BH is rapidly spinning ($\chi_{\rm BH} \geq 0.3$). In contrast, a strong KN signal potentially observable from future surveys (e.g. VRO/LSST) is produced in the 1-family scenario for a wider region of the parameter space, and even for non-rotating BHs ($\chi_{\rm BH} = 0$) when $M_{\rm BH} = 4M_{\odot}$ and $M_{\rm NS} = 1.2M_{\odot}$. We also provide a fit that allows for the calculation of the unbound mass from the observed KN magnitude, without running timely and costly radiative transfer simulations. Findings presented in this paper will be used to interpret light curves anticipated during the fourth observing run (O4), of the advanced LIGO, advanced Virgo and KAGRA interferometers and thus to constrain the EoS of NS matter.


INTRODUCTION
Neutron stars (NSs) comprise the densest stable matter of any observable object known to date.They lie on the threshold between ultra-dense stellar matter and a black hole (BH).Uncovering the structure and nature of matter at this threshold by establishing a single unified equation of state (EoS) has proven challenging; it remains one of the great enigmas of modern astronomy.The 2017 discovery of a transient event identified as the electromagnetic (EM) counterpart of the gravitational wave (GW) produced by the double NS merger GW170817 marked year zero of the multi-messenger GW era (Abbott et al. 2017); it divulged the potential of such mergers to unveil a richer picture of the structure and nature of matter inside compact objects such as NSs.
The EM transient comprises several observable components: gamma ray bursts (GRBs) powered by relativistic jets, the associated GRB afterglow from the interaction of the jet with the interstellar medium, and KNe.KNe are produced when matter ejected by double NS and BH-NS mergers (i.e.merger ejecta) undergoes rapid neutron capture (-process) nucleosynthesis, producing unsta-★ E-mail: lwpm20@outlook.comble neutron-rich elements (Metzger 2019).The rapid beta decay of these heavy elements powers the epic explosion known as a KN.Of the aforementioned EM components, KNe's approximately isotropic nature and relatively long emission timescales (days) make them the most attractive of the available candidates for observation of GWdetected mergers.Encoded within KN observables (light curves, spectra and polarisation) detectable from Earth is information on the physical conditions under which they were produced.EM transients and their associated KN will thus be key to placing further constraints on the EoS of NS matter.
As alluded to above, matter can be ejected during double NS and BH-NS mergers.The properties of this merger ejecta (e.g.mass and velocity) are dependent on both the binary properties of the merger system (masses and spins of the binary components) and the EoS of NS matter.The merger ejecta in turn determines the characteristics of the KN.Using this relation between the progenitor and the KN, the EoS of NS matter can be probed if the binary properties of the merger system are known (these can be identified from the GW preceding the KN).
The 1-family versus 2-families scenario refers to two theories on how compact stars exist in nature.The 1-family scenario assumes the existence of only one family of compact stars.These compact Radius (km)
stars could be NSs (made only of nucleons), hadronic stars (HSs, in which hyperons and/or various resonances appear above a threshold density) or hybrid stars (in which deconfined quarks are present above a certain density); and which it is, is currently unknown.The 1-family scenario assumes that matter inside compact stars is described by a unique EoS, and that the ground state of matter is hadronic, namely 56 Fe.If the 1-family scenario holds true in nature, the EoS must be relatively stiff in order to satisfy observational constraints imposed by analysis of the pulsar PSR J0740+6620 studied by NICER (The Neutron Star Interior Compositions ExploreR, Raaĳmakers et al. (2021a), Riley et al. (2021), Miller et al. (2021)) and by the existence of massive compact stars as PSR J09562-0607 (Romani et al. 2022), see Fig. 1.A stiff EoS describes a star with a large increase in pressure for a given increase in density.This results in stars with relatively large radii and low compactness that are expected to eject more material in a BH-NS merger, due to their large tidal deformability, resulting in a brighter KN.The 2-families scenario refers to the possible coexistence of two families of compact stars in nature: HSs and strange quark stars (QSs).It describes a phenomenology where hadronic matter is metastable, whilst strange quark matter is absolutely stable; this naturally occurs when the Bodmer-Witten hypothesis is considered (Bodmer 1971;Witten 1984).
The scenario was developed to justify the possible existence of compact stars with radii smaller than those compatible with the 1family scenario (Drago et al. 2014a).Indeed, it has been shown that the radius of a 1.4  ⊙ star cannot be smaller than about (11-11.5)km in the absence of a strong quark deconfinement phase transition (Most et al. 2018;Coughlin et al. 2019;Radice & Dai 2019;Kashyap et al. 2022) .Within the 2-families scenario this limit can be violated: HSs can reach significantly smaller radii (but not very large masses), while the most massive compact stars correspond to QSs.HSs here are described by a soft EoS, due to the inexorable formation of hyperons and delta resonances at densities higher than those seen in the 1-family scenario (Drago et al. 2014b).This results in HSs of smaller radii, and thus greater compactness than in the 1-family scenario.The 2-families scenario can be discounted if it is shown that all compact stars of mass ∼ 1.4 ⊙ have radii ≳ 11.5km (Di Clemente et al. 2022b).Notice that, although QSs can reach very large masses (Bombaci et al. 2021), numerical simulations in Kluźniak & Lee (2002) suggest that they do not produce ejecta (because of the sticking properties of strange quark matter Bauswein et al. (2010)) and in turn a KN following a merger.Therefore, within the 2-families scenario, only mergers involving HSs and a BH are investigated in this paper.More extensive numerical simulations will be needed in order to confirm that no matter is ejected in a BH-QS merger.
In this paper, two EoSs from the 1-family scenario, and a soft Eos from the 2-families scenario are investigated; these EoSs were selected based on their compatibility with recent observational constraints, including those placed by NICER and GW170817, see Fig. 1.The two EoSs of the 1-family scenario span a significant fraction of the available range in the mass-radius plane, if only one family is considered.The soft EoS of the 2-families scenario is compatible with small radii as the ones obtained from many X-ray sources (not displayed in Fig. 1) and discussed e.g. in fig. 4 of Özel & Freire (2016).Interestingly, a very recent re-analysis of PSR J0030+0451 (Vinciguerra et al. 2023) concluded that object can have either a mass of 1.70 +0.18  −0.19  ⊙ and a very large radius 14.44 +0.88 −1.05 km (probably not compatible with other astrophysical data and not displayed in Fig. 1) or a mass ∼ 1.4 ⊙ and a rather small radius, compatible with the two most soft hadronic EoSs (AP3 and SFHo) scenario and also with a QS.
By looking at Fig. 1 one is tempted to conclude that all the available constraints on masses and radii can be satisfied by the QS branch, including the recent analysis of HESS J1731-347, indicating a subsolar mass and a small radius.Indeed we have shown that QSs can account for the existence of compact stars having very small (Di Clemente et al. 2022a) or very large masses (Bombaci et al. 2021).On the other hand, it is unlikely that all compact stars are QSs: it is well known that magnetar oscillations pose challenges to QSs (Watts & Reddy 2007).Also, the analysis of the energy released by the SN1987A indicates a binding energy perfectly compatible with that of a NS (Pagliaroli et al. 2009).Instead, in order to satisfy the limit on the mass of HESS J1731-347, one needs to use the significantly larger value of the binding energy of a QS (Di Clemente et al. 2022a).
In this work, KN light curves are predicted from BH-NS merger simulations to aid in constraining the EoS of compact star matter.A similar analysis was carried out in Di Clemente et al. (2022b); we build on this analysis by improving the modelization of the KN using the 3D radiative transfer code POSSIS (Bulla 2019(Bulla , 2023)).

METHODS
In this section, we present the method steps taken to predict KN light curves for various BH-NS merger systems and EoSs.Three EoSs were investigated; two stiff EoSs (DD2, Hempel & Schaffner-Bielich (2010), and AP3, Akmal et al. (1998), both hadronic) from the 1family scenario, and one soft EoS (SFHo+HΔ, Drago et al. (2014b), also hadronic) from the 2-families scenario.KN light curves were predicted from these EoSs for various BH-NS merger systems by implementing the methods described in this section; a diagram of these methods steps can be found in Fig. 2.

Calculating the compactness and tidal deformability from an EoS
The Tolman-Oppenheimer-Volkoff (TOV) equation (Oppenheimer & Volkoff 1939) was solved using the chosen EoSs to calculate the mass ( NS ) and radius ( NS ) of a NS for a given range of central energy densities ( c );  NS and  NS were calculated for 150 linearly spaced   values between 1.4 × 10 14 and 2.5 × 10 15 g cm −3 .The  c values were chosen so that the calculations for  NS and  NS covered their full range of possible values within current observational constraints (Nättilä et al. 2017;Abbott et al. 2018;Riley et al. 2019;Miller et al. 2021;Doroshenko et al. 2022).We note that the EoSs investigated are only marginally compatible (90% confidence interval) with Doroshenko et al. (2022); see Di Clemente et al. (2022a) for other possible astrophysical paths to form such a compact object that better satisfies these constraints.Calculating  NS and  NS across all parameter space allowed for accurate interpolation of the  NS - NS relation of NSs, and thus compactness,  NS = NS / NS , across all parameter space.Fig. 1 shows the  NS - NS for each EoS, together with current observational constraints.The tidal deformability, Λ NS , was calculated using equation (1) below and those in Hinderer (2008) using results obtained from the TOV equation for each of the EoSs discussed above.Λ NS is a dimensionless quantity defined as where  2 is the Love number (Hinderer 2008);  2 is another dimensionless quantity that describes the rigidity of a star.

Calculating merger ejecta properties
During a merger, matter is ejected if the NS is tidally disrupted.This occurs when the separation of the binary components ( tidal ), at which tidal forces become strong enough to disrupt the star, is where  NS ,  BH ,  NS are the NS mass, BH mass and NS radius respectively, and  ISCO is solely a function of black hole spin ( BH ) and is defined in equations (5-7) (Foucart 2012).In cases where this condition is met, merger ejecta produced by BH-NS binary mergers can broadly be categorised into two components: dynamical ejecta and those ejected from the bound disk (Nakar 2020).The calculation of some key properties of these components, which directly determine the brightness of the resultant KN, is discussed below.
In our analysis, 36 different BH-NS merger systems were investigated, each with a unique combination of binary properties, i.e.  NS ,  BH and  BH ; details of their values can be found in Table 1.Merger ejecta properties were calculated for every merger system for each EoS in turn, resulting in 108 calculations in total (36 merger systems ×3 EoS).Calculations of the ejecta properties involved the use of fitting formulae from the semi-analytical models (Foucart et al. 2018;Barbieri et al. 2020;Kawaguchi et al. 2016).As well as being functions of  NS ,  BH and  BH , these fitting formulae also require  NS and Λ NS , which are computed as described in Section 2.1.The values of the binary properties were chosen based on observational constraints (Abbott et al. 2018;Raaĳmakers et al. 2021a;Abbott et al. 2021) and previous research on BH-NS mergers in the 1-family versus 2-families scenario, to ensure that many of the combinations would produce merger ejecta and thus a KN (Di Clemente et al. 2022b).
The ejecta properties calculated from the fitting formulae were the mass of the dynamical ejecta ( dyn ), wind ejecta ( wind ) and the average velocity of the dynamical ejecta ( dyn ); this ejecta, all gravitationally unbound from the merger remnant, directly determines the brightness of the resultant KN.To determine these properties, several other quantities must first be determined.The total ejected mass,  tot , the sum of  dyn and the mass of the bound disk ( disk ), is defined as where , ,  and  are fitting parameters defined in Foucart et al. (2018), and  = (15Λ NS ) −1/5 is a function of the NSs Λ NS .The symmetric mass ratio is a function of the binary mass component  =  NS / BH .RISCO =  ISCO / NS is the dimensionless normalised radius of the innermost stable circular orbit, and is defined as where and The total baryonic mass of the NS,   NS , is the mass contained within baryonic matter, and is defined as Equations (3-8), in addition to equations (9-13) were solved to calculate  dyn ,  wind and  dyn . dyn was calculated from where  1 ,  2 ,  3 ,  4 ,  1 and  2 are fitting parameters defined in Kawaguchi et al. (2016).In equations ( 3) and ( 9), RISCO was calculated using the parallel spin component  BH, ∥ =  BH cos(l tilt ).Here, l tilt is the angle between  BH and the total angular momentum; it was assumed to be 0 for simplicity.However, Foucart (2020) suggests that using  BH in place of  BH, ∥ in the fitting formulae provides equally accurate predictions.Numerical simulations of BH-NS mergers in Foucart et al. (2019) suggest that  dyn cannot exceed where  1 is a fraction that varies based on the mass ratio of the binary components;  1 = 0.35 was assumed based on the values of  NS and  BH used in this work.We note that the fitting formulae from Kawaguchi et al. ( 2016) and Foucart et al. (2018) are valid within a range of  NS values that encompasses those obtained with the SFHo+HΔ EoS in the two family scenario 1 .It is important to note that microphysics, e.g. the 1 Below we find that material is ejected for one system with a NS with mass of 1.4  ⊙ , a value that is higher than 1.35  ⊙ , the maximum value allowed by the range of  NS in Foucart et al. (2018).However, we verified that results are only slightly affected by adopting 1.35 instead of 1.4  ⊙ , i.e. no system produces a kilonova except the same one, with ejecta masses and kilonova light curves that are similar (Δ ∼ 0.01  ⊙ and Δmag ≲ 0.5) .
presence of hyperons and delta resonances in the EoS SFHo+HΔ, is crucial in determining the macroscopic parameters describing a compact star, as Λ NS ,  NS and also the threshold mass in a NS-NS merger.Ejecta properties, from numerical simulations of mergers, have been shown, at least in the case of NS-NS systems, to depend only on these three parameters (Hotokezaka et al. 2011;De Pietri et al. 2019).Once they are given, the amount of mass dynamically ejected and the mass of the disk can be estimated from the fitting formulae.On the other hand, microphysics can still play a crucial role when transport properties as viscosity and thermal conductivity are tested in the numerical simulations, a problem which has not yet been systematically investigated.
disk was calculated from the difference between  tot and  dyn , Since KNe originate from material that is gravitationally unbound from the merger remnant, the bound  disk does not directly contribute.However, a fraction of this disk can be ejected as a wind via different processes (Nakar 2020).The unbound wind mass was calculated from where  2 is a fraction between 0.2 and 0.4 (Nakar 2020).There is uncertainty surrounding the true value of  2 due to multiple complex mechanisms simultaneously determining the evolution of disk ejecta;  2 = 0.4 was used following Kawaguchi et al. (2020), where BH-NS merger simulations similar to those undertaken here (see section 2.3) were executed.
The average velocity of the dynamical mass is a function of the mass ratio, and was calculated from where  is the speed of light.Of the initial 36 merger systems, 23 produced ejecta and thus a KN for at least 1 of the 3 EoSs in question; stiffer EoSs (from the 1-family scenario) were more likely to produce a KN than a softer EoS (from the 2-families scenario).Across all 3 EoSs, 37 of the 108 possible combinations of binary properties and EoS (4 NS ×3 BH ×3 BH ×3EoS) produced ejecta; 20 for DD2, 13 for AP3, and 4 for SFHo+HΔ.Calculated values of  dyn ,  wind and  dyn for each of these 37 merger systems can be found in Table 2.

Ejecta model
Calculated ejecta properties ( dyn ,  wind and  dyn ) were used as input to the 3D radiative transfer (RT) code POSSIS (POlarization Spectral Synthesis In Supernovae, (Bulla 2019(Bulla , 2023))).POSSIS models the evolution of EM radiation through space; from its creation at an event such as a merger, its interaction with matter via absorption and scattering processes, to an observer on Earth.It was used in this work to predict KN light curves.The latest version of POSSIS was used (Bulla 2023), using 10 6 Monte Carlo quanta (each quanta representing a photon packet) and grid adaptions constructed to simulate the aftermath of a BH-NS merger.The grid describes the composition and architecture of ejecta properties following a merger, setting the stage for evolution of the Monte Carlo quanta that mimic EM radiation produced within the ejecta.A 3D cartesian grid is used, with a resolution of 100 cells per dimension, i.e. a total of 100 × 100 × 100 = 10 6 cells.The ejecta morphology is constructed following the BH-NS grid from Kawaguchi et al. ( 2020), with density for the wind and dynamical ejecta respectively.Here, ,  and Θ() are the radius, time and angular dependence respectively.Whilst spherical symmetry is assumed for the wind ejecta, the dynamical ejecta are largely confined to the merger plane, with an angular dependence of the density distribution described by where  is the angle from the polar axis.The assumed density distribution can be seen in Fig. 3 for one model in the grid.Ejecta are assumed to undergo homologous expansion throughout: a constant velocity, , for each cell within the grid, with a positive linear relation between  and .Hence ejecta at larger radii travel at higher velocities.The system is thus described in velocity space, as it better represents the expansion of the system (see Fig. 3).The wind ejecta are restricted to low velocities and within a range that is fixed for all simulations ( wind,min = 0.025 and  wind,max = 0.1), while the dynamical ejecta reach higher velocities ( dyn,min = 0.1) with a different extension ( dyn,max ) for each simulation depending on the value of  dyn , equations (13 and 14).
Values calculated for  dyn ,  wind and  dyn were used as input to the grid, alongside fixed values for  wind,min and  wind,max , equation ( 14).Fixed values for the electron fractions of the dynamical ejecta ( e,dyn ) and wind ejecta ( e,wind ) were also used as input:  e,dyn was assumed to be lanthanide rich material and was assigned a value of 0.1, whilst  e,wind was assigned a flat distribution between 0.2 and 0.3.All fixed values were taken from the BH-NS grid in Kawaguchi et al. (2020).

Radiative transfer simulations
Simulations were run in POSSIS for each of the 37 combinations of binary properties and EoS that produced ejecta.Outputs from POSSIS contained information on the flux of the KN, which was used to produce light curves.The flux was provided at 100 logarithmically spaced time intervals from 0.1 to 30 days after the merger.It was provided for 11 observation angles equally spaced in cosine between cos  obs = 0 • (face-on/polar) to cos  obs = 90 • (edge-on/merger plane), for a fixed azimuthal angle phi in the merger plane (exploiting the axial symmetry of the models).The spectral time series from 0.52 to 30 days were used to construct KN light curves in the broad-band filters u, g, r and z; these filters encompass the wavelength range of radiation expelled in KNe from early to late times (Metzger 2019).

RESULTS AND DISCUSSION
The absolute magnitudes of predicted KNe presented in Fig. 4, Fig. 5, Fig. 6 are overlayed on observational limits of the Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST) (Ivezić et al. 2019), as reported in Chase et al. (2022).The limits are specific to each broad-band filter.Limits for 100 and 200 Mpc are shown as this is the distance range where the majority of mergers involving a NS were detected in O3 (Kasliwal et al. 2020).First light from LSST is expected in mid-2024 and its operation to overlap with O5.

KN light curves
Predicted KN light curves are analysed for EoSs in the 1-family versus 2-families scenario.3 merger systems (of the 23 that produced ejecta for at least 1 of the 3 EoSs in question) with the following combination of binary properties are discussed: all with  NS = 1.2 ⊙ ,  BH = 6 ⊙ , but varying  BH (0, 0.3 and 0.6).These values were selected for discussion as they all fall within uncertainty ranges of the binary properties identified from the GW signal of BH-NS merger GW200115 detected during O3 in 2020 (Abbott et al. 2021); no KN was detected (Anand et al. 2021).The plots are displayed from 0.5 to 30 days with a logarithmic time scale.Light curves for the remaining 20 merger systems that produced ejecta and thus a KN are shown in Appendix A.
Fig. 4 shows that a stiffer EoS increases the likelihood of ejecta being produced by the merger and therefore, produces a brighter KN; these results are consistent with findings in Barbieri et al. (2020).For this reason, the EoSs from the 1-family scenario (DD2, AP3) are more likely to produce KNe in all broad-band filters than the EoS from the 2-families scenario (SFHo+HΔ).
A relation between  BH and the KN can be deduced from Fig. 4; higher spins result in brighter KNe for a given EoS.This is consistent with findings in Di Clemente et al. (2022b), where higher spins are shown to eject more mass during a merger.It is also consistent with the condition required for the production of merger ejecta,  tidal >  ISCO (Foucart 2012); for constant values of  NS ,  BH and  NS (i.e.constant  tidal ), the condition is more easily met for higher  BH .
For all light curves displayed in Fig. 4, a relation between the brightness of the KN and the observation angle can be identified; an observer viewing the system face-on (along the polar axis) sees a brighter KN.This is due to the angular dependence of  dyn in BH-NS mergers, expressed in equation ( 14) and seen in Fig. 3.The asymmetry of the dynamical ejecta results in larger quantities of lanthanide-rich material forming about the merger plane than closer to the polar axis.This lanthanide rich material is associated with high opacities which slows the EM radiation traversing it, resulting in fainter light curves observed in the merger plane.
Once LSST becomes operational, light curves in Fig. 4 could be used to aid in identifying whether the 1-family or 2-families scenario holds true in nature, and in turn constrain the EoS of NS matter.By first identifying the binary properties of a BH-NS merger ( NS ,  BH ,  BH ), calculating the distance to the merger, and constraining the observation angle (all from the GW signal that precedes the KN), the presence or absence of a KN signal would indicate which EoS best describes the NS.Notice anyway that GW observational data from LIGO-Virgo provide a fairly accurate measurement of the chirp mass (LIGO Scientific Collaboration et al. 2015;Acernese et al. 2015), but not an equally accurate measurement of the spin and individual masses of the components of a merger.Also, when performing the data analysis, the values of the masses and spins turn out to be strongly correlated.In a previous paper we have analyzed the impact of these observational uncertainties and we have shown how they affect the probability of observing a KN signal (Di Clemente et al. 2022b).We have shown that, if the "simple" analysis not accounting of the observational uncertainties indicates a strong KN signal (when compared with the sensitivity of LSST), the "complete" analysis shows a high probability that the KN signal can indeed be observed.Similarly, if the signal is estimated to be very weak, the complete analysis suggests a very low probability of detection.In other terms, the present observational errors associated with the GW signal do not spoil the possibility of reaching a conclusion concerning the EoS if the signal is strong.

Relation between absolute magnitude of a KN and the mass of unbound material ejected from the BH-NS merger
Further analysis is carried out on the absolute magnitude of KNe in the  band ( g ):  g is plotted as a function of the gravitationally unbound ejected mass from the merger ( unbound =  dyn +  wind ). g is averaged over all observing angles.Results for all 37 combinations of merger system and EoS that produced ejecta are combined on a single plot (Fig. 5).
As expected, a positive correlation is identified between  g and  unbound (a similar correlation exists in other broad-band filters); an increase in ejecta mass leading to an increase in r-process material surrounding the remnant, resulting in brighter KNe.The correlation was found to be a tight logarithmic fit, consistent across all epochs studied between 1 and 6 days (see Fig. 5 for results between 1 and 3 days), albeit with increasing scatter with time.Curves were fitted to data at 1, 2 and 3 days post-merger; KNe will likely be too faint for detection at later times according to the LSST limits seen in Fig. 5, so no fits were produced at later times.The following logarithmic fits were found where  1d g ,  2d g ,  3d g are  g at 1, 2 and 3 days post-merger respectively, and  unbound is the total unbound ejecta mass from the merger in solar masses.
These logarithmic fits are useful to rapidly calculate  unbound from the observed magnitude   .In this way, it is possible to check the compatibility of an observed KN signal with a specific EoS without having to run a time-consuming RT simulation such as POSSIS.

Relation between absolute magnitude of a KN and EoS
We further investigate the relation between the absolute magnitude of KNe in both the  ( g ) and  bands ( r ) and the EoS describing the NS (Fig. 6); it is again shown that stiffer EoSs produce a brighter KN for a given merger system.This relation is true for all EoSs (DD2, AP3 and SFHo+HΔ) across every merger system investigated, again consistent with findings in Barbieri et al. (2020).
As with previous plots discussed, Fig. 6 could also be used to constrain the EoS of NS matter by aiding in identifying which of the two scenarios hold true in nature.By first identifying the binary properties of a BH-NS merger from a detected GW, the  g or  r of the KN would be predicted using these plots for the identified set of binary properties, for each EoS.The  g /  r measured from the KN signal following the GW would be compared to  g /  r values predicted here for each EoS; EoSs predicting vastly different values to those measured would be ruled out.This would in turn constrain the EoS of NS matter.
For instance, Fig. 6 shows that an EoS from the two families scenario (SFHo+HΔ) is much less likely to produce a KN; a KN was only produced in 4 of the merger systems investigated.These merger systems comprised of low mass BHs (≤ 6 ⊙ ), low mass NSs (≤ 1.4 ⊙ ) and high  BH (≥ 0.3) (see Table 2 for details).On the contrary, strong KN signals are seen at larger BH and NS masses and lower  BH in the 1-family scenario (DD2 and AP3).Additionally, a strong KN signal is seen in the 1-family scenario for merger systems involving a non-rotating BH ( BH = 0) when  BH = 4 ⊙ and  NS = 1.2 ⊙ ; this result is consistent with findings in Di Clemente et al. (2022b).For example, the detection of a KN by LSST, from a merger involving a non-rotating BH of mass  BH = 4 ⊙ and a NS of mass  NS = 1.2 ⊙ at a distance ≤ 100 Mpc, would indicate that a NS from the 1-family scenario was responsible for its formation (since a KN should not be present for the 2-families scenario).The constraining power can be further improved with the Nancy Grace Roman Space Telescope (Spergel et al. 2015), which compared to LSST will allow to reach deeper magnitudes in the near infrared (∼ 2.5 magnitude in the −band, see Chase et al. 2022) and extend to even longer wavelengths (Andreoni et al. 2024).

CONCLUSIONS AND LIMITATIONS
KN light curves and their absolute magnitudes were investigated for 36 unique BH-NS merger systems and 3 EoSs.Analysis showed that stiffer EoSs result in brighter KNe, consistent with simulations in Di Clemente et al. (2022b) and Barbieri et al. (2020).This provides a concrete way to test the two-families scenario, since the hadronic family in that scheme corresponds to a very soft EoS and therefore to a very weak KN signal.We have shown that the KN signal produced in the 2-families scenario is so weak that it should be possible to distinguish between the two scenarios even by using the GW data collected in O4 or O5 and the sensitivity of future observational facilities as e.g.VRO/LSST and the Nancy Grace Roman Space Telescope .A caveat is in order: we have excluded the possibility of producing a KN signal from a BH-QS merger, because of the results of the simulation in Kluźniak & Lee (2002).Notice anyway that in that simulation the BH was not rotating and we have shown that the amount of mass ejected is strongly reduced in that case.
The SFHo+HΔ (2-family) and the AP3 (1-family) EoSs are characterized by differences of ∼1 and ∼1.5 km for 1.2 and 1.4  ⊙ NSs, respectively.As shown in Fig. 7, these correspond to differences of ≳1 mag in g-band magnitudes for the systems considered in this work.These differences are large enough to be distinguished with kilonova observations in the future.However, EoSs where the differences between a 1-family and a 2-family scenario are smaller would be more challenging to tell apart with kilonovae.For instance, two EoSs with a difference in radius of ≲ 500m for a 1.2  ⊙ NS can lead to g-band light curves that differ of ≲ 0.5mag at 1 day and that would be hard to distinguish with kilonova observations.In principle, one can hope to distinguish also between various EoSs belonging to the 1-family scenario.This possibility is highly dependent on the ability to constrain binary properties (especially  BH ) from the preceding GW signal.Due to the limitations of current technologies, constraining the EoS from Fig. 4 is still difficult, unless the signal is expected to be particularly strong.With improved sensitivities of the advanced LIGO (LIGO Scientific Collaboration et al. 2015), advanced Virgo (Acernese et al. 2015) and KAGRA (Kagra Collaboration et al. 2019) interferometers in future runs and the construction of new interferometers such as LIGO-India, the ability to identify binary properties will be greatly improved, enabling much tighter constraints on the EoS of NS matter.
By combining results from all merger systems and EoSs a positive logarithmic correlation was found between the  g and  unbound ; more ejecta mass leading to an increase in r-process rich material surrounding the remnant, resulting in brighter KNe.The correlation was consistent across all epochs studied between 1 and 6 days postmerger, albeit with an increase in scatter with time.Logarithmic fits were provided at 1, 2 and 3 days and they can be used to significantly simplify the analysis of future events.
The main limitation of this work lies in the fitting formulae used to calculate merger ejecta properties, which are known to introduce biases (Raaĳmakers et al. 2021b); assumptions were made in calculations of  dyn,max and  wind .However, an average uncertainty of ∼15% on the total ejecta mass, as estimated in Foucart et al. (2018), translates into an error of only ∼0.2 mag in the g-band light curve at 1 day according to Equation 16 we derived in our work .Limited work exists in the literature that explores the properties of ejecta mass from BH-NS mergers, making accurate assumptions in the fitting formulae challenging.The value of  1 used in the calculation of  dyn max was chosen based on simulations in Foucart et al. (2019), which focus on the near equal-mass regime (where  NS ≃  BH ) of BH-NS mergers, unlike in this paper.This may have subsequently resulted in errors in predictions of KN brightness.Similarly, the value of  2 used to calculate  wind is unlikely to be precise; the multiple mechanisms simultaneously determining the evolution of the disk ejecta are currently poorly understood, making accurate predictions of  wind difficult.Future work will likely benefit from rapidly evolving interest and research in BH-NS mergers ahead of their anticipated detection in O4 and future runs.

Figure 2 .
Figure 2. Summary of method steps.Simulation steps (in green) computations performed in each of the methods sections (2.1, 2.2, 2.3 and 2.4); the inputs and outputs of these computations are shown in blue.Future steps beyond the scope of this paper are shown in pink.

Figure 3 .
Figure 3. 2D representation of the 3D density distribution grid used in POSSIS: central circle showing the  wind (0.025 − 0.1) and lobes about the merger plane showing  dyn (0.1 −  dyn,max ).This grid illustrates the density distribution from the binary property combination  NS = 1.2 ⊙ ,  BH = 6 ⊙ and  BH = 0.3 for the EoS AP3; values of the densities in each cell vary for different combinations of binary properties and EoS, but the distribution remains the same.

Figure 4 .
Figure 4. KN light curves produced by 3 merger systems (1 merger per row) are shown in the u, g, r and z broad-band filters.All merger systems possess  NS = 1.2 ⊙ ,  BH = 6 ⊙ , but varying  BH .Light curves are shown for DD2, AP3 and SFHo+HΔ.Sequential colourmaps for light curves of each EoS indicate observation angle: colours ranging from dark to light (e.g.dark orange to light orange) indicate an observation angle ranging from a polar axis (face-on view) to the merger plane (edge-on view).Vertical lines indicate 1 day post-merger, and horizontal lines indicate LSST limits; magnitudes detectable at 100 and 200 Mpc are coloured in two shades of grey.

Figure 5 .
Figure 5.  g of KNe at 1, 2 and 3 days post-merger as a function of  unbound for all merger systems and EoSs investigated.Each data point at a given epoch represents a merger system and EoS combination.Black solid lines are curves fitted to data at each epoch.Horizontal lines indicate LSST limits in the g band; magnitudes detectable at 100 and 200 Mpc are coloured in two shades of grey.Only datapoints where  g < 0 were considered when fitting the curve, since cases where  g > 0 correspond to faint and noisy light curves that decrease the accuracy of the fit at detectable magnitudes.

Figure 6 .
Figure 6.Absolute magnitude in the  and  bands (top and bottom rows respectively) of KNe at 1 day post merger as a function of EoS for various merger systems.Panels a and c show magnitudes seen in the merger plane (cos  obs = 90 • ) whilst panel b and d show those from the polar axis (cos  obs = 0 • ).The stiffness of the EoS increase along the x-axis of each plot.Solid lines join merger systems with identical binary properties for each EoS.Horizontal lines indicate relevant LSST limits in the  and  bands; magnitudes detectable at 100 and 200 Mpc are coloured in two shades of grey.

MBH = 4 MBH = 4 Figure 7 .
Figure7.Absolute g-band magnitudes at 1 day after the merger as a function of NS radius for systems with  NS = 1.2  ⊙ (left panel) and  NS = 1.4  ⊙ (right panel) and varying BH masses and spins.Curves are computed by estimating the total ejected mass  tot from Equation 3(Foucart et al. 2018) and its connection with the g-band magnitude at 1 day from Equation16.The purple and green vertical lines mark the NS radii for the SFHo+HΔ (2-family) and the AP3 (1-family) EoSs.Magnitudes detectable by VRO/LSST at 100 and 200 Mpc are coloured in two shades of grey.

Table 2 .
Calculated values of dyn ,  wind and  dyn for each of the 37 merger systems that produced ejecta.NS ,  BH ,  dyn and  wind are expressed in units of solar masses, whilst  dyn is expressed in units of the speed of light, .NS ,  BH ,  BH  dyn