Investigating the distribution of double neutron stars and unconventional component mass

ABSTRACT

1 INTRODUCTION

Double neutron stars (DNSs) are important astronomical objects that can be used as natural laboratories for testing the theory of general relativity (GR), as well as studying gravitational waves (GWs; Burgay et al. 2003; Lyne et al. 2004; Stairs 2004; Abbott et al. 2017; Zhu & Ashton 2020; Boersma et al. 2021; Burgay, Perrodin & Possenti 2021). The characteristics of the mass distribution of DNSs can allow us to study the formation and evolution of their binary progenitor stars (Bhattacharya & van den Heuvel 1991; Miller 2002; Lorimer 2008; Tauris et al. 2013; Tauris, Langer & Podsiadlowski 2015; Tauris et al. 2017; van den Heuvel 2019). In particular, this reveals the dissimilar supernova explosion mechanisms and diverse evolutionary histories of DNSs (see, e.g. Podsiadlowski et al. 2005; Ferdman et al. 2013; Andrews et al. 2015; van den Heuvel 2017; van Leeuwen et al. 2020), which contain imprints from the late evolution of massive stars in binary systems (Dewi, Podsiadlowski & Sena 2006). The mass transfer of the binary evolution is useful for understanding the characteristics of the merging DNSs (Andrews 2020), and studying their Galactic birth rate will allow us to probe the DNS population. At present, through radio observations, 25 pairs of DNSs have been detected, 13 of which have measured masses for both NS components. Generally, the mass of the first-formed recycled NS is higher than that of the second-formed companion. However, some unconventional cases with reversed mass ratios exist (van den Heuvel 2017), such as PSR J0509+38011 (Lynch et al. 2018), PSR B1534+12 (Fonseca, Stairs & Thorsett 2014), and PSR J1757-1854 (Cameron et al. 2018). van den Heuvel (2017) hints at the possibility that the recycled NS has a lower mass than the non-recycled NS, because the former is a millisecond pulsar and is therefore spun up by the companion and yet has a lower mass than the companion.

The pioneering work on DNS mass statistics was initiated with the aforementioned systems. Through Bayesian inference, Finn (1994) found that the masses of DNS components mainly fall between 1.3 and 1.6 M_⊙. Then, a later statistical study of 19 systems, including also both white dwarfs and NSs and DNSs, showed that the NS masses follow a narrow Gaussian distribution (Thorsett & Chakrabarty 1999). To compare the two-component masses of DNSs, Nice, Stairs & Kasian (2008) analysed the masses of six DNSs with pairs of recycled and non-recycled components and found that the recycled pulsars tend to be more massive. Fostered by the increasing sample size, several ideas emerged that the formation mechanism of NSs was closely related to their mass distribution (Wong, Willems & Kalogera 2010; Kiziltan et al. 2013; Andrews et al. 2015). The NSs in DNS systems have evolved through different evolutionary paths, leading to distinctive signatures such as a different mass distribution peak and mass cut-off values. The distinct evolutionary formation scenarios are put forward to explain the bimodal distribution for 14 well-measured DNSs (Schwab, Podsiadlowski & Rappaport 2010). The statistical analysis of 18 DNS masses revealed that their average mass M = 1.32 ± 0.14 M_⊙ is systematically lower than the typical mass value of other pulsars (such as NSs and white dwarfs). This hints that the formation or evolutionary history of DNSs should differ from those of other binary systems (Zhang et al. 2011; Miller & Miller 2015; Özel & Freire 2016). Recently, statistical tests, such as Anderson–Darling and Mann–Whitney–Wilcoxon tests, have confirmed the bimodal distribution of the DNS mass, which can be attributed to the different formation mechanisms of the first- and second-born NSs in the pair (see Yang et al. 2019). In addition, the total DNS mass has also been studied. For example, Huang et al. (2018) analysed the total masses of 15 known Galactic DNS systems and suggested that the total mass distribution of DNSs follows a two-component Gaussian mixture model (GMM). However, Keitel (2019) re-studied the data used by Huang et al. (2018) through various statistical tests, and showed that a single statistical test cannot overcome the overfitting. So, he recommended that multiple tests should be considered when using the GMM for small samples of data and confirmed that a two-component mixture is consistent with the data. Farrow, Zhu & Thrane (2019) performed the Bayesian inference using a sample of 12 Galactic DNSs with precise component mass. Their results suggest that the mass parameters of recycled NSs follow a bimodality distribution, and uniform mass distribution for slow NSs. Moreover, Zhang et al. (2019) applied a Bayesian method and carried out Markov chain Monte Carlo simulations to study 15 Galactic DNS systems. They obtained the distribution range of the total masses, component masses, mass ratios, and chirp masses of DNS systems, and also explored the relationship between the mass distribution of DNSs and their merger sources. With the increasing sample of DNSs that have well-measured masses, it emerges that the non-recycled NSs could be more massive than their recycled pulsar companions. There are five such DNS systems, including PSR J0509+3801, PSR B1534+12 and PSR J1757−1854, and the corresponding data parameters are shown in Table 1. Particularly of note is PSR J1906+0746, which has been identified as a young pulsar with a characteristic age of about 113 kyr (Lorimer et al. 2006), and it is concluded that PSR J1906+0746 is a non-recycled companion with 1.291 M_⊙ based on orbital analysis (van Leeuwen et al. 2015).

