The Mass Distribution of Galactic Double Neutron Stars

The breakthrough discovery of GW170817 was the first ever detection of gravitational waves arising from a neutron star (NS) merger, 44 Mpc from Earth in the galaxy NGC4993 (Abbott et al. 2017; Coulter et al. 2017). More binary NS events are expected when Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO; Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2015) detectors resume observing with improved sensitivity this year (Abbott et al. 2018). In our own Galaxy, the first double neutron star (DNS) system, PSR B1913+16, was discovered in 1974 (Hulse & Taylor 1975). The number of known Galactic DNS systems has increased to 18 so far (see Zhu et al. 2018, and references therein), with three new discoveries published in 2018 (Cameron et al. 2018; Lynch et al. 2018; Stovall et al. 2018).

The mass distribution of Galactic DNSs has been studied by numerous groups. Özel et al. (2012) and Kiziltan et al. (2013) fit the mass measurements for nine binaries to a Gaussian distribution. This Gaussian model was used by Taylor & Gair (2012) to investigate the prospect of constraining cosmological parameters using gravitational-wave observations alone, and by Zhu et al. (2013) to compute the gravitational-wave background from a population of NS mergers. Pejcha et al. (2012) demonstrated that the observed mass distribution contains information about the supernova explosion mechanism. Özel & Freire (2016) updated the initial Gaussian fit using measurements of two additional systems and found the most likely values of the mean and width to be μ = 1.33 M_⊙ and σ = 0.09 M_⊙.

Since the discovery of GW170817, the mass distribution of Galactic DNSs has seen applications in a wider range of research. Specifically, the conventional model of Özel & Freire (2016) was used (1) as an astrophysical prior to place constraints on the NS equation of state using data from GW170817 (De et al. 2018), (2) to argue that GW170817 is unlikely to have come from a population of DNSs like those found in our Galaxy (Pankow 2018), (3) to investigate the possibility of detecting a population of postmerger signals from NS mergers (Yang et al. 2018), and (4) as a prior to argue that fast radio bursts may be associated with NS mergers resulting in prompt collapse (Paschalidis & Ruiz 2018). Additionally, Ma et al. (2018) adopted the single-Gaussian models derived in Özel et al. (2012) and Kiziltan et al. (2013) to investigate the prospect of forming supramassive NSs from NS mergers. Zhu et al. (2018) developed a framework to infer DNS population properties through gravitational-wave observations, assuming the masses and spin periods of Galactic DNSs are representative of the merging population.

In light of the new discoveries of Galactic DNSs from radio astronomy and the ever-expanding catalog of gravitational-wave events, we revisit the problem of inferring the mass distribution of Galactic DNS systems. Apart from using a sample of observations that is nearly twice the size of earlier data sets, we expand on previous studies in several ways. First, we enlarge the space of models so that the recycled NS and its nonrecycled companion do not necessarily follow the same mass distribution. Second, we consider different functional forms for the mass distribution to include uniform, Gaussian, and two-Gaussian distributions. Third, we account for the correlation between mass measurements for pulsars and their companions, which appears to be ignored in previous work.

This paper is organized as follows. In Section 2, we review the mass measurements of Galactic DNS systems. In Section 3, we describe the formalism for Bayesian inference and model selection. In Section 4, we present and discuss our analysis results. Finally, we conclude in Section 5.

Among 18 Galactic DNS systems, 17 have reported measurements of the binary total mass, 12 of which have masses measured for both component stars; see Table 3 in Zhu et al. (2018) and references therein for details. In this work, we divide the NSs in these binaries into two categories: recycled NSs and nonrecycled ("slow") NSs, for which the masses are labeled as m_r and m_s, respectively.

In the standard isolated binary formation channel (e.g., Tauris et al. 2017), a recycled NS is the first born object that gets spun up to ∼10–100 ms through an accretion/recycling process, whereas the second born NS quickly (∼Myr) spins down to a period of ${ \mathcal O }(1)$ second after its birth. The double pulsar PSR J0737−3039A/B is an excellent example, with pulsars A and B representing recycled and slow NSs, respectively (Burgay et al. 2003; Lyne et al. 2004). All pulsars in the remaining 16 binaries are recycled NSs, except PSR J1906+0746, which is a young slow pulsar (van Leeuwen et al. 2015). Therefore, the companion of PSR J1906+0746 is included in the category of recycled NSs. Although the possibility of a massive white dwarf companion cannot be ruled out from radio observations, we take PSR J1906+0746 to be in a DNS system, as its orbital characteristics and mass estimates are consistent with other Galactic DNSs, especially the double pulsar (Yang et al. 2017).

Note that PSR B2127+11C and PSR J1807−2500B are located in globular clusters, whereas the rest are in the Galactic disk. It has been hypothesized that these two DNSs were formed from a different evolutionary scenario than the one described above (Prince et al. 1991; Lynch et al. 2012). In addition, the companion of PSR J1807−2500B could also be a massive white dwarf. Nevertheless, we include them in the analysis. As can be seen in Figure 1, their masses are consistent with other DNS systems. Therefore, we do not expect such a choice to significantly bias our results.

Zoom In Zoom Out Reset image size

Table 1 lists the measurements or constraints of the binary total mass M_T, mass of the recycled (m_r), and slow (m_s) NS for 17 binaries. To aid the comparison with gravitational-wave measurements, we also list values for the binary chirp mass, defined as

$$\begin{eqnarray}&&{{ \mathcal M }}_{c}=\displaystyle \frac{{\left({m}_{1}{m}_{2}\right)}^{3/5}}{{\left({m}_{1}+{m}_{2}\right)}^{1/5}},\end{eqnarray} \tag{ 1 }$$

