Population Properties of Neutron Stars in the Coalescing Compact Binaries

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 (Jin et al. 2015; 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. 2017, 2020a, 2020b, 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. 2020, 2021; Baxter 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; Farr & Chatziioannou 2020; 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 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. 2018, 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.

2.1. Selected Events

2.2. Models

For the hierarchical population inference, here we consider three typical models for NS mass distribution and two for BH mass distribution. The first NS mass function model is the double Gaussian scenario,

$$\begin{eqnarray}&&\begin{array}{l}\pi ({m}_{\mathrm{NS}}| {m}_{\min },{m}_{\max },{\mu }_{1},{\sigma }_{1},{\mu }_{2},{\sigma }_{2},{r}_{1})\\ =\ {r}_{1}{{ \mathcal N }}_{1}({m}_{\mathrm{NS}}| {\mu }_{1},{\sigma }_{1})/{{\rm{\Phi }}}_{1}+(1-{r}_{1}){{ \mathcal N }}_{2}({m}_{\mathrm{NS}}| {\mu }_{2},{\sigma }_{2})/{{\rm{\Phi }}}_{2},\\ \ \ \mathrm{for}\,{m}_{\mathrm{NS}}\in ({m}_{\min },{m}_{\max }).\end{array}\end{eqnarray} \tag{ 1 }$$

The second is the single Gaussian model,

$$\begin{eqnarray}&&\begin{array}{l}\pi ({m}_{\mathrm{NS}}| {m}_{\min },{m}_{\max },\mu ,\sigma )\\ =\,{ \mathcal N }({m}_{\mathrm{NS}}| \mu ,\sigma )/{\rm{\Phi }},\,\mathrm{for}\,{m}_{\mathrm{NS}}\in ({m}_{\min },{m}_{\max }).\end{array}\end{eqnarray} \tag{ 2 }$$

And the third is the uniform scenario,

$$\begin{eqnarray}&&\begin{array}{l}\pi ({m}_{\mathrm{NS}}| {m}_{\min },{m}_{\max })\\ =\,1/({m}_{\max }-{m}_{\min }),\,\mathrm{for}\,{m}_{\mathrm{NS}}\in ({m}_{\min },{m}_{\max }).\end{array}\end{eqnarray} \tag{ 3 }$$

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,

$$\begin{eqnarray}&&\begin{array}{l}\pi ({m}_{\mathrm{BH}}| {m}_{\mathrm{low}},{m}_{\mathrm{up}},\mu ,\sigma )\\ =\,{ \mathcal N }({m}_{\mathrm{BH}}| \mu ,\sigma )/{\rm{\Phi }}^{\prime} ,\,\mathrm{for}\,{m}_{\mathrm{BH}}\in ({m}_{\mathrm{low}},{m}_{\mathrm{up}}).\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\begin{array}{l}\pi ({m}_{\mathrm{BH}}| {m}_{\mathrm{low}},{m}_{\mathrm{up}},\alpha )\\ =\,\displaystyle \frac{{m}_{\mathrm{BH}}^{-\alpha }}{{\displaystyle \int }_{{m}_{\mathrm{low}}}^{{m}_{\mathrm{up}}}{m}^{-\alpha }{dm}},\,\mathrm{for}\,{m}_{\mathrm{BH}}\in ({m}_{\mathrm{low}},{m}_{\mathrm{up}}).\end{array}\end{eqnarray} \tag{ 5 }$$

In the above equations, m_NS and m_BH are the mass of NS and BH, while Φ₁, Φ₂, Φ, and ${\rm{\Phi }}^{\prime}$ are the normalization constants, and other parameters are described in Table 1.

Table 1. Distribution of the Priors for Hierarchical Bayesian Inference

Description	Parameters	NS Mass Distribution Priors
		Double Gaussian	Single Gaussian	Uniform
minimum mass of NSs	${m}_{\min }[{M}_{\odot }]$	U(0.9,1.3)	U(0.9,1.3)	U(0.9,1.3)
maximum mass of NSs	${m}_{\max }[{M}_{\odot }]$	U(1.7,2.9)	U(1.7,2.9)	U(1.7,2.9)
mean of first Gaussian	μ1[M⊙]	U(1, 2.2)	–	−
mean of second Gaussian	μ2[M⊙]	U(1, 2.2)	–	−
standard deviation of first Gaussian	σ1[M⊙]	U(0.01, 1)	–	−
standard deviation of second Gaussian	σ2[M⊙]	U(0.01, 1)	–	−
fraction of first Gaussian	r1	U(0, 1)	–	−
mean of single Gaussian	μ[M⊙]	–	U(1, 2.2)	−
standard deviation of single Gaussian	σ[M⊙]	–	U(0.01, 1)	−
fraction of BNS population	rBNS	U(0, 1)	U(0, 1)	U(0, 1)
constraint	–	${m}_{\min }\lt {\mu }_{1}\lt {\mu }_{2}\lt {m}_{\max }$	${m}_{\min }\lt \mu \lt {m}_{\max }$	−
	–	σ1 < σ2	−	−
		BH Mass Distribution Priors
		Truncated Power Law		Truncated Gaussian
low boundary of BH mass distribution	mlow[M⊙]	U(3,7)		U(3,7)
high boundary of BH mass distribution	mup[M⊙]	U(8,30)		U(8,30)
negative spectral index of power law	α	U( − 4, 12)		−
mean of Gaussian for BH	μBH[M⊙]	–		U(5, 15)
standard deviation of Gaussian for BH	σBH[M⊙]	–		U(0.5, 10)
constraint	–	–		mlow < μBH < mup
		Spin Priors
mean of spin for NS	${\mu }_{{\rm{a}}}^{\mathrm{NS}}$		U(0, 0.05)
standard deviation of spin for NS	${\sigma }_{{\rm{a}}}^{\mathrm{NS}}$		U(0.001, 0.05)
mean of spin for BH	${\mu }_{{\rm{a}}}^{\mathrm{BH}}$		U(0, 0.99)
standard deviation of spin for BH	${\sigma }_{{\rm{a}}}^{\mathrm{BH}}$		U(0.01, 0.25)
width of spin misalignment for field mergers	σt		U(0.01, 4. )
fraction of field mergers	ζ		U(0, 1)

Note. Here, "U" means the uniform distribution.

Download table as:  ASCIITypeset image

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

$$\begin{eqnarray}&&\pi (z| \zeta ,{\sigma }_{{\rm{t}}})=\zeta { \mathcal N }(z| 1,{\sigma }_{{\rm{t}}})/{{\rm{\Phi }}}_{{\rm{z}}}+(1-\zeta )/2,\,\mathrm{for}\,z\in (-1,1),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\begin{array}{l}\pi ({a}_{\mathrm{BH}}| {\mu }_{{\rm{a}}}^{\mathrm{BH}},{\sigma }_{{\rm{a}}}^{\mathrm{BH}})\\ \quad =\ { \mathcal N }({a}_{\mathrm{BH}}| {\mu }_{{\rm{a}}}^{\mathrm{BH}},{\sigma }_{{\rm{a}}}^{\mathrm{BH}})/{{\rm{\Phi }}}_{{\rm{a}}}^{\mathrm{BH}},\,\mathrm{for}\,{a}_{\mathrm{BH}}\in (0,1),\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\begin{array}{l}\pi ({a}_{\mathrm{NS}}| {\mu }_{{\rm{a}}}^{\mathrm{NS}},{\sigma }_{{\rm{a}}}^{\mathrm{NS}})\\ \quad =\ { \mathcal N }({a}_{\mathrm{NS}}| {\mu }_{{\rm{a}}}^{\mathrm{NS}},{\sigma }_{{\rm{a}}}^{\mathrm{NS}})/{{\rm{\Phi }}}_{{\rm{a}}}^{\mathrm{NS}},\,\mathrm{for}\,{a}_{\mathrm{NS}}\in (0,{a}_{\max }^{\mathrm{NS}}),\end{array}\end{eqnarray} \tag{ 8 }$$

where $z=\cos {\theta }_{\mathrm{1,2}}$ (θ_1,2 is the tilt angle between component spin and binary's orbital angular momentum), Φ_z, ${{\rm{\Phi }}}_{{\rm{a}}}^{\mathrm{BH}}$, and ${{\rm{\Phi }}}_{{\rm{a}}}^{\mathrm{NS}}$ 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}_{\max }^{\mathrm{NS}}=0.05$. 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,

$$\begin{eqnarray}&&\pi ({m}_{\mathrm{NS}},{a}_{\mathrm{NS}},z| {{\boldsymbol{\Lambda }}}_{{\bf{NS}}})=\pi ({m}_{\mathrm{NS}}| {{\boldsymbol{\Lambda }}}_{{\bf{NS}}}^{{\bf{m}}})\pi ({a}_{\mathrm{NS}}| {{\boldsymbol{\Lambda }}}_{{\bf{NS}}}^{{\bf{a}}})\pi (z| {{\boldsymbol{\Lambda }}}_{{\bf{z}}}),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\pi ({m}_{\mathrm{BH}},{a}_{\mathrm{BH}},z| {{\boldsymbol{\Lambda }}}_{{\bf{BH}}})=\pi ({m}_{\mathrm{BH}}| {{\boldsymbol{\Lambda }}}_{{\bf{BH}}}^{{\bf{m}}})\pi ({a}_{\mathrm{BH}}| {{\boldsymbol{\Lambda }}}_{{\bf{BH}}}^{{\bf{a}}})\pi (z| {{\boldsymbol{\Lambda }}}_{{\bf{z}}}),\end{eqnarray} \tag{ 10 }$$

where ${{\boldsymbol{\Lambda }}}_{{\bf{NS}}}^{{\bf{m}}}$ (${{\boldsymbol{\Lambda }}}_{{\bf{BH}}}^{{\bf{m}}}$), ${{\boldsymbol{\Lambda }}}_{{\bf{NS}}}^{{\bf{a}}}$ (${{\boldsymbol{\Lambda }}}_{{\bf{BH}}}^{{\bf{a}}}$), and Λ_z are the population parameters for mass distribution of NS (BH), spin magnitude of NS (BH), and spin orientations, respectively. And ${{\boldsymbol{\Lambda }}}_{{\bf{NS}}}={{\boldsymbol{\Lambda }}}_{{\bf{NS}}}^{{\bf{m}}}\cup {{\boldsymbol{\Lambda }}}_{{\bf{NS}}}^{{\bf{a}}}\cup {{\boldsymbol{\Lambda }}}_{{\bf{z}}}$ and ${{\boldsymbol{\Lambda }}}_{{\bf{BH}}}={{\boldsymbol{\Lambda }}}_{{\bf{BH}}}^{{\bf{m}}}\cup {{\boldsymbol{\Lambda }}}_{{\bf{BH}}}^{{\bf{a}}}\cup {{\boldsymbol{\Lambda }}}_{{\bf{z}}}$. Since the secondary objects of all the binaries are NSs, the model for them is

$$\begin{eqnarray}&&\pi ({m}_{2},{a}_{2},{z}_{2}| {{\boldsymbol{\Lambda }}}_{{\bf{2}}})=\pi ({m}_{\mathrm{NS}},{a}_{\mathrm{NS}},{z}_{2}| {{\boldsymbol{\Lambda }}}_{{\bf{NS}}}).\end{eqnarray} \tag{ 11 }$$

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

$$\begin{eqnarray}&&\begin{array}{l}\pi ({m}_{1},{a}_{1},{z}_{1}| {{\boldsymbol{\Lambda }}}_{{\bf{1}}})={r}_{\mathrm{BNS}}\pi ({m}_{\mathrm{NS}},{a}_{\mathrm{NS}},{z}_{1}| {{\boldsymbol{\Lambda }}}_{{\bf{NS}}})\\ \ +(1-{r}_{\mathrm{BNS}})\pi ({m}_{\mathrm{BH}},{a}_{\mathrm{BH}},{z}_{1}| {{\boldsymbol{\Lambda }}}_{{\bf{BH}}}).\end{array}\end{eqnarray} \tag{ 12 }$$

So the synthesis models for the NSBH and BNS binaries are represented by

$$\begin{eqnarray}&&\begin{array}{l}\pi (\theta | {\boldsymbol{\Lambda }})=\pi ({m}_{2},{a}_{2},{z}_{2},{m}_{1},{a}_{1},{z}_{1}| {\boldsymbol{\Lambda }})\\ =\,\pi ({m}_{2},{a}_{2},{z}_{2}| {{\boldsymbol{\Lambda }}}_{{\bf{2}}})\pi ({m}_{1},{a}_{1},{z}_{1}| {{\boldsymbol{\Lambda }}}_{{\bf{1}}}),\end{array}\end{eqnarray} \tag{ 13 }$$

where θ is the source parameters for each event, and Λ = Λ₁ ∪ Λ₂ = Λ_NS ∪ Λ_BH ∪ {r_BNS} is the population parameter for synthesis models. Additionally, in order to test whether the distribution of NSs observed via GWs is consistent with those observed electromagnetically in the Galaxy, we fix the parameters of double Gaussian for NSs as ${m}_{\min }=1{M}_{\odot },\,{m}_{\max }=2.25{M}_{\odot },\,{\mu }_{1}=1.36{M}_{\odot }$, σ₁ = 0.09M_⊙, μ₂ = 1.9M_⊙, σ₂ = 0.5M_⊙, andr₁ = 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, ⁴ 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.

2.3. 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

$$\begin{eqnarray}&&{ \mathcal L }(\{d\}| {\boldsymbol{\Lambda }})\propto \prod _{i=1}^{{N}_{\det }}\displaystyle \frac{\int { \mathcal L }({d}_{i}| {\theta }_{i})\pi ({\theta }_{i}| {\boldsymbol{\Lambda }})d{\theta }_{i}}{{\xi }_{i}({\boldsymbol{\Lambda }})},\end{eqnarray} \tag{ 14 }$$

where ξ_i (Λ) means the detection fraction, and the single-event likelihood ${ \mathcal L }({d}_{i}| {\theta }_{i})$ 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

$$\begin{eqnarray}&&{ \mathcal L }(\{d\}| {\boldsymbol{\Lambda }})\propto \prod _{i=1}^{{N}_{\det }}\displaystyle \frac{1}{{\xi }_{i}({\boldsymbol{\Lambda }})}\displaystyle \frac{1}{{n}_{i}}\sum _{k}^{{n}_{i}}\displaystyle \frac{\pi ({\theta }_{i}^{k}| {\boldsymbol{\Lambda }})}{{\pi }_{\phi }({\theta }_{i}^{k})},\end{eqnarray} \tag{ 15 }$$

where ${\pi }_{\phi }({\theta }_{i}^{k})$ 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. ⁵ 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).

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 (M_b ) as ${{ \mathcal B }}_{b}^{a}={Z}_{a}/{Z}_{b}$, 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 ratio ${{ \mathcal O }}_{b}^{a}=({Z}_{a}/{Z}_{b})\cdot ({\pi }_{a}/{\pi }_{b})$; we assume the prior odds π_a and π_b for the models are equal, so that ${{ \mathcal O }}_{b}^{a}={{ \mathcal B }}_{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 power-law model has a preference over the truncated Gaussian by $\mathrm{ln}{ \mathcal B }=\sim 0.5\mbox{--}1.5$.

Table 2. Logarithm Bayes Factors of the NS and BH Mass Distribution Models

$\mathrm{ln}{ \mathcal B }$	NS Mass Distribution Models
BH Mass Distribution Models	Double Gaussian	Single Gaussian	Uniform	Galactic
Truncated Gaussian	0	0.21	−0.15	1.37
Truncated Power Law	0.78	0.63	0.82	2.68

Note. The values are relative to the evidence of double Gaussian (NS mass distribution model) + truncated Gaussian (BH mass distribution model).

Download table as:  ASCIITypeset image

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_⊙, 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}_{-1.52}^{+0.96}{M}_{\odot }$ (${4.61}_{-1.37}^{+1.32}{M}_{\odot }$) by the truncated power-law (truncated Gaussian) model; however, the high-mass cutoff of BHs in NSBHs is constrained poorly by either model (i.e., ${17.41}_{-7.89}^{+11.00}{M}_{\odot }$ and ${14.26}_{-5.21}^{+13.59}{M}_{\odot }$ for truncated power-law and truncated Gaussian models). The spectrum index of truncated power law is ${5.95}_{-4.97}^{+4.70}$, which is consistent with the constraints of Abbott et al. (2021c), while the posterior support pushes to higher values. Interestingly, the mass distribution of BHs obtained by the truncated Gaussian, i.e., ${\mu }_{\mathrm{BH}}={6.97}_{-1.66}^{+2.92}{M}_{\odot }$ and ${\sigma }_{\mathrm{BH}}={3.47}_{-2.30}^{+5.39}{M}_{\odot }$, is consistent with that constrained by the observations of X-ray binaries in the Galaxy (Özel et al. 2010). The fraction of BNS mergers in the NSBH and BNS population is ${0.68}_{-0.31}^{+0.20}$ (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}_{-240}^{+490}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ obtained by Abbott et al. 2021b) and NSBH (i.e., ${130}_{-69}^{+142}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ obtained by Abbott et al. 2021a; see also Li et al. 2017 for an NSBH merger rate ≥ 100 Gpc⁻³ yr⁻¹ based on the kilonova modeling).

