The mass distribution of Galactic double neutron stars

TL;DR

This study revisits the mass distribution of Galactic double neutron stars using Bayesian inference with an expanded sample, suggesting possible differences between recycled and non-recycled neutron star masses and exploring various distribution models.

Contribution

It introduces a Bayesian framework with expanded models to analyze neutron star mass distributions, considering separate distributions for recycled and non-recycled stars.

Findings

Support for different mass distributions for recycled and non-recycled neutron stars (Bayes factor 9).

Most probable model involves a two-Gaussian distribution for recycled stars and a uniform distribution for non-recycled stars.

Approximately 20 precise measurements are needed to confidently distinguish distribution differences.

Abstract

The conventional wisdom, dating back to 2012, is that the mass distribution of Galactic double neutron stars is well-fit by a Gaussian distribution with a mean of $1.33 M_\odot$ and a width of $0.09 M_\odot$. With the recent discovery of new Galactic double neutron stars and GW170817, the first neutron star merger event to be observed with gravitational waves, it is timely to revisit this model. In order to constrain the mass distribution of double neutron stars, we perform Bayesian inference using a sample of 17 Galactic double neutron stars effectively doubling the sample used in previous studies. We expand the space of models so that the recycled neutron star need not be drawn from the same distribution as the non-recycled companion. Moreover, we consider different functional forms including uniform, single-Gaussian, and two-Gaussian distributions. While there is insufficient data to draw firm conclusions, we find positive support (a Bayes factor of 9) for the hypothesis that recycled and non-recycled neutron stars have distinct mass distributions. The most probable model---preferred with a Bayes factor of 29 over the conventional model---is one in which the recycled neutron star mass is distributed according to a two-Gaussian distribution and the non-recycled neutron star mass is distributed uniformly. We show that precise component mass measurements of $\approx 20$ double neutron stars are required in order to determine with high confidence (a Bayes factor of 150) if recycled and non-recycled neutron stars come from a common distribution. Approximately $60$ are needed in order to establish the detailed shape of the distributions.

Tables