and mass ratio q = m₂/m₁ ≤ 1.

Table 1.  Mass Measurements of Galactic DNS Systems

Pulsar Name	MT (M⊙)	mr (M⊙)	ms (M⊙)	${{ \mathcal M }}_{c}$ (M⊙)	q	Pb (day)	Tc (Gyr)	References
Systems Will Merge within a Hubble Time
J1946+2052	2.50(4)	<1.35	>1.17	(1.05, 1.11)	(0.68, 1)	0.078	0.046	(1)
J1756−2251	2.56999(6)	1.341(7)	1.230(7)	1.1178(3)	0.92(1)	0.320	1.656	(2)
J0737−3039A/B	2.58708(16)	1.3381(7)	1.2489(7)	1.1253(1)	0.933(1)	0.102	0.086	(3)
J1906+0746	2.6134(3)	1.322(11)	1.291(11)	1.1372(2)	(0.956, 1)	0.166	0.308	(4)
B1534+12	2.678463(4)	1.3330(2)	1.3455(2)	1.165870(2)	0.9907(3)	0.421	2.734	(5)
B2127+11C	2.71279(13)	1.358(10)	1.354(10)	1.18043(8)	(0.975, 1)	0.335	0.217	(6)
J1757−1854	2.73295(9)	1.3384(9)	1.3946(9)	1.18930(4)	0.960(1)	0.184	0.076	(7)
J0509+3801	2.805(3)	1.34(8)	1.46(8)	1.215(5)	(0.793, 1)	0.380	0.574	(8)
B1913+16	2.828378(7)	1.4398(2)	1.3886(2)	1.230891(5)	0.9644(3)	0.323	0.301	(9)
J1913+1102	2.875(14)	1.64(4)	1.24(4)	1.239(8)	0.76(4)	0.206	0.475	(10)
Systems Will Not Merge within a Hubble Time
J1807−2500B	2.57190(73)	1.3655(21)	1.2064(21)	1.1169(3)	0.883(3)	9.957	1044	(11)
J1518+4904	2.7183(7)	1.41(8)	1.31(8)	1.181(5)	(0.794, 1)	8.634	8832	(12)
J0453+1559	2.733(4)	1.559(5)	1.174(4)	1.175(2)	0.753(5)	4.072	1453	(13)
J1411+2551	2.538(22)	<1.64	>0.92	(1.05, 1.11)	(0.57, 0.95)	2.616	466	(14)
J1811−1736	2.57(10)	<1.75	>0.91	(1.02, 1.17)	(0.58, 0.95)	18.78	1794	(15)
J1829+2456	2.59(2)	<1.36	>1.25	(1.08, 1.14)	(0.65, 1)	1.176	55	(16)
J1930−1852	2.59(4)	<1.32	>1.30	(1.07, 1.15)	(0.58, 0.96)	45.06	∼105	(17)

Note. The binary total mass (M_T), the masses of the recycled NS (m_r) and the slow NS (m_s), binary chirp mass (${{ \mathcal M }}_{c}$), mass ratio (q), binary orbital period (P_b), and coalescence time T_c. Figures in parentheses are 1σ uncertainties in the last quoted digit of the mass measurement. For five systems without component mass measurements, we list the 99% confidence upper and lower limits for m_r and m_s, derived from their measured M_T and mass functions: 0.268184 (J1946+2052), 0.1223898 (J1411+2551), 0.128121 (J1811−1736), 0.29413 (J1829+2456), and 0.34690765 (J1930−1852). For ${{ \mathcal M }}_{c}$ and for q, if the distribution is non-Gaussian or extended, 90% confidence intervals are given.

References. (1) Stovall et al. (2018); (2) Ferdman et al. (2014); (3) Kramer et al. (2006); (4) van Leeuwen et al. (2015); (5) Fonseca et al. (2014); (6) Jacoby et al. (2006); (7) Cameron et al. (2018); (8) Lynch et al. (2018); (9) Weisberg et al. (2010); (10) Ferdman & PALFA Collaboration (2018); (11) Lynch et al. (2012); (12) Janssen et al. (2008); (13) Martinez et al. (2015); (14) Martinez et al. (2017); (15) Corongiu et al. (2007); (16) Champion et al. (2005); (17) Swiggum et al. (2015). In addition, m_r and m_s values of PSR J1518+4904 were taken from Tauris et al. (2017) who cite a private communication with G. Janssen.

Download table as:  ASCIITypeset image

Figure 1 shows the joint distribution of (m_r, m_s) for 17 Galactic DNS systems, along with the mass constraints of GW170817 (Abbott et al. 2019). Figure 2 shows the marginalized distributions of m_r and m_s. We briefly describe these mass distributions below; details of NS mass measurements through pulsar timing can be found in Stairs (2003).

Zoom In Zoom Out Reset image size

As can be seen in Figure 1, the measurements of m_r and m_s are anticorrelated due to the constraint on M_T. In the left plot are 12 systems for which precise component mass measurements are available, whereas on the right are another 5 binaries with only total mass measurements. For the 12 systems, (m_r, m_s) follows a bivariate normal distribution³ with a covariance of $({\sigma }_{T}^{2}-{\sigma }_{r}^{2}-{\sigma }_{s}^{2})/2$, where σ_T, σ_r, and σ_s are the measurement errors quoted in parentheses in Table 1 for M_T, m_r, and m_s, respectively.

For the five systems shown in the right panel of Figure 1, the joint distribution of (m_r, m_s) is determined by measurements of M_T and the mass function