Zoom In Zoom Out Reset image size

Table 3. Inferred Population Parameters of the NS and BH Mass Distribution Models

	Parameters	NS Mass Distribution Models Assuming Truncated Power Law BH Mass Distribution Model
		Double Gaussian	Single Gaussian	Uniform
	${m}_{\min }({M}_{\odot })$	${1.12}_{-0.18}^{+0.15}$	${1.12}_{-0.18}^{+0.14}$	${1.16}_{-0.19}^{+0.12}$
	${m}_{\max }({M}_{\odot })$	${2.13}_{-0.29}^{+0.67}$	${2.10}_{-0.28}^{+0.71}$	${1.94}_{-0.16}^{+0.25}$
	μ1(μ)(M⊙)	${1.35}_{-0.22}^{+0.25}$	${1.42}_{-0.26}^{+0.36}$	−
	μ2(M⊙)	${1.68}_{-0.33}^{+0.36}$	–	−
	σ1(σ)(M⊙)	${0.28}_{-0.23}^{+0.40}$	${0.47}_{-0.25}^{+0.42}$	−
	σ2(M⊙)	${0.67}_{-0.39}^{+0.29}$	–	−
	r1	${0.58}_{-0.46}^{+0.37}$	–	−
	rBNS	${0.68}_{-0.31}^{+0.20}$	${0.67}_{-0.30}^{+0.20}$	${0.67}_{-0.31}^{+0.21}$
		BH Mass Distribution Models Assuming NS Mass Distribution Model of Double Gaussian
		Truncated Power Law		Truncated Gaussian
	mlow(M⊙)	${5.55}_{-1.52}^{+0.96}$		${4.61}_{-1.37}^{+1.32}$
	mup(M⊙)	${17.41}_{-7.89}^{+11.00}$		${14.26}_{-5.21}^{+13.59}$
	α	${5.95}_{-4.97}^{+4.70}$		−
	μBH(M⊙)	–		${6.97}_{-1.66}^{+2.92}$
	σBH(M⊙)	–		${3.47}_{-2.30}^{+5.39}$