Table 1. Table 1: Mass measurements of Galactic DNS systems: the binary total mass ( M T subscript 𝑀 𝑇 M_{T} ), the masses of the recycled NS ( m r subscript 𝑚 𝑟 m_{r} ) and the slow NS ( m s subscript 𝑚 𝑠 m_{s} ), binary chirp mass ( ℳ c subscript ℳ 𝑐 \mathcal{M}_{c} ), mass ratio ( q 𝑞 q ), binary orbital period ( P b subscript 𝑃 𝑏 P_{b} ) and coalescence time T c subscript 𝑇 𝑐 T_{c} . Figures in parentheses are 1- σ 𝜎 \sigma 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 subscript 𝑚 𝑟 m_{r} and m s subscript 𝑚 𝑠 m_{s} , derived from their measured M T subscript 𝑀 𝑇 M_{T} and mass functions: 0.268184 (J1946+2052), 0.1223898 (J1411+2551), 0.128121 (J1811 − - 1736), 0.29413 (J1829+2456), 0.34690765 (J1930 − - 1852). For ℳ c subscript ℳ 𝑐 \mathcal{M}_{c} and for q 𝑞 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 ( 2017 ) ; (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 subscript 𝑚 𝑟 m_{r} and m s subscript 𝑚 𝑠 m_{s} values of PSR J1518+4904 were taken from Tauris et al. ( 2017 ) who cites a private communication with G. Janssen.

Pulsar Name	$M_T$ ($M_{\odot }$)	$m_r$ ($M_{\odot }$)	$m_s$ ($M_{\odot }$)	${ℳ}_c$ ($M_{\odot }$)	$q$	$P_b$ (day)	$T_c$ (Gyr)	References
Systems will merge within a Hubble time
J1946+2052	2.50(4)		$>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.886(1)	1.65(5)	1.24(5)	1.242(8)	0.75(5)	0.206	0.473	(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)		$>0.92$	(1.05, 1.11)	(0.57, 0.95)	2.616	466	(14)
J1811$-$1736	2.57(10)		$>0.91$	(1.02, 1.17)	(0.58, 0.95)	18.78	1794	(15)
J1829+2456	2.59(2)		$>1.25$	(1.08, 1.14)	(0.65, 1)	1.176	55	(16)
J1930$-$1852	2.59(4)		$>1.30$	(1.07, 1.15)	(0.58, 0.96)	45.06	$\sim {10}^5$	(17)

Table 2. Table 2: Prior ranges for three functional forms of NS mass distribution. Mass parameters are in unit of M ☉ subscript 𝑀 ☉ M_{\sun} .

Form	Parameter	Min	Max
Uniform	$m^l$	0.8	2
	$m^u$	$m^l$	2
Gaussian	$\mu$	0.8	2
	$\sigma$	0.005	0.5
Two-Gaussian	${\mu }_1$	0.8	2
	${\mu }_2$	${\mu }_1$	2
	${\sigma }_1$	0.005	0.5
	${\sigma }_2$	0.005	0.5
	$\alpha$	0	1

Table 3. Table 3: For all Galactic DNS systems: relative evidences (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} . On average, we find 𝒵 ℬ / 𝒵 𝒜 = 8.9 superscript 𝒵 ℬ superscript 𝒵 𝒜 8.9 \mathcal{Z^{B}}/\mathcal{Z^{A}}=8.9 .

Subhypothesis		${𝒵}_{\text{u}}^{𝒜}$	${𝒵}_{\text{g}}^{𝒜}$	${𝒵}_{\text{t}}^{𝒜}$
	RE	0.027	0.597	0.376
${𝒵}_{\text{uu}}^{ℬ}$	0.197
${𝒵}_{\text{ug}}^{ℬ}$	0.011
${𝒵}_{\text{ut}}^{ℬ}$	0.002
${𝒵}_{\text{gu}}^{ℬ}$	0.087
${𝒵}_{\text{gg}}^{ℬ}$	0.004
${𝒵}_{\text{gt}}^{ℬ}$	0.001
${𝒵}_{\text{tu}}^{ℬ}$	0.643		28.6
${𝒵}_{\text{tg}}^{ℬ}$	0.040
${𝒵}_{\text{tt}}^{ℬ}$	0.015

Table 4. Table 4: As Table 3 but for Merging Galactic DNS systems. In this case, we find 𝒵 ℬ / 𝒵 𝒜 = 3.3 superscript 𝒵 ℬ superscript 𝒵 𝒜 3.3 \mathcal{Z^{B}}/\mathcal{Z^{A}}=3.3 .

Subhypothesis		${𝒵}_{\text{u}}^{𝒜}$	${𝒵}_{\text{g}}^{𝒜}$	${𝒵}_{\text{t}}^{𝒜}$
	RE	0.087	0.486	0.427
${𝒵}_{\text{uu}}^{ℬ}$	0.155
${𝒵}_{\text{ug}}^{ℬ}$	0.008
${𝒵}_{\text{ut}}^{ℬ}$	0.002
${𝒵}_{\text{gu}}^{ℬ}$	0.041
${𝒵}_{\text{gg}}^{ℬ}$	0.002
${𝒵}_{\text{gt}}^{ℬ}$	0.001
${𝒵}_{\text{tu}}^{ℬ}$	0.672		13.9
${𝒵}_{\text{tg}}^{ℬ}$	0.086
${𝒵}_{\text{tt}}^{ℬ}$	0.034

Table 5. Table 5: As Table 3 but for Nonmerging Galactic DNS systems. In this case, we find 𝒵 ℬ / 𝒵 𝒜 = 1.7 superscript 𝒵 ℬ superscript 𝒵 𝒜 1.7 \mathcal{Z^{B}}/\mathcal{Z^{A}}=1.7 .

Subhypothesis		${𝒵}_{\text{u}}^{𝒜}$	${𝒵}_{\text{g}}^{𝒜}$	${𝒵}_{\text{t}}^{𝒜}$
	RE	0.705	0.204	0.091
${𝒵}_{\text{uu}}^{ℬ}$	0.667	4.8
${𝒵}_{\text{ug}}^{ℬ}$	0.089
${𝒵}_{\text{ut}}^{ℬ}$	0.041
${𝒵}_{\text{gu}}^{ℬ}$	0.092
${𝒵}_{\text{gg}}^{ℬ}$	0.013
${𝒵}_{\text{gt}}^{ℬ}$	0.007
${𝒵}_{\text{tu}}^{ℬ}$	0.075
${𝒵}_{\text{tg}}^{ℬ}$	0.011
${𝒵}_{\text{tt}}^{ℬ}$	0.005

Equations

$${\mathcal{M}}_c=\frac{(m_1{m}_2{)}^{3/5}}{(m_1+m_2{)}^{1/5}},$$

$$\begin{array}{rl}p(m_r,m_s)& \propto \exp \left[-\frac{(m_r+m_s-{\hat{M}}_T{)}^2}{2{\sigma }_T^2}\right]\\ & \times {\int }_{-1}^1\delta (f-\hat{f})\text{d}(\cos i)\\ & \propto \frac{\exp \left[-\frac{(m_r+m_s-{\hat{M}}_T{)}^2}{2{\sigma }_T^2}\right]}{\frac{3m_s^3}{(m_r+m_s{)}^2}{\left[\frac{\hat{f}(m_r+m_s{)}^2}{m_s^3}\right]}^{\frac{1}{3}}{\left[1-\frac{{\hat{f}}^{2/3}(m_r+m_s{)}^{4/3}}{m_s^2}\right]}^{\frac{1}{2}}},\end{array}$$

$$f=\frac{(m_c\sin i{)}^3}{M_T^2}={\left(\frac{2\pi }{P_b}\right)}^2\frac{(a\sin i{)}^3}{G},$$

$$p(m_r^i,m_s^i|d^i),$$

$$\pi (m_r^i|{\Lambda }_r,H)$$

$$\pi (m_s^i|{\Lambda }_s,H).$$

$$p(\Lambda |\{d\})=$$

$$\pi (m_r^i|{\Lambda }_r)\pi (m_s^i|{\Lambda }_s)\pi ({\Lambda }_r)\pi ({\Lambda }_s).$$

$${\mathcal{Z}}_{\Lambda }\equiv$$

$$\mathcal{L}(d^i|m_r^i,m_s^i)\pi (m_r^i|{\Lambda }_r)\pi (m_s^i|{\Lambda }_s)\pi ({\Lambda }_r)\pi ({\Lambda }_s).$$

$$\mathcal{L}(d^i|m_r^i,m_s^i)=\frac{{\mathcal{Z}}_0^i}{{\pi }_0(m_r^i,m_s^i)}p(m_r^i,m_s^i|d^i),$$

$${\mathcal{Z}}_0\equiv \int d{m}_r\int d{m}_s\mathcal{L}(d|m_r,m_s){\pi }_0(m_r,m_s),$$

$$p(\Lambda |\{d\})=$$

$$\frac{\pi (m_r^i|{\Lambda }_r)\pi (m_s^i|{\Lambda }_s)}{{\pi }_0(m_r^i,m_s^i)}\pi ({\Lambda }_r)\pi ({\Lambda }_s).$$

$$\pi (m|\Lambda )=\pi (m|\{\mu ,\sigma \})=\frac{1}{\sigma \sqrt{2\pi }}\exp \left[-{\left(\frac{m-\mu }{\sqrt{2}\sigma }\right)}^2\right],$$

$$\pi (m|\Lambda )$$

$$\mathcal{O}=\frac{P(H_1|\{d\})}{P(H_2|\{d\})}=\frac{{\mathcal{Z}}_{\Lambda }(H_1)P(H_1)}{{\mathcal{Z}}_{\Lambda }(H_2)P(H_2)},$$

$${\text{BF}}_2^1=\frac{{\mathcal{Z}}_{\Lambda }(H_1)}{{\mathcal{Z}}_{\Lambda }(H_2)}$$

$${\mathcal{Z}}^{\mathcal{A}}=\frac{1}{3}\left({\mathcal{Z}}_{\text{u}}^{\mathcal{A}}+{\mathcal{Z}}_{\text{g}}^{\mathcal{A}}+{\mathcal{Z}}_{\text{t}}^{\mathcal{A}}\right),$$

$${\mathcal{Z}}^{\mathcal{B}}=\frac{1}{9}\left({\mathcal{Z}}_{\text{uu}}^{\mathcal{B}}+{\mathcal{Z}}_{\text{ug}}^{\mathcal{B}}+{\mathcal{Z}}_{\text{ut}}^{\mathcal{B}}+\cdots +{\mathcal{Z}}_{\text{tt}}^{\mathcal{B}}\right),$$

$${\text{RE}}_{\text{u}}^{\mathcal{A}}=\frac{{\mathcal{Z}}_{\text{u}}^{\mathcal{A}}}{({\mathcal{Z}}_{\text{u}}^{\mathcal{A}}+{\mathcal{Z}}_{\text{g}}^{\mathcal{A}}+{\mathcal{Z}}_{\text{t}}^{\mathcal{A}})}.$$

$$p_{{\Lambda }_r}(m_r)=$$

$$p_{{\Lambda }_s}(m_s)=$$

$$p_{\Lambda }(M_T)=\int \text{d}{m}_r{p}_{{\Lambda }_r}(m_r)p_{{\Lambda }_s}(M_T-m_r).$$

$$p_{\Lambda }(q)=$$

$$+\int \text{d}{m}_s{p}_{{\Lambda }_r}(q{m}_s)p_{{\Lambda }_s}(m_s)m_s].$$