$$\begin{eqnarray}&&\begin{array}{l}p({m}_{r},{m}_{s})\propto \exp \left[-\displaystyle \frac{{\left({m}_{r}+{m}_{s}-{\hat{M}}_{T}\right)}^{2}}{2{\sigma }_{T}^{2}}\right]\\ \quad \times \,{\displaystyle \int }_{-1}^{1}\delta (f-\hat{f})d(\cos i)\\ \quad \propto \,\displaystyle \frac{\exp \left[-\tfrac{{\left({m}_{r}+{m}_{s}-{\hat{M}}_{T}\right)}^{2}}{2{\sigma }_{T}^{2}}\right]}{\tfrac{3{m}_{s}^{3}}{{\left({m}_{r}+{m}_{s}\right)}^{2}}{\left[\tfrac{\hat{f}{\left({m}_{r}+{m}_{s}\right)}^{2}}{{m}_{s}^{3}}\right]}^{\tfrac{1}{3}}{\left[1-\tfrac{{\hat{f}}^{2/3}{\left({m}_{r}+{m}_{s}\right)}^{4/3}}{{m}_{s}^{2}}\right]}^{\tfrac{1}{2}}},\end{array},\end{eqnarray} \tag{ 2 }$$

where $\hat{f}$ and ${\hat{M}}_{T}$ are measured mean values of the mass function and total mass, respectively. The mass function is defined as

$$\begin{eqnarray}&&f=\displaystyle \frac{{\left({m}_{c}\sin i\right)}^{3}}{{M}_{T}^{2}}={\left(\displaystyle \frac{2\pi }{{P}_{b}}\right)}^{2}\displaystyle \frac{{\left(a\sin i\right)}^{3}}{G},\end{eqnarray} \tag{ 3 }$$

where P_b, a, and i are the orbital period, semimajor axis, and orbital inclination angle, respectively; m_c is the mass of the pulsar's companion, which is replaced with m_s in Equation (2) because for the systems in question, the companion is a slow NS. As in Alsing et al. (2018), we ignore the measurement uncertainty of mass function to derive the last line in Equation (2).

3.1. Constraining the Shape of NS Mass Distributions with Hyperparameters

Our starting point is the measurement of NS masses from radio or gravitational-wave observations.⁴ These previously published constraints on NS masses are posterior distributions, which we denote as

$$\begin{eqnarray}&&p({m}_{r}^{i},{m}_{s}^{i}| {{\rm{d}}}^{i}),\end{eqnarray} \tag{ 4 }$$

where i refers to the ith binary. Formally, this is the joint posterior on the recycled NS mass m_r and the slow NS mass m_s given some data d, as illustrated in Figure 1.

We use the posterior distributions from published NS mass constraints to study the population properties of NS using hierarchical Bayesian inference. We introduce priors for m_r and m_s, which are conditional on hyperparameters $({{\rm{\Lambda }}}_{r},{{\rm{\Lambda }}}_{s})\in {\rm{\Lambda }}$

$$\begin{eqnarray}&&\pi ({m}_{r}^{i}| {{\rm{\Lambda }}}_{r},H),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\pi ({m}_{s}^{i}| {{\rm{\Lambda }}}_{s},H).\end{eqnarray} \tag{ 6 }$$

The hyperparameters (Λ_r, Λ_s) describe the shape of the (m_r, m_s) priors given a hypothesis H. For example, later we consider a hypothesis in which $\pi ({m}_{r}| {{\rm{\Lambda }}}_{r},H)$ is a Gaussian distribution in m_r. In this case, Λ_r consists of hyperparameters for the mean ${\mu }_{r}\in {{\rm{\Lambda }}}_{r}$ and width ${\sigma }_{r}\in {{\rm{\Lambda }}}_{r}$ of the distribution.

Our first goal is to calculate (hyper)posterior distributions for (Λ_r, Λ_s). Given a data set with N binary mass measurements, the hyperposterior distribution is given by (Thrane & Talbot 2019)

$$\begin{eqnarray}&&\begin{array}{r}\begin{array}{l}p({\rm{\Lambda }}| \left\{{\rm{d}}\right\})=\displaystyle \frac{1}{{{Z}}_{{\rm{\Lambda }}}}\displaystyle \prod _{i}^{N}\displaystyle \int {dm}_{r}^{i}\displaystyle \int {dm}_{s}^{i}\,{L}({{\rm{d}}}^{i}| {m}_{r}^{i},{m}_{s}^{i})\\ \,\pi ({m}_{r}^{i}| {{\rm{\Lambda }}}_{r})\pi ({m}_{s}^{i}| {{\rm{\Lambda }}}_{s})\pi ({{\rm{\Lambda }}}_{r})\pi ({{\rm{\Lambda }}}_{s}).\end{array}\end{array}\end{eqnarray} \tag{ 7 }$$

Here π(Λ_r) and π(Λ_s) are the priors of our hyperparameters, which we take to be flat. The term ${{ \mathcal Z }}_{{\rm{\Lambda }}}$ is the (hyper) evidence

$$\begin{eqnarray}&&\begin{array}{r}\begin{array}{l}{{Z}}_{{\rm{\Lambda }}}\equiv \displaystyle \int d{{\rm{\Lambda }}}_{s}\displaystyle \int d{{\rm{\Lambda }}}_{r}\,\displaystyle \prod _{i}^{N}\displaystyle \int {dm}_{r}^{i}\displaystyle \int {dm}_{s}^{i}\\ \,{L}({{\rm{d}}}^{i}| {m}_{r}^{i},{m}_{s}^{i})\pi ({m}_{r}^{i}| {{\rm{\Lambda }}}_{r})\pi ({m}_{s}^{i}| {{\rm{\Lambda }}}_{s})\pi ({{\rm{\Lambda }}}_{r})\pi ({{\rm{\Lambda }}}_{s}).\end{array}\end{array}\end{eqnarray} \tag{ 8 }$$