In this study, we analysed the masses of pulsars (first-formed) and their companion NSs (second-formed) in 13 pairs of DNS systems by employing the GMM that is based on the Akaike information criterion (AIC) and the Bayesian information criterion (BIC). We made several improvements compared with previous studies, as follows. (i) Regarding the method, following the strategy and research ideas of Keitel (2019), we used the corrected AIC (AICc) and BIC methods that prove to be especially useful for small samples, although Farrow et al. (2019) did use the statistical test method known as the Bayes factor (BF). (ii) Regarding the results, our results showed that the mass distributions for recycled and non-recycled NSs prefer a double-Gaussian and single-Gaussian distribution, respectively. This is different from the results of Farrow et al. (2019) who presented a uniform distribution for non-recycled NS mass. (iii) Furthermore, we have calculated the percentage of DNSs with unconventional component mass, and we discuss the implications.

The paper is organized as follows. In Section 2, we investigate the mass distribution of the DNSs with GMMs, with a suite of statistical tests. Then, the probability of discovering unconventional DNS systems is given by the Monte Carlo sampling based on above model selection as described in Section 3. In Section 4, we discuss the birth rate of various DNS systems. In Section 5, we present a short conclusion.

2 METHOD

In order to investigate the distribution of the masses of DNSs, we employed the GMM, which allows us to distinguish the mass distribution types for recycled and non-recycled NSs. Then, we used several well-established statistical tests to validate the results, including the information standards (AICc and BIC) and the Bayesian evidence ratio, by which we can check the reliability of the mass distribution types. Untip un now, 25 pairs of DNS systems have been discovered, as listed in Table 1; the major data are taken from the Australia Telescope National Facility (ATNF) Pulsar Catalogue (see Manchester et al. 2005). Table 1 gives the physical parameters of the component masses of the recycled (M_p) and non-recycled (M_c) NSs, the total mass (M_t), orbital period (P_orb) and eccentricity (e) of the binaries, the spin period (P_s), spin period derivative (⁠|$\dot{P}$|⁠) and characteristic age. However, only 13 pairs of DNSs have a well-measured mass for the recycled and non-recycled NSs, so these are taken as our statistical sample.

In order to achieve the GMM fitting and model selection for the data set of DNSs, we assume that the masses of recycled pulsars and the non-recycled companions are independent. From previous studies, all known DNSs tend to follow a skewed distribution rather than a Gaussian distribution. There exist some non-negligible correlations between M_p and M_c for each DNS (i.e. they are not independent; Kiziltan et al. 2013). We take this into account in our analysis of the 13 pairs of DNSs with precise mass measures, which have a remarkably narrow mass range. Here, we use two independent statistical tests, similar to Keitel (2019), as depicted in the following subsections: (i) the information criteria (AIC and BIC; see Section 2.1); (ii) the Bayesian evidence ratios (see Section 2.2), which proved to be more robust than a simple likelihood ratio test that can suffer from overfitting, especially for a small group of samples, as discussed by Keitel (2019).

2.1 AICc and BIC