Note. The values represent the medians and symmetric 90% credible intervals of the parameters.

Download table as:  ASCIITypeset image

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. The distribution of BHs' dimensionless spin magnitudes is constrained to ${\mu }_{{\rm{a}}}^{\mathrm{BH}}={0.12}_{-0.10}^{+0.18}$ and ${\sigma }_{{\rm{a}}}^{\mathrm{BH}}={0.11}_{-0.08}^{+0.11}$, 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.

Zoom In Zoom Out Reset image size

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)) with ${m}_{\min }=1{M}_{\odot }$, ${m}_{\max }=2.25$, μ₁ = 1.36M_⊙, σ₁ = 0.09M_⊙, μ₂ = 1.9M_⊙, σ₂ = 0.5M_⊙, and r = 0.65. The BH mass distribution is described by the truncated power-law model (see Equation (4)) with m_low = 5M_⊙, m_up = 10M ⊙, and α = 3. The underlying spin distribution of binaries is described by Equations (6), (7), and (8), with ζ = 0.8, σ_t = 0.2, ${\mu }_{{\rm{a}}}^{\mathrm{BH}}=0.1$, ${\sigma }_{{\rm{a}}}^{\mathrm{BH}}=0.1$, ${\mu }_{{\rm{a}}}^{\mathrm{NS}}=0.005$, and ${\sigma }_{{\rm{a}}}^{\mathrm{NS}}=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₁, m₂, a₁, a₂, z₁, and z₂) are drawn from the given distribution as described above. Additionally, ϕ₁₂ and ϕ_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 network signal-to-noise ratio threshold ρ_th = 12 and the design sensitivity noise curves ⁶ (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_⊙, 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 $\mathrm{ln}{ \mathcal B }\sim 39$ (compared with the double Gaussian model), and the double Gaussian model is more favored than the single Gaussian model by $\mathrm{ln}{ \mathcal B }\sim 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 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 Figure 4 (right). However, the spin properties of binaries are not well constrained, where ${\mu }_{a}^{\mathrm{NS}}={0.02}_{-0.004}^{+0.003}$, ${\sigma }_{a}^{\mathrm{NS}}={0.004}_{-0.002}^{+0.003}$, ${\mu }_{a}^{\mathrm{BH}}={0.168}_{-0.022}^{+0.019}$, ${\sigma }_{a}^{\mathrm{BH}}={0.071}_{-0.016}^{+0.024}$, ${\sigma }_{t}={0.87}_{-0.14}^{+0.20}$, and $\zeta ={0.93}_{-0.19}^{+0.06}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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., $\mathrm{ln}{ \mathcal B }=\mathrm{ln}{Z}_{\mathrm{with}\,\mathrm{spin}}-\mathrm{ln}{Z}_{\mathrm{no}\,\mathrm{spin}}=0.5$ for double Gaussian NS mass distribution model + truncated power-law BH mass distribution model. ⁷ 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 $\mathrm{ln}{ \mathcal B }\sim 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)\propto {q}^{{\beta }_{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 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 understand the equation of state of dense matter (Chatziioannou & Farr 2020; Wysocki et al. 2020).

This work was supported in part by NSFC under grant Nos. 11921003, 11703098, and 12073080, the Chinese Academy of Sciences via the Strategic Priority Research Program (grant No. XDB23040000), and Key Research Program of Frontier Sciences (No. QYZDJ-SSW-SYS024). This research has made use of data and software obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration, and the Virgo Collaboration. LIGO is funded by the U.S. National Science Foundation. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN), and the Dutch Nikhef, with contributions by Polish and Hungarian institutes.

Software: Bilby (Ashton et al. 2019, https://git.ligo.org/lscsoft/bilby/), Dynesty (Speagle 2020, https://github.com/joshspeagle/dynesty), PyMultiNest (Buchner 2016, https://github.com/JohannesBuchner/PyMultiNest), PyCBC (Biwer et al. 2019; Nitz et al. 2021, https://github.com/gwastro/pycbc).