The variable ${L}({{\rm{d}}}^{i}| {m}_{r}^{i},{m}_{s}^{i})$ is the likelihood function of the data given $({m}_{r}^{i},{m}_{s}^{i})$. We do not have direct access to this likelihood, but it is related to the posterior distribution

$$\begin{eqnarray}&&{L}({{\rm{d}}}^{i}| {m}_{r}^{i},{m}_{s}^{i})=\displaystyle \frac{{{Z}}_{0}^{i}}{{\pi }_{0}({m}_{r}^{i},{m}_{s}^{i})}p({m}_{r}^{i},{m}_{s}^{i}| {{\rm{d}}}^{i}),\end{eqnarray} \tag{ 9 }$$

where ${\pi }_{0}({m}_{r},{m}_{s})$ is the initial prior used to derive the mass posterior (assumed here to be flat). Meanwhile,

$$\begin{eqnarray}&&{{Z}}_{0}\equiv \int {dm}_{r}\int {dm}_{s}\,{L}({\rm{d}}| {{\rm{m}}}_{{\rm{r}}},{{\rm{m}}}_{{\rm{s}}}){\pi }_{0}({{\rm{m}}}_{{\rm{r}}},{{\rm{m}}}_{{\rm{s}}}),\end{eqnarray} \tag{ 10 }$$

is the initial evidence. Without access to the raw radio data d, we do not know ${{ \mathcal Z }}_{0}$, but this factor can be ignored since it will ultimately cancel, via either normalization of the hyperposterior or the construction of a Bayesian evidence ratio. For the sake of readability, we therefore set ${{ \mathcal Z }}_{0}^{i}=1$ henceforth.

Plugging Equation (9) into (7), we obtain

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}p({\rm{\Lambda }}| \left\{{\rm{d}}\right\}) & = & \displaystyle \frac{1}{{{Z}}_{{\rm{\Lambda }}}}\displaystyle \prod _{i}^{N}\displaystyle \int {dm}_{r}^{i}\displaystyle \int {dm}_{s}^{i}\,p({m}_{r}^{i},{m}_{s}^{i}| {{\rm{d}}}^{i})\\ & & \times \displaystyle \frac{\pi ({m}_{r}^{i}| {{\rm{\Lambda }}}_{r})\pi ({m}_{s}^{i}| {{\rm{\Lambda }}}_{s})}{{\pi }_{0}({m}_{r}^{i},{m}_{s}^{i})}\pi ({{\rm{\Lambda }}}_{r})\pi ({{\rm{\Lambda }}}_{s}).\end{array}\end{array}\end{eqnarray} \tag{ 11 }$$

Using Equation (11), we construct posteriors on the hyperparameters that describe the shape of NS mass distributions.

In this work, we employ the pymultinest implementation (Buchner et al. 2014) of the MultiNest algorithm (Feroz et al. 2009) for the stochastic sampling of posteriors and the evidence calculation. Our codes are publicly available on GitHub, https://github.com/NicholasFarrow/GalacticDNSMass.

3.1.1. Functional Forms of NS Mass Distributions

We consider three functional forms for the distributions of m_r and m_s. The first is a uniform distribution. It comprises two parameters, the lower and upper bound: $\left\{{m}_{r}^{l},{m}_{r}^{u}\right\}\in {{\rm{\Lambda }}}_{r}$, and $\left\{{m}_{s}^{l},{m}_{s}^{u}\right\}\in {{\rm{\Lambda }}}_{s}$.

The second form is a Gaussian distribution, which is commonly used in previous studies on the NS mass distributions (e.g., Özel et al. 2012; Kiziltan et al. 2013). In this case, the probability distribution is given by

$$\begin{eqnarray}&&\pi (m| {\rm{\Lambda }})=\pi (m| \left\{\mu ,\sigma \right\})=\displaystyle \frac{1}{\sigma \sqrt{2\pi }}\exp \left[-{\left(\displaystyle \frac{m-\mu }{\sqrt{2}\sigma }\right)}^{2}\right],\end{eqnarray} \tag{ 12 }$$

where μ and σ correspond to the mean and width of the distribution, respectively.

A natural extension to the Gaussian form is that the distribution has two distinct peaks, leading to the two-Gaussian distribution:

$$\begin{eqnarray}&&\begin{array}{l}\pi (m| {\rm{\Lambda }})=\pi (m| \left\{{\mu }_{1},{\sigma }_{1},{\mu }_{2},{\sigma }_{2},\alpha \right\})=\displaystyle \frac{\alpha }{{\sigma }_{1}\sqrt{2\pi }}\\ \quad \times \,\exp \left[-{\left(\displaystyle \frac{m-{\mu }_{1}}{\sqrt{2}{\sigma }_{1}}\right)}^{2}\right]+\displaystyle \frac{1-\alpha }{{\sigma }_{2}\sqrt{2\pi }}\exp \left[-{\left(\displaystyle \frac{m-{\mu }_{2}}{\sqrt{2}{\sigma }_{2}}\right)}^{2}\right].\end{array}\end{eqnarray} \tag{ 13 }$$

Here α and 1 − α give the weight of the first (with mean μ₁ and width σ₁) and second (with mean μ₂ and width σ₂) peak, respectively.

We adopt uniform priors for all parameters. Table 2 lists the prior ranges, which apply to both m_r and m_s. The Gaussian and two-Gaussian distributions are truncated at 0.8 and 2 M_⊙ to reflect these prior ranges. Note that the precisely measured NS masses are well within this range, so our results depend only weakly on the choice of this prior range.

Table 2.  Prior Ranges for Three Functional Forms of NS Mass Distribution

Form	Parameter	Min	Max
Uniform	ml	0.8	2
	mu	ml	2