According to previous studies, by assuming that a sample of data points with errors (x_n) come from a Gaussian with mean μ_n and width σ_n, the improved likelihood function with N_comp (which is the number of components in the mixture model) and measurement errors can be expressed as (Mackay 2003; Ivezić et al. 2014; Keitel 2019)

where C_k ∈ [0, 1] is the component weight. In order to fit GMMs with N_comp = 1, 2, 3 components, we use the two independent Python packages (Keitel 2019), sklearngmm (Pedregosa et al. 2011) and xdgmm (Vanderplas et al. 2012; Ivezić et al. 2014; Holoien, Marshall & Wechsler 2017). Parameter estimates are consistent using sklearngmm and xdgmm, with the difference that sklearngmm supports basic multicomponent GMM fitting without measurement errors, but the full error treatment is important in assessing statistical robustness. The AICc for small sample sizes (Akaike 1974; Hurvich & Tsai 1989) and the BIC (Schwarz 1978) are described, respectively, as (see also Keitel 2019)

where N_coeffs is the number of coefficients of the model. For GMMs with N_comp = 1, 2, 3 components, respectively, the values of AICc and BIC are also collected in Table 2.

In general, lower values of AICc or BIC indicate a preferred model, and the significance of preference is judged by the differences between models (Burnham & Anderson 2002, 2004; Liddle 2007), that is, ΔIC. By convention, we here adopted the scale of Jeffreys (1961): ΔIC > 5 denotes ‘strong’ and ΔIC > 10 denotes ‘decisive’. The main results that we are interested in are those from xdgmm because this package can handle observational errors. For a more visual comparison, in Fig. 1 we provide scatter plots of ΔIC, which are quoted with respect to N_comp = 1. For the mass of the first-formed recycled pulsars, the AICc and BIC from xdgmm show that it slightly prefers a single- or double-component model (i.e. ΔAICc ≈ −0.5, ΔBIC ≈ −5.0, for N_comp = 2) and very strongly rejects a third-component model (i.e. ΔAICc ≈ 19.7, ΔBIC ≈ 0.81, for N_comp = 3). The inferences from sklearngmm agree with the previous xdgmm results, with slightly higher ΔIC values. Moreover, it is indecisive between the N_comp = 1 and N_comp = 2 cases, and it rules out the N_comp = 3 case decisively, based on the inferences of the AICc and BIC. For the mass of the non-recycled companions, the analysis process is similar to that of the mass of recycled pulsars; the lower values of ΔIC (i.e. ΔAICc ≈ 8.7, ΔBIC ≈ 4.2, for N_comp = 1) suggest that the GMM with N_comp = 1 is the optimal model within N_comp = 1−3 (see also the bottom panel of Fig. 1).

Figure 1.

Comparisons of ΔIC values for GMMs with N_comp = 1, 2, 3 applied to recycled pulsars (top panel) and non-recycled companions (bottom panel). The differences are quoted with respect to the N_comp = 1 value, and the negative ΔIC is preferred.

Open in new tabDownload slide

2.2 Bayesian evidence ratios

The component mass of a single DNS can be regarded as following a Gaussian distribution, approximately, based on previous studies (Kiziltan et al. 2013; Antoniadis et al. 2016). Each measurement of NS mass can be expressed as m_i ∼ N(μ, σ²). According to the Bayesian theorem, the posterior distribution for the parameter m_θ of the model |$\mathcal {H}$| can be defined as (refer to Jaynes 2003; Gregory 2005; Antoniadis et al. 2016)

where |$P\left(m_\theta \, |\, \, \mathcal {H},\mathcal {I}\right)$| is the prior distribution, |$P\left(m_i\, |\, m_\theta ,\mathcal {H},\mathcal {I}\right)$| is the likelihood function, and |$\mathcal {I}$| represents all other contextual information. Also, |$P\left(m_i\, |\, \, \mathcal {H},\mathcal {I}\right)$| is the Bayesian evidence for the model |$\mathcal {H}$| (i.e. the probability of the data by the model), given as

