Population properties of neutron stars in the coalescing compact binaries

We perform a hierarchical Bayesian inference to investigate the population properties of the coalescing compact binaries involving at least one neutron star (NS). With the current observation data, we can not rule out either of the Double Gaussian, Single Gaussian and Uniform NS mass distribution models, although the mass distribution of the Galactic NSs is slightly preferred by the gravitational wave (GW) observations. The mass distribution of black holes (BHs) in the neutron star-black hole (NSBH) population is found to be similar to that for the Galactic X-ray binaries. Additionally, the ratio of the merger rate densities between NSBHs and BNSs is estimated to be about 3 : 7. The spin properties of the binaries, though constrained relatively poor, play nontrivial role in reconstructing the mass distribution of NSs and BHs. We find that a perfectly aligned spin distribution can be ruled out, while a purely isotropic distribution of spin orientation is still allowed.


INTRODUCTION
Recently, the observations of gravitational wave (GW) signals from two neutron star-black hole (NSBH) coalescences, i.e., GW200105 and GW200115, have been reported by LIGO/Virgo/KAGRA Collaboration (LVKC; Abbott et al. 2021a). These two events are also the first confident observations of NSBH binaries in the Universe, although previously the modeling of the kilonova emission of hybrid GRB 060614 was also strongly in favor of an NSBH merger origin (Yang et al. 2015;Jin et al. 2015). Since the first successful detection of GW (Abbott et al. 2016) in 2015, about ∼ 50 compact binary coalescences (CBCs) have been formally reported (Abbott et al. 2021b), including ∼ 45 events from binary black hole (BBH) coalescences and several mergers involving at least one neutron star (e.g., Abbott et al. 2017Abbott et al. , 2020a. Very recently, The LIGO Scientific Collaboration et al. (2021) releases a deep extended catalog, which reports another 8 new events. With a rapidly increasing sample of GW events, our understanding of the population properties of the stellar BHs in the Universe has been advanced (e.g., Talbot & Thrane 2017;Abbott et al. 2019a;Fishbach & Holz 2020;Abbott et al. 2021c;Wang et al. 2021b;Li et al. 2021;Kimball et al. 2021;Tiwari & Fairhurst 2021;Kimball et al. 2021;Galaudage et al. 2021b), and some formation/evolution processes of the compact binaries are being revealed (e.g., Kimball et al. 2020Kimball et al. , 2021Safarzadeh & Wysocki 2021;Baxter et al. 2021;Tang et al. 2021b;Mapelli et al. 2021;Wang et al. 2021a). However, due to the limited observation of NSBH and BNS mergers up to The corresponding author: yzfan@pmo.ac.cn (Y.Z.F) arXiv:2108.06986v2 [astro-ph.HE] 28 Oct 2021 now, the population properties of these kinds of coalescence systems are hard to probe. Moreover, Tang et al. (2020) found that the misclassification of BBH into NSBH makes it more challenging to reconstruct the mass function of the BHs.
Thanks to the observation of Galactic radio pulsars, a large number of NS masses have been measured (Valentim et al. 2011;Antoniadis et al. 2016;Özel & Freire 2016;Alsing et al. 2018;Farrow et al. 2019;Rocha et al. 2019;Shao et al. 2020a;Galaudage et al. 2021a;. The mass distribution of these Galactic NSs has been investigated and a double gaussian model with two peaks at ∼ 1.35M and ∼ 1.9M are found to be able to well reproduce the data (Alsing et al. 2018;Shao et al. 2020a). Nevertheless, it is very interesting to investigate the mass distribution of NSs in the Universe, and the GW observations provide the unprecedented valuable chance. Recently, Landry & Read (2021) tried to reconstruct the mass distribution of this extragalactic population of NSs and found that the mass function seems to be more consistent with a uniform distribution than the bimodal Galactic population. In their hierarchical inferences, only the masses of compact objects have been taken into account, and the spin information of the compact binaries has not been included. However, as found in the literatures (e.g. Baird et al. 2013;Creswell et al. 2018;Pratten et al. 2020), for the GW events involving at least one BH, there is mass-spin degeneracy. Additionally, the measurements of both the component masses and misaligned spins are of prime importance in understanding the origin and evolution of astrophysical compact binaries (Rodriguez et al. 2016;Qin et al. 2018Qin et al. , 2019. Binaries born in isolated evolution are expected to form with nearly aligned spins (Kalogera 2000), and after considering the supernova kicks, the compact object binaries tend to still have small misalignments (Rodriguez et al. 2016;Gerosa et al. 2018;Wysocki et al. 2018). On the contrary, binaries originating from dynamical capture are expected to have isotropic distribution for their spin orientations (Vitale et al. 2017;Rodriguez et al. 2016). By surveying the spin properties of the BBH population, Abbott et al. (2021c) find that the current BBH merger events likely originate from both formation channels (see also Wong et al. 2021;Zevin et al. 2021).
In this work, we perform a hierarchical Bayesian inference to investigate the population properties of NSs detected by LIGO and Virgo with the information of their masses and spins. In Sec. 2, we introduce the data and the models used for inference, and in Sec. 3 we present the results. In Sec. 4, we carry out simulations with a mock population and evaluate the feasibility for reconstruction of the NS mass distribution with one hundred of BNSs/NSBHs detected in the design sensitivity run of Advanced LIGO/Virgo, and Sec. 5 is our conclusion and discussion.

Selected events
Our hierarchical analysis is based on the GW events involving at least one NS, which consist of four published confident events, including GW170817, GW190425, GW200105 and GW200115 (Abbott et al. 2017(Abbott et al. , 2020a(Abbott et al. , 2021a, and one marginal event, GW190426 152155. Though the nature of GW190425 and GW190426 152155 are not fully solved yet Li et al. 2020), the former (latter) is more likely a BNS (NSBH) (Abbott et al. 2020a(Abbott et al. , 2021b. We do not include GW190814 in our analysis, because the NS nature of the secondary object in GW190814 is inconsistent with either the maximum mass of nonrotating NS (Nathanail et al. 2021;Tang et al. 2021a), determined by the multi-messenger analyses of GW170817/GRB 170817A/AT2017gfo (Ruiz et al. 2018;Rezzolla et al. 2018;Shibata et al. 2019;Shao et al. 2020b;Fan et al. 2020) or the constraints obtained from energetic heavy-ion collisions (Fattoyev et al. 2020). Though a rapidly rotating NS may be allowed to have the mass of the secondary object in GW190814, it will coalesce to a BH before merger due to the rotational instabilities (Shao et al. 2020b;Biswas et al. 2021). Therefore, our data set includes two BNS events and three NSBH events (or candidate), and the parameter estimation results for each event are adopted from Abbott et al. (2019bAbbott et al. ( , 2021b 1 . Since we need the spin information (including spin magnitude and orientation) for each compact object, all the posterior samples should be estimated by the precession waveforms. So in this work we adopt the "IMRPhenomPv2NRT lowSpin posterior", "PhenomPNRT-LS", "C01:PhenomXPHM low spin", "C01:PhenomXPHM low spin", and "PrecessingSpinIMRHM" posterior samples for GW170817, GW190425, GW200105, GW200115, and GW190426 152155, respectively.
Motivated byÖzel et al. (2010) and Abbott et al. (2021c), we introduce a Truncated Gaussian and a Truncated Power Law for the mass distribution of BHs, In the above equations, m NS and m BH are the mass of NS and BH, while Φ 1 , Φ 2 , Φ, and Φ are the normalization constants, and other parameters are described in Table 1.
The parameterization for the component (BH or NS) spin magnitudes and tilts is similar to the Default model in Abbott et al. (2021c). The spin tilt of each component in a binary is assumed to be independently drawn from the same underlying distribution, see Eq. (6) below. For simplicity, we use two truncated gaussians to fit the dimensionless spin magnitudes of BHs and NSs (i.e., a BH and a NS ), respectively. Therefore, the model of spin reads π(a BH |µ BH a , σ BH a ) = N (a BH |µ BH a , σ BH a )/Φ BH a , for a BH ∈ (0, 1), π(a NS |µ NS a , σ NS a ) = N (a NS |µ NS a , σ NS a )/Φ NS a , for a NS ∈ (0, a NS max ), where z = cos θ 1,2 (θ 1,2 is the tilt angle between component spin and binary's orbital angular momentum), Φ z , Φ BH a , and Φ NS a are the normalization constants. The second term of Eq. (6) represents the probability density of systems that have an isotropic spin distribution, and ζ represents the mixing fraction of the field mergers. Note that, NSs are able to achieve larger spin tilts since they are lower mass and the orbital plane of the binary can theoretically be tilted a larger degree due to the impact of the supernova (Rodriguez et al. 2016). For simplicity, we assume the same tilt angle distribution for both BHs and NSs, since the currently limited data are not able to recognize such difference between BHs and NSs. The highest dimensionless NS spin implied by pulsar-timing observations of binaries that merge within a Hubble time is 0.04 (Stovall et al. 2018;Zhu et al. 2018), therefore, it's reasonable to set a NS max = 0.05. All the description of the parameters and the priors of each model are summarized in Table 1 Combining the mass and spin models, we get the synthesis models for NSs and BHs, where Λ m NS (Λ m BH ), Λ a NS (Λ a BH ), and Λ z are the population parameters for mass distribution of NS (BH), spin magnitude of NS (BH), and spin orientations, respectively. And Since the secondary objects of all the binaries are NSs, the model for them is π(m 2 , a 2 , z 2 |Λ 2 ) = π(m NS , a NS , z 2 |Λ NS ).
While the primary objects can be NSs or BHs, we define r BNS as the fraction of BNS populations, then the model for primary objects becomes So the synthesis models for the NSBH and BNS binaries are represented by π(θ|Λ) = π(m 2 , a 2 , z 2 , m 1 , a 1 , z 1 |Λ) = π(m 2 , a 2 , z 2 |Λ 2 )π(m 1 , a 1 , z 1 |Λ 1 ), where θ is the source parameters for each event, and Λ = Λ 1 ∪ Λ 2 = Λ NS ∪ Λ BH ∪ {r BNS } is the population parameters for synthesis models. Additionally, in order to test whether the distribution of NSs observed via GW is consistent with that observed electromagnetically in the Galaxy, we fix the parameters of Double Gaussian for NSs as m min = 1M , m max = 2.25M , µ 1 = 1.36M , σ 1 = 0.09M , µ 2 = 1.9M , σ 2 = 0.5M , and r 1 = 0.65 (i.e., the median values of the parameters taken from Shao et al. (2020a), where the NS mass distribution was obtained based on the observed systems involving at least one NS in the Galaxy). Thus, there are 4 NS mass distribution models (i.e., Double Gaussian 2 , Single Gaussian, Uniform, and Galactic distribution) × 2 BH mass distribution models (i.e., Truncated Power Law and Truncated Gaussian) × 1 spin model. In other words we have 8 synthesis models.

Hierarchical likelihood
With the data of the observed events {d} that described in Sec. 2.1, and the population models described in Sec. 2.2, we perform a hierarchical Bayesian inference following Abbott et al. (2021c) and Thrane & Talbot (2019). For the given data {d} from N det GW detections, the likelihood of Λ can be expressed as where ξ i (Λ) means the detection fraction, and the single-event likelihood L(d i |θ i ) can be estimated using the posterior samples described in Sec. 2.1 (see Abbott et al. (2021c) for detail), then the Eq. (14) takes the form of where π φ (θ k i ) is the default prior used for parameter estimation of each individual event, and n i is the size of posterior samples. For the estimation of ξ i (Λ), we use a Monte Carlo integral over detected injections as introduced in Abbott et al. (2021c) and Tiwari (2018). And the injection campaigns are made by injecting simulated signal to the LIGO-Virgo detectors with noise curves 3 . We define the threshold of the injected signal being detected as the network signal-to-noise ratio greater than 12. For the hierarchical Bayesian inference, we use the Bilby package (Ashton et al. 2019) and Dynesty sampler (Speagle 2020).

RESULTS
In this section we report the results obtained by the hierarchical inference using synthesis models described in Sec. 2. We calculate the Bayes factor for each model (M a ) relative to the model of Double Gaussian + Truncated Gaussian where Z a and Z b are the evidences for the model a and b. Formally, the correct metric to compare two models is the odds ratio O a b = (Z a /Z b ) · (π a /π b ), we assume the prior odds π a and π b for the models are equal, so that O a b = B a b . With the logarithm Bayes factors summarized in Table 2, we can conclude that all three NS mass function models are comparable to fit the data. Nevertheless, it is interesting to see that the Galactic model (i.e., the fixed Double Gaussian model) is slightly preferred than the others. It is certainly too early to conclude that the mass distribution of NSs with GW is the same as that of the Galactic NSs, and much more data are needed to draw a robust conclusion. As for BH mass distribution, we find that the Truncated Power Law model has a preference over the Truncated Gaussian by ln B =∼ 0.5 − 1.5. The inferred population parameters for NS and BH mass distribution models are presented in Table 3. The results for NS mass distribution are inferred by the Double Gaussian, Single Gaussian, and Uniform models, assuming a BH mass distribution model of Truncated Power Law, and the results for BH mass distribution are obtained using Truncated Power Law and Truncated Gaussian models, assuming a NS mass distribution model of Double Gaussian. Actually, the results assuming other BH mass distribution model/ NS mass distribution models are rather similar. In order to build more intuition, we plot the mass distributions of BHs/NSs as shown in Fig. 1. For the Double Gaussian model of NS mass distribution, it is clear that there is a peak around 1.35M , which is consistent with the first peak found in the mass distribution of Galactic NSs, but there is no peak in the higher mass range (around 1.9M , for example). As for the Single Gaussian model, there is no significant peak. In comparison to the Uniform model, there is a moderate bulge in the low mass range (i.e., ∼ 1.2M − 1.6M ). The low mass cutoff of the BHs is constrained to 5.55 +0.96 −0.31 (assuming Truncated Power Law + Double Gaussian model), and it is rather similar in all the synthesis models. This estimated fraction r BNS is consistent with the merger rate densities of BNS (i.e., 320 +490 −240 Gpc −3 yr −1 obtained by Abbott et al. (2021b)) and NSBH (i.e., 130 +142 −69 Gpc −3 yr −1 obtained by Abbott et al. (2021a); see also Li et al. (2017) for a NSBH merger rate ≥ 100 Gpc −3 yr −1 based on the kilonova modeling). 3 We take the "O3actual"PSD, which is available at https://dcc.ligo.org/LIGO-T2000012/public, for the events except GW170817. Since GW170817 was detected in O2, the noise curve can be approximated by the Early High Sensitivity curve in Abbott et al. (2018), and can be found at https://git.ligo.org/lscsoft/bilby/-/tree/master/bilby/gw/detector/noise curves   Shao et al. (2020a). In each panel, the shaded region represents the 90% credible interval. As for the spin properties of the binaries, the posterior distribution obtained by the Truncated Power Law + the Double Gaussian model is shown in Fig. 2, and the results from other models are similar. Though the distribution of misalignment is not well constrained, it is quite clear that a perfectly aligned spin distribution (σ t = 0, ζ = 1) has been ruled out (see also Abbott et al. (2021c) for the same conclusion for BBHs), while a purely isotropic distribution of spin orientation (σ t > 2) is still allowed by the current data. The distribution of BHs' dimensionless spin magnitudes is constrained to µ BH a = 0.12 +0.18 −0.10 and σ BH a = 0.11 +0.11 −0.08 , but the distribution for NSs is poorly constrained. We do not expect to determine the population properties of spins of these kinds of binaries with currently limited observations, while the correlation between the spin parameters and the component masses of the compact binaries will influence the inference of the NS mass distribution with the GW data (Baird et al. 2013;Creswell et al. 2018;Pratten et al. 2020), so it is helpful to incorporate the spin information into our population inference.
We expect to have dozens of BNS/NSBH detections by the end of the fourth observing run of LIGO/Virgo/KAGRA network, and the number of events may reach one hundred within the duration of their fifth observing run (Abbott et al. 2018). Thus it is interesting to investigate whether we can determine the mass distribution of the NS via GWs in the next few years. In this section, we perform our analysis on mock GW detections generated from a prespecified underlying population. We assume the underlying NS mass distribution as that of the Galactic distribution, i.e., the Double Gaussian model (see Eq. (1)) with m min = 1M , m max = 2.25, µ 1 = 1.36M , σ 1 = 0.09M , µ 2 = 1.9M , σ 2 = 0.5M , and r = 0.65. The BH mass distribution is described by the Truncated Power Law model (see Eq. (4)) with m low = 5M , m up = 10M , and α = 3. The underlying spin distribution of binaries is described by Eq. (6), Eq. (7), and Eq. (8), with ζ = 0.8, σ t = 0.2, µ BH a = 0.1, σ BH a = 0.1, µ NS a = 0.005, and σ NS a = 0.01. Additionally, the fraction of BNSs in the surveyed population (including both BNSs and NSBHs) is assumed to 0.7. In generating mock detections, we assume that the underlying population follows a uniform in comoving volume and source frame time merger rate, with isotropic sky positions and inclinations. The true values of (m 1 , m 2 , a 1 , a 2 , z 1 , z 2 ) are drawn from the given distribution as described above. Additionally, φ 12 and φ jl (i.e., azimuthal angle between the spins of the two components, azimuthal angle between the total binary angular momentum and the orbital angular momentum) are set to be uniform in (0, 2π). The IMRPhenomXPHM waveform (London et al. 2018;Chatziioannou et al. 2017;Khan et al. 2019;García-Quirós et al. 2021) is used for generating simulated signals, and the identification of a GW signal is based on a network signal-to-noise ratio threshold ρ th = 12 and the design sensitivity noise curves 4 (Abbott et al. 2018).
To obtain the posteriors for each 'detected' event, we apply the user-friendly software Bilby (Ashton et al. 2019) with sampler Pymultinest (Buchner 2016) for parameter estimation. The sampling priors for single event parameter estimation are set as following: detector frame component masses are uniform in (0.9, 30)M , spin magnitudes are uniform in (0, 1)/(0, 0.05) for BHs/NSs, the spin orientations are isotropic, the priors of distance and inclination angle are set assuming the binaries are uniform in comoving volume, source frame time, and with isotropic directions. In order to accelerate the inference, we fix the other extrinsic parameters (i.e., the sky coordinate, the polarization angle, the coalescing phase, and the coalescing time) as their injected values.
We find with 100 'detected' events generated from the simulated population, we can rule out the Uniform distribution for NS mass by ln B ∼ 39 (compared with the Double Gaussian model), and the Double gaussian model is more favored than the Single Gaussian model by ln B ∼ 19 . The posteriors of hyperparameters for the Double Gaussian model are presented in Fig. 3, it is clear that the first peak of the NS mass distribution is well determined, while the second peak is ambiguous. The minimum and maximum masses of NS are both well constrained, additionally, the BNS fraction is constrained within uncertainty of 0.13 at 90% credible level. The distributions of NS mass modeled by the Double Gaussian and the Single Gaussian are shown in Fig. 4 (left). We can see that for the Double Gaussian model, the first peak is fairly significant, while the second peak is not obvious and it looks like a tail of the first peak. As for the Single Gaussian model, it fails to figure out the features of the injected population as is indicated by the Bayes factor mentioned above. The BH mass distribution is reconstructed by the Truncated Power Law, and all the hyperparameters are constrained well as presented in Fig. 4

CONCLUSION AND DISCUSSION
We perform hierarchical population inferences for the GW binaries involving at least one NS, and construct several synthesis models that include both the component masses and spin properties of the BNS and NSBH binaries, since there is a correlation between the component masses and the spin properties of compact binaries (Baird et al. 2013;Pratten et al. 2020;Abbott et al. 2021a). We have also performed inferences without spin information for comparison, and find that the corresponding evidence is slightly smaller than the inferences with spin information, e.g., ln B = ln Z with spin − ln Z no spin = 0.5 for Double Gaussian NS mass distribution model + Truncated Power Law BH mass distribution model 5 . The results are similar to those reported in Table 2, but the peak of Double Gaussian model in the spin-uninformative case is not as significant as that in the spin-informative case. Therefore, it is worthwhile to take the spin information into consideration when inferring the mass distribution of the compact binaries, though we do not aim to determine the spin properties of these kinds of binaries by currently limited observations. Our main conclusion is that the three NS mass distribution models (i.e., Double Gaussian, Single Gaussian, and Uniform) can not be reliably distinguished by the current GW data. Our result is similar to Landry & Read (2021), where it is found that the Flat (in this work we call the Uniform) and Bimodal (i.e., the Double Gaussian in this work) models are equally preferred using the Akaike information criterion (Akaike 1981). Nevertheless, there is a significant peak in the mass distribution obtained in our Double Gaussian model, and the location of the peak resembles that of Galactic NSs (Shao et al. 2020a). We further find that the mass distribution of Galactic NSs (Shao et al. 2020a) is slightly more preferred than the other models (by ln B ∼ 1.9), which indicates that the NSs observed via GW signals may have the similar mass distribution as the NSs identified by electromagnetic observation in the Galaxy. This conclusion is different from Landry & Read (2021), where the authors find that the mass distribution of the NSs observed via GWs is unimodal, predicting far fewer low-mass and moderately more high-mass NSs in the population. This difference may be caused by the inclusion of the spin information of compact binaries into our hierarchical population inference. Another reason is that Landry & Read (2021) use a uniform distribution for BH masses with a pairing function p(q) ∝ q βq for NSBHs, while we use the models for BH masses as indicated byÖzel et al. (2010). Both the BH mass distribution models indicate a narrow distribution located in ∼ 5 − 10M , which is consistent with the BHs in the Galactic X-ray binaries (Özel et al. 2010). As for the population of the compact binaries involving at least one NS, a fraction of ∼ 70% is BNSs, and this result agrees with the merger rate densities of BNS and NSBH estimated by Abbott et al. (2021b) and Abbott et al. (2021a). Due to the currently limited sample and the measurement uncertainties of NS spins, the distribution of NS spin magnitudes is not well constrained. From the inferred result of spin orientations, we conclude that a perfectly aligned spin distribution can be ruled out, but a purely isotropic distribution of spin orientation is still allowed by the GW data.
Encouragingly, LIGO/Virgo/KAGRA will start their O4 and O5 observing runs in the next few years, at that time the sensitivities of GW detectors will be significantly improved (Abbott et al. 2018), and the detection number of such binaries like NSBHs and BNSs will exceed one dozen in the duration of O4 and may reach one hundred by the end of O5. By performing simulations, we find that the mass distribution of NS can be well determined with a total number of 100 events (including NSBHs and BNSs) via GW observations. We do not characterize the mass ratios of the binaries in our models currently, which may be a potential improvement in determining population properties of NSs. For instance, very recently, Tang et al. (2021b) find that there is a mild evidence for a mass correlation among the two components of the low mass ratio binaries. Therefore, if such a feature exists, the population properties of NSs in the coalescing compact binaries can be characterized more accurately. With the observation of NSs in binary systems via GW, we can study the properties of NSs in their final moments, while the radio pulsar observations of Galactic NSs enable us to study the properties of NSs in their midlife. Combining measurements from GW and radio, provides us with a more complete understanding of compact binaries (involving at least one NS) from formation to merger (Galaudage et al. 2021a). Additionally, the population difference between the extragalactic NSs and the galactic NSs may have an impact on the determination of the origin of the heavy elements (Wang et al. 2017;Chen et al. 2021).
What's more, the reliable measurement of NS mass distribution, enables us to better understand the equation of state of dense matter Wysocki et al. 2020).