Gaussian	μ	0.8	2
	σ	0.005	0.5
Two-Gaussian	μ1	0.8	2
	μ2	μ1	2
	σ1	0.005	0.5
	σ2	0.005	0.5
	α	0	1

Note. Mass parameters are in units of M_☉.

Download table as:  ASCIITypeset image

3.2. Model Selection

The second goal of this work is to differentiate hypotheses of NS mass distributions. Using Bayes' theorem, we can compute the ratio of posterior probability between the two hypotheses H₁ and H₂, as follows:

$$\begin{eqnarray}&&{O}=\displaystyle \frac{P({H}_{1}| \left\{{\rm{d}}\right\})}{P({H}_{2}| \left\{{\rm{d}}\right\})}=\displaystyle \frac{{{Z}}_{{\rm{\Lambda }}}({H}_{1})P({H}_{1})}{{{Z}}_{{\rm{\Lambda }}}({H}_{2})P({H}_{2})},\end{eqnarray} \tag{ 14 }$$

where P(H₁) and P(H₂) are the prior probability for hypotheses H₁ and H₂, respectively; ${{ \mathcal Z }}_{{\rm{\Lambda }}}({H}_{1})$ and ${{ \mathcal Z }}_{{\rm{\Lambda }}}({H}_{2})$ are the Bayesian evidence for hypotheses H₁ and H₂, respectively, see Equation (8). The quantity defined above is usually called the odds ratio. Assuming equal prior probability, P(H₁) = P(H₂), model selection is performed by computing the Bayes factor (BF)

$$\begin{eqnarray}&&{\mathrm{BF}}_{2}^{1}=\displaystyle \frac{{{ \mathcal Z }}_{{\rm{\Lambda }}}({H}_{1})}{{{ \mathcal Z }}_{{\rm{\Lambda }}}({H}_{2})}.\end{eqnarray} \tag{ 15 }$$

We follow the BF interpretation outlined by Kass & Raftery (1995): $0\lt {\mathrm{BF}}_{2}^{1}\lt 3$ indicates that the support for H₁ is "worth not more than a bare mention", $3\lt {\mathrm{BF}}_{2}^{1}\lt 20$ indicates positive support, $20\lt {\mathrm{BF}}_{2}^{1}\lt 150$ means the data strongly favor H₁, and ${\mathrm{BF}}_{2}^{1}\gt 150$ indicates very strong support. We use BF for small BFs and its natural logarithm, ln(BF), for large values. We choose a threshold for BF of 150, corresponding to ln(BF) = 5, as required for confident model selection.

Here we wish to compare two hypotheses: ${ \mathcal A }$—recycled and slow NSs follow an identical mass distribution, and ${ \mathcal B }$—they are drawn from two distinct populations. For hypothesis ${ \mathcal A }$, there are three possibilities: m_r and m_s follow (1) a uniform distribution, (2) a Gaussian distribution, or (3) a two-Gaussian distribution. In contrast, there are nine possibilities for hypothesis ${ \mathcal B }$, with m_r and m_s each following any one of the three distributions. The evidence for hypotheses ${ \mathcal A }$ and ${ \mathcal B }$ is

$$\begin{eqnarray}&&{{ \mathcal Z }}^{{ \mathcal A }}=\displaystyle \frac{1}{3}\left({{ \mathcal Z }}_{{\rm{u}}}^{{ \mathcal A }}+{{ \mathcal Z }}_{{\rm{g}}}^{{ \mathcal A }}+{{ \mathcal Z }}_{{\rm{t}}}^{{ \mathcal A }}\right),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{{ \mathcal Z }}^{{ \mathcal B }}=\displaystyle \frac{1}{9}\left({{ \mathcal Z }}_{\mathrm{uu}}^{{ \mathcal B }}+{{ \mathcal Z }}_{\mathrm{ug}}^{{ \mathcal B }}+{{ \mathcal Z }}_{\mathrm{ut}}^{{ \mathcal B }}+...+{{ \mathcal Z }}_{\mathrm{tt}}^{{ \mathcal B }}\right),\end{eqnarray} \tag{ 17 }$$

where u, g, and t in the subscripts denote uniform, Gaussian, and two-Gaussian distributions, respectively. Here the notation works as follows: ${{ \mathcal Z }}_{{\rm{u}}}^{{ \mathcal A }}$ is the evidence for the case that m_r and m_s are drawn independently from an identical uniform distribution, whereas ${{ \mathcal Z }}_{\mathrm{uu}}^{{ \mathcal B }}$ is the evidence for the case that m_r and m_s follow two uniform distributions with different parameters.

To quantitatively distinguish the support for each subhypothesis within the two hypotheses (three subhypotheses for ${ \mathcal A }$, and nine subhypotheses for ${ \mathcal B }$), we also compute the relative evidence (RE), defined as the ratio between the evidence of each subhypothesis and the total evidence of the hypothesis.⁵ For example, the RE for subhypothesis u under the hypothesis ${ \mathcal A }$ is

$$\begin{eqnarray}&&{\mathrm{RE}}_{{\rm{u}}}^{{ \mathcal A }}=\displaystyle \frac{{{ \mathcal Z }}_{{\rm{u}}}^{{ \mathcal A }}}{({{ \mathcal Z }}_{{\rm{u}}}^{{ \mathcal A }}+{{ \mathcal Z }}_{{\rm{g}}}^{{ \mathcal A }}+{{ \mathcal Z }}_{{\rm{t}}}^{{ \mathcal A }})}.\end{eqnarray} \tag{ 18 }$$

3.3. Posterior Predictive Distributions

Using posterior distributions of hyperparameters (Λ_r, Λ_s), we can derive a posterior predictive distribution (PPD)