In order to achieve GMM fitting and to evaluate |$\mathcal {Z}_\mathcal {H}$| for different GMMs, for the purpose of model selection, we used the Python package cpnest (Skilling 2004; Veitch et al. 2017) similar to Keitel (2019). In our sampling process, we adopted |$N_{\, \mathrm{live}}=1024$| sample live points. After obtaining the convergent results of cpnest sampling, we used the Python package corner developed by Foreman-Mackey (2016) to plot the posterior distributions for the parameters, as shown in Appendix  A; their estimated values are also listed in Table 2.

Because the Bayesian evidence inference needs to use the same prior ranges for the corresponding parameters between the models (Jaynes 2003; Gregory 2005; Heavens 2009), our input priors for all parameters were selected to be broad and uniform over their corresponding intervals: C_k ∈ [0, 1], ρ_k ∈ [−1, 1], μ_k ∈ [0, 80], σ_k ∈ [0, 30], and constraints enforce ∑_kC_k = 1, C_{k + 1} ≥ C_k ∀ k for all models.

The model selection was performed by computing the Bayesian evidence ratio, also called the BF,

By convention, we adopted the criterion of Kass & Raftery (1995) to judge the preference for model |$\mathcal {H}_1$| over |$\mathcal {H}_2$|⁠: |$1< \mathrm{BF}^{\mathcal {H}_1}_{\mathcal {H}_2} < 3$| as ‘weak’, |$3< \mathrm{BF}^{\mathcal {H}_1}_{\mathcal {H}_2} < 20$| as ‘positive’, |$20< \mathrm{BF}^{\mathcal {H}_1}_{\mathcal {H}_2} < 150$| as ‘strong’, and |$\mathrm{BF}^{\mathcal {H}_1}_{\mathcal {H}_2} > 150$| as ‘very strong’.

By viewing the fitting results of the recycled NS mass distributions for the N_comp = 1, 2, 3 cases (see Table 2), and checking the values of the BF between the models, we obtained the following results: |$\mathrm{BF}^2_1 = \mathcal {Z}_2/\mathcal {Z}_1 \approx 3.56 \in \lbrace 3,5\rbrace$|⁠, |$\mathrm{BF}^2_3 = \mathcal {Z}_2/\mathcal {Z}_3 \approx 1.85 \in \lbrace 3,5\rbrace$|⁠. This implies that the mass distribution of pulsars slightly prefers a double component among the three models according to BF inference. Similarly, considering the companion mass distribution for the N_comp = 1, 2, 3 cases, the estimates are roughly consistent with the previous results of sklearngmm and xdgmm, and the results are more in favour of the single Gaussian model among the three cases based on the BF inference (⁠|$3 < \mathrm{BF}^1_2 = \mathcal {Z}_1/\mathcal {Z}_2 \approx 4.95 < 20$|⁠, |$3 < \mathrm{BF}^1_3 = \mathcal {Z}_1/\mathcal {Z}_3 \approx 12.42 < 20$|⁠).

3 RESULTS

As introduced above (see Section 2), it is yet to be determined which of the N_comp = 1 and N_comp = 2 cases better describes the distribution of component mass (M_p and M_c) of DNSs (see Fig. 2). However, the N_comp = 3 case can be confidently excluded, while applying the inferences of AICc, BIC, and Bayesian evidence to the M_p data set. It is suggested that the GMM with N_comp = 1 is the best-fitting model for the mass distribution of companions (M_c). It is worth noting that the mass distribution for DNSs hints at dissimilar supernova explosion mechanisms and different evolutionary histories (see, e.g. Podsiadlowski et al. 2004; Yang et al. 2019). The minimum mass magnitude of the remnant NS after a supernova explosion is shown to be about 1.17 M_⊙ based on the advanced stellar evolution calculations and neutrino-radiation hydrodynamics simulations for core-collapse supernova explosions. This is consistent with an asymmetric system of the recent detected DNS, PSR J0453+1559, with a value of about 1.174 M_⊙ (Martinez et al. 2015).

We randomly sampled DNS masses based on the model selection in Section 2, and the details of the Monte Carlo procedure are described below: 100 000 sample points (DNS) were selected from the data set of the first-formed NSs and their companions, in which the mass distribution of the first-formed NS follows a double-Gaussian distribution (also see Fig. 3), with two peaks around 1.34 ± 0.03 M_⊙ and 1.50 ± 0.26 M_⊙ respectively, while the mass distribution of the second-formed NS prefers a single Gaussian with a mean of 1.29 M_⊙ and a width of 0.08 M_⊙. The results of Fig. 4 showed that the probability of an unconventional mass DNS was about 27.7 per cent. In addition, it is also noticed that unconventional DNSs with q-factor less than 1 account for 27.7 per cent of the total sample, with the averaged (minimum/maximum) q-factor of 1.07 (0.78/1.62), which is similar to the observational value of 1.15 (0.92/1.33).

