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 gravitational-wave (GW) observation data, we can rule out none of the double Gaussian, single Gaussian, and uniform NS mass distribution models, though a specific double Gaussian model inferred from the Galactic NSs is found to be slightly more preferred. The mass distribution of black holes (BHs) in the neutron star–black hole (NSBH) population is found to be similar to that in the Galactic X-ray binaries. Additionally, the ratio of the merger rate densities between NSBHs and BNSs is estimated to be ∼3:7. The spin properties of the binaries, though constrained relatively poorly, play a 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. To evaluate the feasibility of reliably determining the population properties of NSs in the coalescing compact binaries with upcoming GW observations, we perform simulations with a mock population. We find that with 100 detections (including BNSs and NSBHs) the mass distribution of NSs can be well determined, and the fraction of BNSs can also be accurately estimated.


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 Collaborations (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). Since the first successful detection of GWs (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. , 2020aAbbott et al. , 2020bAbbott et al. , 2020c. Very recently, The LIGO Scientific Collaboration et al. (2021) released a deep extended catalog, which reports another eight 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. , 2021Baxter et al. 2021;Safarzadeh & Wysocki 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 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 e and ∼ 1.9M e 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 an unprecedented valuable opportunity. 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 literature (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 (Rodriguez et al. 2016;Vitale et al. 2017). 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 Section 2, we introduce the data and the models used for inference, and in Section 3 we present the results. In Section 4, we carry out simulations with a mock population and evaluate the feasibility for reconstruction of the NS mass distribution with 100 BNSs/NSBHs detected in the design sensitivity run of Advanced LIGO-Virgo, and Section 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 is 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 multimessenger 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. 3 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.

Models
For the hierarchical population inference, here we consider three typical models for NS mass distribution and two for BH mass distribution.
In the above equations, m NS and m BH are the mass of NS and BH, while Φ 1 , Φ 2 , Φ, and F¢ 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 Equation (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 the spin reads where z cos 1,2 q = (θ 1,2 is the tilt angle between component spin and binary's orbital angular momentum), Φ z , a BH F , and are the normalization constants. The second term of Equation (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 a 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 is reasonable to set a 0.05 max NS = . All the descriptions 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, 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 a z m a z m a z m a z , , , , , , 2 2 2 1 1 1 where θ is the source parameters for each event, and 5M e , 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 four NS mass distribution models (i.e., double Gaussian, 4 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 eight synthesis models.  Note. Here, "U" means the uniform distribution.

Hierarchical Likelihood
With the data of the observed events {d} that are described in Section 2.1, and the population models described in Section 2.2, we perform a hierarchical Bayesian inference following Abbott et al. (2021c) and Thrane & Talbot (2019). For the given data {d} from the N det GW detections, the likelihood of Λ can be expressed as where ξ i (Λ) means the detection fraction, and the single-event likelihood d i i q ( | ) can be estimated using the posterior samples described in Section 2.1 (see Abbott et al. 2021c for details), then Equation (14) takes the form of d n  Tiwari (2018). And the injection campaigns are made by injecting simulated signal to the LIGO-Virgo detectors with noise curves. 5 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  and Dynesty sampler (Speagle 2020).

Results
In this section we report the results obtained by the hierarchical inference using synthesis models described in Section 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 evidence for the models a and b. Formally, the correct metric to compare two models is the odds we assume the prior odds π a and π b for the models are equal, so that b a b a =   . 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 more preferred than the others. It is certainly too early to conclude that the mass distribution of NSs with GWs 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 powerlaw model has a preference over the truncated Gaussian by ln 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 an NS mass distribution model of double Gaussian. Actually, the results assuming other BH mass distribution/NS mass distribution models are rather similar. In order to build more intuition, we plot the mass distributions of BHs/NSs as shown in Figure 1. For the double Gaussian model of NS mass distribution, it is clear that there is a peak around 1.35M e , 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 e , 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 e − 1.6M e ). The low-mass cutoff of the BHs is constrained to   , is consistent with that constrained by the observations of X-ray binaries in the Galaxy (Özel et al. 2010 Li et al. 2017 for an NSBH merger rate 100 Gpc −3 yr −1 based on the kilonova modeling).
As for the spin properties of the binaries, the posterior distribution obtained by the truncated power law + the double Gaussian model is shown in Figure 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. , 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.

Simulation
We expect to have dozens of BNS/NSBH detections by the end of the fourth observing run of the LIGO-Virgo-KAGRA network, and the number of events may reach 100 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 NSs 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 Equation (1) , μ 1 = 1.36M e , σ 1 = 0.09M e , μ 2 = 1.9M e , σ 2 = 0.5M e , and r = 0.65. The BH mass distribution is described by the truncated power-law model (see Equation (4)) with m low = 5M e , m up = 10M e, and α = 3. The underlying spin distribution of binaries is described by Equations (6) 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 , and z 2 ) are drawn from the given distribution as described above. Additionally, f 12 and f jl (i.e., azimuthal angle between the spins of the two components and azimuthal angle between the total binary angular momentum and the orbital angular momentum) are set to be uniform in (0, 2π). The IMRPhenomXPHM waveform (Chatziioannou et al. 2017;London et al. 2018;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 The right column shows the NS mass distributions obtained by double Gaussian, single Gaussian, and uniform models, respectively, and their partner model for BH mass is the truncated power law. The purple shaded region represents the mass distribution of the Galactic NSs as obtained by Shao et al. (2020a). In each panel, the shaded region represents the 90% credible interval.
network signal-to-noise ratio threshold ρ th = 12 and the design sensitivity noise curves 6 (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 follows: detector frame component masses are uniform in (0.9, 30)M e , 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 39  (compared with the double Gaussian model), and the double Gaussian model is more favored than the single Gaussian model by ln 19  . The posteriors of hyperparameters for the double Gaussian model are presented in Figure 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 an 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 Figure 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

Discussion and Conclusion
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., Z Z ln ln ln 0.5 with spin no spin = -=  for double Gaussian NS mass distribution model + truncated power-law BH mass distribution model. 7 The results are similar to those reported in Table 2, but the peak of the 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) cannot 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 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 µ b ( ) 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 e , 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 100 by the end of O5. By performing simulations, we find that the mass distribution of NSs 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) found that there is 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 GWs, 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 GWs 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). Furthermore, the reliable measurement of NS mass distribution enables us to better . Left: mass distributions of NSs inferred from the mock GW data with the double Gaussian and single Gaussian models. The black curves are the mass distribution of the injected population, and the shaded region represents the 90% credible interval. Right: posterior distribution for the BH hyperparameter parameters inferred from the mock GW data. The contours represent 50% and 90% credible bounds, respectively, and the solid lines stand for the injected hyperparameter values.