$$\begin{eqnarray}&&{p}_{{{\rm{\Lambda }}}_{r}}({m}_{r})=\int d{{\rm{\Lambda }}}_{r}\,p({{\rm{\Lambda }}}_{r}| \left\{{\rm{d}}\right\})\pi ({m}_{r}| {{\rm{\Lambda }}}_{r}),\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{p}_{{{\rm{\Lambda }}}_{s}}({m}_{s})=\int d{{\rm{\Lambda }}}_{s}\,p({{\rm{\Lambda }}}_{s}| \left\{{\rm{d}}\right\})\pi ({m}_{s}| {{\rm{\Lambda }}}_{s}).\end{eqnarray} \tag{ 20 }$$

The PPD is our best guess for the updated prior on (m_r, m_s) given the data {d}. In some instances, it is useful to convert into a PPD for related variables such as M_T = m_r+m_s. This is accomplished through convolution

$$\begin{eqnarray}&&{p}_{{\rm{\Lambda }}}({M}_{T})=\int {{dm}}_{r}{p}_{{{\rm{\Lambda }}}_{r}}({m}_{r}){p}_{{{\rm{\Lambda }}}_{s}}({M}_{T}-{m}_{r}).\end{eqnarray} \tag{ 21 }$$

The PPD for mass ratio $q\in (0,1]$ can be obtained as follows:

$$\begin{eqnarray}\begin{array}{rcl}{p}_{{\rm{\Lambda }}}(q) & = & \displaystyle \frac{1}{2}\left[\displaystyle \int {{dm}}_{r}{p}_{{{\rm{\Lambda }}}_{s}}({{qm}}_{r}){p}_{{{\rm{\Lambda }}}_{r}}({m}_{r}){m}_{r}\right.\\ & & \left.+\displaystyle \int {{dm}}_{s}{p}_{{{\rm{\Lambda }}}_{r}}({{qm}}_{s}){p}_{{{\rm{\Lambda }}}_{s}}({m}_{s}){m}_{s}\right].\end{array}\end{eqnarray} \tag{ 22 }$$

4.1. Model Selection

We organize our results three ways. In Table 3, we present the BFs for all 17 Galactic DNS systems (merging and nonmerging). In Tables 4 and 5, we present the BFs for the 10 merging and 7 nonmerging systems, respectively. Here the merging (nonmerging) category includes those systems that will (not) merge within a Hubble time (see Table 1). We present hyperposteriors in the Appendix for the best subhypothesis found in Table 3. We highlight several take-away messages:

• 1.  
  The data exhibits positive support for the hypothesis that the recycled NS mass and the slow NS mass are distributed differently. Using the entire sample, we find BF = 9 for hypothesis ${ \mathcal B }$ over ${ \mathcal A }$. The best subhypothesis (${{ \mathcal Z }}_{\mathrm{tu}}^{{ \mathcal B }}$) is a two-Gaussian distribution for m_r and a uniform distribution for m_s, with an RE of 64%. We find it is strongly favored (BF = 29) over the subhypothesis ${{ \mathcal Z }}_{{\rm{g}}}^{{ \mathcal A }}$ (the conventional-wisdom model), which states that m_r and m_s are drawn from an identical Gaussian distribution. The second and third best subhypotheses are ${{ \mathcal Z }}_{\mathrm{uu}}^{{ \mathcal B }}$ and ${{ \mathcal Z }}_{\mathrm{gu}}^{{ \mathcal B }}$, having an RE of 20% and 9%, respectively.
• 2.  
  Some subhypotheses are strongly ruled out. For example, the combination of Gaussian m_r and two-Gaussian m_s, or uniform m_r and two-Gaussian m_s, yield an RE of ∼0.1%.
• 3.  
  The results do not change significantly if we analyze merging binaries as a special group. This is unsurprising because 9 out of 12 systems with precise component mass measurements are in this category.
• 4.  
  Data from nonmerging binaries is inconclusive. Presently, the BFs in the nonmerging table are too small to draw even preliminary conclusions. However, future observations could potentially highlight differences in the merging and nonmerging populations, for example, due to differences in the accretion process (Tauris et al. 2017).

Table 3.  For All Galactic DNS Systems: Relative Evidence (RE) for Subhypotheses in Hypotheses ${ \mathcal A }$ and ${ \mathcal B }$, and the BF between the Best Subhypothesis in ${ \mathcal A }$ and the Best Subhypothesis in ${ \mathcal B }$

Subhypothesis		${{ \mathcal Z }}_{{\rm{u}}}^{{ \mathcal A }}$	${{ \mathcal Z }}_{{\rm{g}}}^{{ \mathcal A }}$	${{ \mathcal Z }}_{{\rm{t}}}^{{ \mathcal A }}$
	RE	0.027	0.597	0.376
${{ \mathcal Z }}_{\mathrm{uu}}^{{ \mathcal B }}$	0.197	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{ug}}^{{ \mathcal B }}$	0.011	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{ut}}^{{ \mathcal B }}$	0.002	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{gu}}^{{ \mathcal B }}$	0.087	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{gg}}^{{ \mathcal B }}$	0.004	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{gt}}^{{ \mathcal B }}$	0.001	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{tu}}^{{ \mathcal B }}$	0.643	⋯	28.6	⋯
${{ \mathcal Z }}_{\mathrm{tg}}^{{ \mathcal B }}$	0.040	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{tt}}^{{ \mathcal B }}$	0.015	⋯	⋯	⋯

Note. On average, we find ${{ \mathcal Z }}^{{ \mathcal B }}/{{ \mathcal Z }}^{{ \mathcal A }}=8.9$.

Download table as:  ASCIITypeset image

Table 4.  As Table 3, but for Merging Galactic DNS Systems