Figure 2.

Comparisons for mass distributions of the recycled and non-recycled NSs based on three model selections (SKLEARNGMM, XDGMM and CPNEST). Posterior predictive distributions of recycled NS mass M_p and non-recycled NS mass M_c are described in the case of N_comp = 1 and N_comp = 2, respectively. The former follows a double-Gaussian distribution, while the latter conforms to the single-Gaussian distribution.

Open in new tabDownload slide

Figure 3.

The mass distributions of recycled and non-recycled NS based on the Monte Carlo simulation. The area of the plot enclosed by the solid (dashed) line shows the double-Gaussian (single-Gaussian) distribution of recycled (non-recycled) NS.

Open in new tabDownload slide

Figure 4.

Histogram of the Monte Carlo DNS sampling q-factor that is defined as the mass ratio of a recycled NS to its companion NS. In the process of Monte Carlo sampling, the mass of the first-formed recycled NS obeys a double-Gaussian distribution, in the range of 1.25−1.62 M_⊙, while the mass of the second-formed non-recycled NS (companion) prefers a single-Gaussian distribution, in the range of 1.17–1.46 M_⊙. The masses of 100 000 pairs of DNSs ware randomly selected based on the above selection conditions.

Open in new tabDownload slide

4 DISCUSSION

In fact, the mass differences between the recycled pulsars and their companions for three systems (PSR J0509+3801, PSR B1534+12, and PSR J1757−1854) are very small. They are only 0.06, 0.01, and 0.06 M_⊙ for PSR J0509+3801, PSR B1534+12, and PSR J1757−1854, respectively. However, for DNSs that have conventional masses, such as PSR J0453+1559 and PSR J1913+1102, which have mass differences of 0.39 and 0.35 M_⊙, respectively, there seems to be a much wider range of masses, where the companion NS mass is smaller than the recycled pulsar mass. The DNS systems with unconventional component masses formed through common-envelope evolution and, during the formation process, they exchange mass more frequently.

We considered more robust model selection criteria, such as the information criteria (including AICc and BIC); the goal is to avoid overfitting by using a pure BF test for a small sample size. For the mass of the first-formed recycled pulsars, the result of AICc and BIC slightly prefers a double-component model compared with the single- and third-component models. For the mass of non-recycled pulsars, the AICc and BIC suggest that the mass distribution prefers a single-Gaussian distribution. The various GMM fitting methods employed here still all agree with consistent results to the small sample. The two-component model of the recycled NSs may imply that they have different mechanisms of formation. Low-mass NSs may form in electron-capture supernovae, and massive NSs (such as PSR J0453+1559 and J1913+1102) are more likely to form in an Fe core-collapse supernovae.

We used the same procedure described by Keitel (2019) to re-analyse the total masses of DNSs. The statistical results are shown in Appendix  B. The results of previous research based on the GMM show that the total mass distribution of DNSs follows a double-Gaussian distribution. Two subpopulations are thought to exist (Huang et al. 2018). However, Keitel (2019) argued that the preference was not obvious and that a single statistical test cannot overcome the overfitting problem, and so he recommended that multiple tests should be considered when using the GMM for small samples. Our result shows that the mass distributions for the total mass of 19 DNSs prefer a single-Gaussian distribution rather than a double-Gaussian distribution. As the number of samples increases, the statistical results tend to the single-Gaussian model. Keitel (2019) proposes that additional DNS discoveries might reveal the true statistical distribution. There is still some uncertainty for small samples of the total mass of DNSs. The question of subpopulations of the total mass distribution of DNSs would need to be re-analysed in the future.

5 CONCLUSIONS

We have presented a possible explanation for the formation of unconventional DNS systems, and we have suggested that the mass distribution of DNSs provides useful hints for understanding various origins of DNSs and the properties of their progenitors. Here we provide a short summary of the paper.