Full text

pulsars: general — stars: neutron — methods: data analysis — gravitational waves

1 Introduction

The breakthrough discovery of GW170817 was the first ever detection of gravitational waves arising from a NS merger, $44\,\mathrm{Mpc}$ from Earth in the galaxy NGC4993 (Abbott et al., 2017; Coulter et al., 2017). More binary neutron star 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 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 $\mu=1.33M_{\odot}$ and $\sigma=0.09M_{\odot}$.

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 datasets, 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.

2 Masses of Galactic DNS

Among 18 Galactic DNS systems, 17 have reported measurements of the binary total mass, and 12 of which have masses measured for both component stars; see Table III 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 $\sim 10-100$ ms through an accretion/recycling process, whereas the second born NS quickly ($\sim$ 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.

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

[TABLE]

and mass ratio $q=m_{2}/m_{1}\leq 1$.

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

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 distribution111This is a reasonable approximation to the original posterior distribution derived from radio pulsar timing observations; see, e.g., Figure 2 in Lynch et al. (2018). with a covariance of $(\sigma_{T}^{2}-\sigma_{r}^{2}-\sigma_{s}^{2})/2$, where $\sigma_{T}$, $\sigma_{r}$ and $\sigma_{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

[TABLE]

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

[TABLE]

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 Bayesian Inference

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 observations222So far there is only one DNS system (GW170817) that is identified via gravitational waves. Its mass constraints are presented for $m_{1}$ and $m_{2}$. Since we are interested in the distribution of $m_{r}$ and $m_{s}$, we focus on radio pulsar observations in this work.. These previously published constraints on NS masses are posterior distributions, which we denote as

[TABLE]

where $i$ refers to the $i$th 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 $(\Lambda_{r},\Lambda_{s})\in\Lambda$

[TABLE]

The hyperparameters $(\Lambda_{r},\Lambda_{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}|\Lambda_{r},H)$ is a Gaussian distribution in $m_{r}$. In this case, $\Lambda_{r}$ consists of hyperparameters for the mean $\mu_{r}\in\Lambda_{r}$ and width $\sigma_{r}\in\Lambda_{r}$ of the distribution.

Our first goal is to calculate (hyper)posterior distributions for $(\Lambda_{r},\Lambda_{s})$. Given a data set with $N$ binary mass measurements, the hyperposterior distribution is given by (Thrane & Talbot, 2019)

[TABLE]

Here $\pi(\Lambda_{r})$ and $\pi(\Lambda_{s})$ are the priors of our hyperparameters, which we take to be flat. The term ${\cal Z}_{\Lambda}$ is the (hyper)evidence

[TABLE]

The variable ${\cal L}(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

[TABLE]

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

[TABLE]

is the initial evidence. Without access to the raw radio data $d$, we do not know ${\cal 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 ${\cal Z}_{0}^{i}=1$ henceforth.

Plugging Equation 9 into 3.1, we obtain

[TABLE]

Using Eq. 3.1, 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: $\{m_{r}^{l},m_{r}^{u}\}\in\Lambda_{r}$, and $\{m_{s}^{l},m_{s}^{u}\}\in\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:

[TABLE]

where $\mu$ and $\sigma$ 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:

[TABLE]

Here $\alpha$ and $1-\alpha$ give the weight of the first (with mean $\mu_{1}$ and width $\sigma_{1}$) and second (with mean $\mu_{2}$ and width $\sigma_{2}$) 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_{\odot}$ 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.

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_{1}$ and $H_{2}$ as follows:

[TABLE]

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

[TABLE]

We follow the BF interpretation outlined by Kass & Raftery (1995): $0<\text{BF}^{1}_{2}<3$ indicates that the support for $H_{1}$ is ‘worth not more than a bare mention’, $3<\text{BF}^{1}_{2}<20$ indicates positive support, $20<\text{BF}^{1}_{2}<150$ means the data strongly favor $H_{1}$, and $\text{BF}^{1}_{2}>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 a 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, or 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

[TABLE]

where u, g, and t in the subscripts denote uniform, Gaussian, and two-Gaussian distributions, respectively. Here the notation works as follows: $\mathcal{Z}^{\mathcal{A}}_{\text{u}}$ is the evidence for the case that $m_{r}$ and $m_{s}$ are drawn independently from an identical uniform distribution, whereas $\mathcal{Z}^{\mathcal{B}}_{\text{uu}}$ 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 hypothesis333This is equivalent to the relative model probability used in Farr et al. (2011), who performed model selection for the mass distribution of stellar-mass black holes.. For example, the RE for subhypothesis u under the hypothesis $\mathcal{A}$ is

[TABLE]

3.3 Posterior Predictive Distributions

Using posterior distributions of hyperparameters $(\Lambda_{r},\Lambda_{s})$, we can derive a posterior predictive distribution (PPD)

[TABLE]

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

[TABLE]

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

[TABLE]

4 Results

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:

• •

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 $\text{BF}=9$ for hypothesis $\mathcal{B}$ over $\mathcal{A}$. The best subhypothesis ($\mathcal{Z}^{\mathcal{B}}_{\text{tu}}$) 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}^{\mathcal{A}}_{\text{g}}$ (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}^{\mathcal{B}}_{\text{uu}}$ and $\mathcal{Z}^{\mathcal{B}}_{\text{gu}}$, having an RE of 20% and 9%, respectively.

• •

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 $\sim 0.1\%$.

• •

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.

• •

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

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 $\mu=1.35M_{\sun}$, $\sigma=0.05M_{\sun}$, and $\mu=1.32M_{\sun}$, $\sigma=0.05M_{\sun}$, 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_{\sun}$ (PSR J1913+1102), the lightest $m_{s}$ at $1.174M_{\sun}$ (PSR J0453+1559), and the most massive $m_{s}$ at $1.46M_{\sun}$ (PSR J0509+3801). Tauris et al. (2017) noticed a $\sim 0.1M_{\sun}$ 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 $\leq 0.02M_{\sun}$ from accretion.

The median values of marginalized posterior distributions under the best subhypothesis $\mathcal{Z}^{\mathcal{B}}_{\text{tu}}$ are $\{\mu_{1},\sigma_{1},\mu_{2},\sigma_{2},\alpha\}$ $=\{1.34M_{\sun},0.02M_{\sun},1.47M_{\sun},0.15M_{\sun},0.68\}$ for $m_{r}$, and $m_{s}^{l}=1.16M_{\sun}$ and $m_{s}^{u}=1.42M_{\sun}$. 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, $\mu=1.46M_{\sun}$ and $\sigma=0.21M_{\sun}$, 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 supernovae444We caution that the mean values of our peaks appear to differ from those in Schwab et al. (2010).. 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) re-analyzed 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_{\sun}$ for both $m_{r}$ and $m_{s}$, and a smaller error of 0.014 $M_{\sun}$ 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 $\lesssim 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 evidence555This is usually the case when there are fewer than 20 measurements., within hypothesis $\mathcal{B}$. The red horizontal line marks the threshold of $\ln{(\text{BF})}=5$ as required for confident model selection. One can see that, on average, $\approx 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.

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_{\sun}$ and a width of 0.09 $M_{\sun}$.

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 do 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_{\sun}$ and $q$ is larger than 0.691 with 99% confidence. In contrast, these bounds become 2.364 $M_{\sun}$ and 0.779 for $M_{T}$ and $q$, respectively, if $m_{r}$ and $m_{s}$ are assumed to follow the conventional Gaussian distribution with $\mu=1.33M_{\sun}$ and $\sigma=0.09M_{\sun}$.

5 Conclusions

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 Array666https://www.skatelescope.org/ 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.