Subhypothesis	RE	${{ \mathcal Z }}_{{\rm{u}}}^{{ \mathcal A }}$	${{ \mathcal Z }}_{{\rm{g}}}^{{ \mathcal A }}$	${{ \mathcal Z }}_{{\rm{t}}}^{{ \mathcal A }}$
		0.087	0.486	0.427
${{ \mathcal Z }}_{\mathrm{uu}}^{{ \mathcal B }}$	0.155	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{ug}}^{{ \mathcal B }}$	0.008	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{ut}}^{{ \mathcal B }}$	0.002	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{gu}}^{{ \mathcal B }}$	0.041	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{gg}}^{{ \mathcal B }}$	0.002	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{gt}}^{{ \mathcal B }}$	0.001	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{tu}}^{{ \mathcal B }}$	0.672	⋯	13.9	⋯
${{ \mathcal Z }}_{\mathrm{tg}}^{{ \mathcal B }}$	0.086	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{tt}}^{{ \mathcal B }}$	0.034	⋯	⋯	⋯

Note. In this case, we find ${{ \mathcal Z }}^{{ \mathcal B }}/{{ \mathcal Z }}^{{ \mathcal A }}=3.3.$

Download table as:  ASCIITypeset image

Table 5.  As Table 3, but for Nonmerging Galactic DNS Systems.

Subhypothesis	RE	${{ \mathcal Z }}_{{\rm{u}}}^{{ \mathcal A }}$	${{ \mathcal Z }}_{{\rm{g}}}^{{ \mathcal A }}$	${{ \mathcal Z }}_{{\rm{t}}}^{{ \mathcal A }}$
		0.705	0.204	0.091
${{ \mathcal Z }}_{\mathrm{uu}}^{{ \mathcal B }}$	0.667	4.8	⋯	⋯
${{ \mathcal Z }}_{\mathrm{ug}}^{{ \mathcal B }}$	0.089	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{ut}}^{{ \mathcal B }}$	0.041	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{gu}}^{{ \mathcal B }}$	0.092	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{gg}}^{{ \mathcal B }}$	0.013	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{gt}}^{{ \mathcal B }}$	0.007	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{tu}}^{{ \mathcal B }}$	0.075	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{tg}}^{{ \mathcal B }}$	0.011	⋯	⋯	⋯
${{ \mathcal Z }}_{\mathrm{tt}}^{{ \mathcal B }}$	0.005	⋯	⋯	⋯

Note. In this case, we find ${{ \mathcal Z }}^{{ \mathcal B }}/{{ \mathcal Z }}^{{ \mathcal A }}=1.7.$

Download table as:  ASCIITypeset image

4.1.1. Discussion

In regard to the difference in the distribution of m_r and m_s, Özel et al. (2012) carried out a fit of the Gaussian distribution separately to the masses of pulsars and companions. They found the most likely parameters are μ = 1.35 M_☉, σ = 0.05 M_☉, and μ = 1.32 M_☉, σ = 0.05 M_☉, for pulsars and companions, respectively. Therefore, they concluded that two members of DNS systems are drawn from the same underlying population. Although PSR J1906+0746, which is a slow pulsar, was mixed up with other recycled pulsars in Özel et al. (2012), the difference between our result and theirs is due to the addition of new DNS systems with masses measured outside that "typical" range. Some notable examples include the most massive m_r at 1.64 M_☉ (PSR J1913+1102), the lightest m_s at 1.174 M_☉ (PSR J0453+1559), and the most massive m_s at 1.46 M_☉ (PSR J0509+3801). Tauris et al. (2017) noticed a ∼0.1 M_☉ difference in the histogram of m_r and m_s. They argued that this discrepancy is due to differences in the distribution of birth masses because recycled NSs in DNS systems are expected to gain ≤0.02 M_☉ from accretion.

The median values of marginalized posterior distributions under the best subhypothesis ${{ \mathcal Z }}_{\mathrm{tu}}^{{ \mathcal B }}$ are {μ₁, σ₁, μ₂, σ₂, α} = {1.34 M_☉, 0.02 M_☉, 1.47 M_☉, 0.15 M_☉, 0.68} for m_r, and ${m}_{s}^{l}=1.16{M}_{\odot }$ and ${m}_{s}^{u}=1.42{M}_{\odot }$. We note that the secondary peak of the two-Gaussian distribution for m_r is consistent with the distribution derived for recycled NSs with mostly white dwarf companions, μ = 1.46 M_☉ and σ = 0.21 M_☉, in Özel et al. (2012). If this bimodal structure is confirmed with the addition of more DNS measurements, it could lend support to the claim by Schwab et al. (2010) that the NS mass spectrum has two peaks: one from Fe core collapse supernovae and one from electron capture supernovae.⁶ Extending our analysis to all NS mass measurements may shed new light on this problem, as well as provide insights into the massive star evolution (Raithel et al. 2018).

Recently, Huang et al. (2018) analyzed the total mass measurements of 15 Galactic DNS systems, and found strong preference for a two-Gaussian distribution against the single-Gaussian model, using a simple likelihood ratio statistic. Keitel (2019) reanalyzed the same data and questioned the validity of such a claim. By performing various statistical tests, Keitel (2019) suggested that there is no clear support for the two-Gaussian model and more data are required to draw firm conclusions. Since we focus on the component mass measurements in this work, our results are not directly comparable to that of Keitel (2019). A common message from our work and Keitel (2019) is that additional DNS discoveries are needed to settle the question of subpopulations.

4.1.2. How Many Observations Are Required to Distinguish Two Populations?