• A complete DNS data parameter table is established, based on which we calculate the merge times of all the detected DNS samples. We found that about half of them will inspiral and merge within the Hubble time.

• We studied and analysed the masses of pulsars (first-formed) and their companion NSs (second-formed) in 13 pairs of DNS systems with classical GMMs, This suggested that the mass parameters of recycled NSs follow a double-Gaussian distribution, in the range 1.25−1.62 M_⊙, and those of the non-recycled NSs prefer a single-Gaussian distribution, within the range 1.17–1.46 M_⊙.

• We randomly selected 100 000 recycled NSs and non-recycled companions using a Monte Carlo method based on the model selection and given distribution space, respectively. After sampling, we found that the fraction of unconventional DNS systems is about 27.7 per cent of the entire DNS population.

Acknowledgement

Special thanks are given to Xinjiang Zhu for help and useful conversations. This research is supported by the Center for Basic Science (Grant No. 11988101), the Guizhou Provincial Science and Technology Foundation (Grant No.[2022]164), the Subsidy project of the National Natural Science Foundation (Grant No. 2021GZJ006, U1938117), the Doctoral Project of Guizhou Education University (Grant No. 2020BS021), the International Partnership Program of Chinese Academy of Sciences No. 114A11KYSB20160008, and the National Key R&D Program of China No. 2016YFA0400702. We thank the anonymous referee for critical comments and suggestions that have significantly improved the quality of the paper.

DATA AVAILABILITY

The data for the DNSs underlying this paper are taken from the ATNF Pulsar Catalogue (Manchester et al. 2005), available at https://www.atnf.csiro.au/research/pulsar/psrcat/. This paper includes archived data obtained through the Parkes Pulsar Data archive on the CSIRO Data Access Portal (http://data.csiro.au).

References

APPENDIX A: THE POSTERIOR DISTRIBUTIONS

Figs A1 and A2 present the posterior distributions of the parameters of the cpnest results for the GMMs with N_comp = 1, 2 components, respectively; the estimates are also listed in Table 2. In each figure, the red vertical solid line indicates the median value and the other two vertical dashed lines encompass the 68 per cent (1σ level) credible interval. Here we present the posterior distributions for the best subhypothesis based on Gaussian mixture models and model selection. The result shows the posterior distributions of masses for the double-Gaussian distribution for recycled NS mass, and the normal distribution for non-recycled NS mass, respectively.

Figure A1.

The posterior distributions of the parameters for the normal distribution for non-recycled NS mass for the GMM with N_comp = 1.

Open in new tabDownload slide

Figure A2.

The posterior distributions of parameters for the double-Gaussian distribution for recycled NS mass for the GMM with N_comp = 2.

Open in new tabDownload slide

APPENDIX B: THE GAUSSIAN MIXTURE MODEL SELECTION OF DNS TOTAL MASSES

We analysed the distribution the total mass of the 19 DNS systems by employing the GMM based on the AIC and the BIC by following the strategy of Keitel (2019). Our analysis suggests that the total mass parameters prefer a single-Gaussian distribution, with the average and standard deviation as 2.65 ± 0.02 M_⊙ (also see Fig. B1).

Figure B1.

The posterior distributions of the total mass of DNSs for the GMM with N_comp = 1.

Open in new tabDownload slide

Figure B2.

Comparisons of ΔIC values for GMMs with N_comp = 1, 2, 3 applied to the total mass data set of DNSs. The differences are quoted with respect to the N_comp = 1 value, and the lower ΔIC is preferred.

Open in new tabDownload slide

Figure B3.

Posterior predictive distributions of total mass M_t. The area of the plot enclosed by the dot-dashed line shows the single-Gaussian distribution.

Open in new tabDownload slide

The GMM (including sklearngmm, xdgmm, and cpnest) results are given in Table B1. The N_comp = 2 and N_comp = 3 cases are not well constrained, and the posterior estimates are very broad. The lower ΔIC values will recommend N_comp = 1 (also see Fig. B2). The distributions of the total mass of DNSs based on the three model selections are shown in Fig. B3. All of the model selection criteria suggest that the posterior predictive distribution of the total mass of 19 DNSs follows a single-Gaussian distribution.