Here we determine how many observations are required to confidently distinguish the two distributions for m_r and m_s, assuming they are indeed drawn from two different parent populations. We generate mock mass measurements of DNS systems assuming that m_r follows a two-Gaussian distribution, and m_s follows a uniform distribution, with parameters taking the median values of marginalized posterior distributions.

The mass measurement uncertainties of DNS systems to be discovered in future pulsar surveys are difficult to predict. Uncertainties of existing measurements vary by orders of magnitude, as can be seen in Figure 1. For simplicity and conservativeness, we assume a measurement error of 0.04 M_☉ for both m_r and m_s, and a smaller error of 0.014 M_☉ for M_T, similar to those of PSR J1913+1102, for all future discoveries. We further assume the posterior probability $p({m}_{r},{m}_{s}| d)$ peaks at injected values of m_r and m_s for all the mock measurements. In reality, measurements do not always produce posteriors that peak at the true parameter values. This would make our calculations somewhat optimistic by a factor of ≲2. Since the assumed mass measurement precision is conservative compared to those in Table 1, we expect our results to be realistic when interpreted as the required number of DNS systems with at least two post-Keplerian parameters measured. What is needed is the discovery of new DNS systems and measurement of their component masses through follow-up radio timing observations.

Figure 3 shows the natural logarithm of BF as a function of the number of DNS discoveries. The solid blue curve with shaded region shows the BF between the hypotheses ${ \mathcal B }$ and ${ \mathcal A }$. The solid black curve with error bars shows the BF between the true subhypothesis (two-Gaussian m_r, and uniform m_s), and the second best subhypothesis, or the best subhypothesis when the true subhypothesis does not produce the highest evidence,⁷ within hypothesis ${ \mathcal B }$. The red horizontal line marks the threshold of $\mathrm{ln}(\mathrm{BF})=5$ as required for confident model selection. One can see that, on average, ≈20 measurements are required to tell that m_r and m_s follow different distributions, and about 60 measurements are required to confidently distinguish the true subhypothesis from the other eight subhypotheses. The blue and black stars show results from existing observations, which are in good agreement with the simulations.

Zoom In Zoom Out Reset image size

4.2. Posterior Predictive Distributions

In this section, we present PPDs for the subhypothesis that best describes masses of all 17 Galactic DNS systems (two-Gaussian m_r and uniform m_s). Figure 4 shows the PPDs for m_r and m_s. For comparison, we also plot the conventional Gaussian distribution with a mean of 1.33 M_☉ and a width of 0.09 M_☉.

Zoom In Zoom Out Reset image size

Distributions shown in Figure 4 can be converted to that for total mass and mass ratio, both of which are important for gravitational-wave data analysis. The distribution of total mass is also a key to telling how often NS mergers produce black holes versus hypermassive NSs. Figure 5 shows the PPDs (lower panels) and their cumulative distribution (upper panels) for M_T and q. The red vertical lines in the upper panel mark the 99% and 90% confidence lower bounds. For example, we find that M_T is above 2.334 M_☉ and q is larger than 0.691 with 99% confidence. In contrast, these bounds become 2.364 M_☉ and 0.779 for M_T and q, respectively, if m_r and m_s are assumed to follow the conventional Gaussian distribution with μ = 1.33 M_☉ and σ = 0.09 M_☉.

Zoom In Zoom Out Reset image size

The number of known DNS systems in our Galaxy has nearly doubled in the past several years. Radio pulsar timing observations provide precise measurements of NS masses for many of these systems. This growing sample allows tighter constraints to be placed on the mass distribution. In this work, we make use of mass measurements of 17 Galactic DNS systems, including precise component mass measurements of 12 systems, and perform Bayesian inference and model selection on the mass distribution. We find tentative evidence for two distinct mass distributions for recycled and slow NSs, and for the bimodality of the recycled NS mass distribution. Furthermore, the conventional model which states that both component NSs follow an identical Gaussian distribution is disfavored with a BF of 29, compared to a two-Gaussian mass distribution for recycled NSs and uniform mass distribution for slow NSs. If recycled NSs have a different mass distribution from slow NSs, this could hint at subtleties in our understanding of supernovae mechanisms and/or accretion processes.

We show that precise measurements of component masses of around 20 DNS systems are required to confidently tell—with ln(BF) > 5—whether or not recycled and slow NSs follow different mass distributions. To establish the detailed shape of the mass distribution, measurements of 60 binaries are needed. This might become possible in the next decade, through ongoing pulsar surveys and new surveys with the Five-hundred-meter Aperture Spherical Telescope (Nan et al. 2011) in China and MeerKAT (Bailes et al. 2018), a precursor for the planned Square Kilometre Array⁸ in South Africa. The method developed here will prove useful for further review of the mass distribution as new DNS systems are discovered and their masses are measured.

In this work we have focused on radio pulsar observations. Our analysis does not take into account selection effects, which means the distributions derived here do not necessarily match the mass distribution of DNSs at birth or merger. The former can be predicted from population synthesis studies. The latter is directly measurable in gravitational-wave observations and will soon be probed with the many new discoveries expected when Advanced LIGO and Advanced Virgo progress toward their design sensitivities in the next few years. We plan to develop a framework that links the three populations, facilitating comparison of gravitational-wave mergers with the Galactic DNS population.

X.Z. and E.T. are supported by ARC CE170100004. E.T. is additionally supported through ARC FT150100281. We thank Matthew Bailes for helpful feedback about selection effects. We also thank the anonymous referee, David Keitel, and Harald Pfeiffer for useful comments on the manuscript.

Here we present the (hyper)posterior distributions for the best subhypothesis found in Section 4.1. Figures 6 and 7 show the posterior distributions of parameters for the two-Gaussian distribution for recycled NS mass m_r, and the uniform distribution for nonrecycled NS mass m_s, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size