THE NEUTRON STAR MASS DISTRIBUTION

The canonical mass limit M_ch ∼ 1.4 M_☉ is the critical mass beyond which the degenerate remnant core of a massive star or a WD will lose gravitational stability and collapse into an NS. This limiting mass is an approximation that is sensitive to several nuclear, relativistic, and geometric effects (see Ghosh 2007; Haensel et al. 2007, for review). In addition to these effects, the variety of evolutionary processes that produce NSs warrant a careful treatment.

A more precise parameterization of the Chandrasekhar mass is

$$\begin{equation} M_{{\rm ch}}=5.83\, Y_{e}^{2}\, {M}_{\odot } , \end{equation} \tag{ 1 }$$

where Y_e = n_p/(n_p + n_n) is the electron fraction. A perfect neutron–proton equality (n_p = n_n) with Y_e = 0.50 yields a critical mass of

$$\begin{equation} M_{{\rm ch}}=1.457\, {M}_{\odot }. \end{equation} \tag{ 2 }$$

However, we have a sufficiently good insight into the processes that affect M_ch. So, we can go beyond the idealized cases and estimate the remnant's expected initial mass more realistically.

The inclusion of more reasonable electron fractions (Y_e < 0.50) yields smaller values for M_ch. General relativistic implications, surface boundary pressure corrections, and the reduction of pressure due to non-ideal Coulomb interactions (e⁻–e⁻ repulsion, ion–ion repulsion, and e⁻–ion attraction) at high densities all reduce the upper limit of M_ch.

On the other hand, the electrons of the progenitor material (i.e., the WD or the core of a massive star) are not completely relativistic. This reduces the pressure leading to an increase in the amount of mass required to reach the gravitational potential to collapse the star. Finite entropy corrections and the effects of rotation will also enhance the stability for additional mass. These corrections, as a result, yield a higher upper limit for M_ch.

The level of impact on the birth masses due to some of these competing effects is not well constrained as the details of the processes are not well understood. An inclusion of the effects that are due to the diversity in the evolutionary processes alone requires a ≈20% correction (for a detailed numerical treatment, see Butterworth & Ipser 1975) and therefore implies a broader mass range, i.e., M_ch ∼ 1.17–1.75 M_☉.

The measured masses, however, are the effective gravitational masses rather than a measure of the baryonic mass content. After applying the quadratic correction term $M_{{\rm baryon}}-M_{{\rm grav}}\approx 0.075\, M_{{\rm grav}}^{2}$ (Timmes et al. 1996), one can obtain

$$\begin{equation} M_{{\rm birth}}\sim 1.08\hbox{--}1.57\, {M}_{\odot } \end{equation} \tag{ 3 }$$

as a viable range for gravitational NS masses, which is believed to encapsulate the range of NS birth masses.

There is considerable evidence that at least some millisecond pulsars have evolved from a first generation of NSs that have accumulated mass and angular momentum from their evolved companion (Alpar et al. 1982; Radhakrishnan & Srinivasan 1982; Wijnands & van der Klis 1998; Markwardt et al. 2002; Galloway et al. 2002, 2005). There is also a line of arguments that supports the possibility of alternative evolutionary processes that may enrich the millisecond pulsar population (Bailyn & Grindlay 1990; Kiziltan & Thorsett 2009).

Possible production channels for isolated millisecond pulsars are mergers of compact primaries or accretion-induced collapse (AIC). In the case where an NS is produced via AIC, the final mass configuration of the remnant is determined by the central density of the progenitor (C–O or O–Ne WD) and the speed at which the conductive deflagration propagates (Woosley & Weaver 1992).

While there are uncertainties on the parameters that describe the ignition and flame propagation, a careful treatment of the physics that tune the transition of an accreting WD yields a unique baryonic mass M_baryon ≈ 1.39 M_☉ for the remnant, which gives a gravitational mass of M_grav ∼ 1.27 M_☉ for NSs produced via AIC (Timmes et al. 1996). There is indirect evidence that the occurrence rate of AICs can be significant (Bailyn & Grindlay 1990).

The physics of these production channels is still not understood well enough to make quantitative predictions of the NS mass distribution produced via these processes. But we can estimate the mass required to spin NSs up to millisecond periods by using timescale and angular momentum arguments.

For low-mass X-ray binaries (LMXBs) accreting at typical rates of $\dot{m}\sim 10^{-3}\,\dot{{M}}_{{\rm Edd}}$, the amount of mass accreted onto an NS in 10¹⁰ yr is Δm ≈ 0.10 M_☉. One can also estimate the amount of angular momentum required to spin the accreting progenitor up to velocities that equal the Keplerian velocity at the co-rotation radius. In order to transfer sufficient angular momentum (L = I × ω) and spin up a normal pulsar (R ≈ 12 km, I ≈ 1.4 × 10⁴⁵ g cm²) to millisecond periods, an additional mass of Δm ≈ 0.20 M_☉ is required. Hence,

$$\begin{eqnarray} &&\Delta m_{{\rm acc}} \approx 0.10\hbox{--}0.20\, {M}_{\odot } \end{eqnarray} \tag{ 4 }$$

will be sufficient to recycle NS primaries into millisecond pulsars.

The mass and the composition of NSs are intricately related. One of the most important empirical clues that would lead to constraints on a wide range of physical processes is the maximum mass of NSs. For instance, secure constraints on the maximum mass provide insight into the range of viable EOSs for matter at supranuclear densities.

A first-order theoretical upper limit can be obtained by numerically integrating the Oppenheimer–Volkoff equations for a low-density EOS at the lowest energy state of the nuclei (Baym et al. 1971). This yields an extreme upper bound to the maximum mass of an NS at M_max ∼ 3.2 M_☉ (Rhoades & Ruffini 1974). For any compact star to stably support masses beyond this limit requires stronger short-range repulsive nuclear forces that stiffen the EOSs beyond the causal limit. For cases in which causality is not a requisite (v → ∞), an upper limit still exists in general relativity (≈5.2 M_☉) that considers uniform density spheres (Shapiro & Teukolsky 1983). However, for these cases, the extremely stiff EOSs that require the sound speed to be superluminal (dP/dρ ⩾ c²) are considered non-physical.

Differentially rotating NSs that can support significantly more mass than uniform rotators can be temporarily produced by binary mergers (Baumgarte et al. 2000). While differential rotation provides excess radial stability against collapse, even for modest magnetic fields, magnetic braking and viscous forces will inevitably bring differentially rotating objects into uniform rotation (Shapiro 2000). Therefore, radio pulsars can be treated as uniform rotators when calculating the maximum NS mass.

While general relativity, along with the causal limit put a strict upper limit on the maximum NS mass at ∼3.2 M_☉, the lower bound is mostly determined by the still-unknown EOS of matter at these densities and therefore is not well constrained. There are modern EOSs with detailed inclusions of nuclear processes such as kaon condensation and nucleon–nucleon scattering that affect the stiffness. These EOSs give a range of 1.5–2.2 M_☉ as the lower bound for the maximum NS mass (Thorsson et al. 1994; Kalogera & Baym 1996). Although these lower bounds for a maximum NS mass are implied for a variety of more realistic EOSs, it is still unclear whether any of these values are favored. Therefore,

$$\begin{eqnarray} &&M_{{\rm max}} \sim 1.5\hbox{--}3.2\, {M}_{\odot } \end{eqnarray} \tag{ 5 }$$

can be considered a secure range for the maximum NS mass value.

Here, we present the statistical model to estimate the NS mass distribution. The approach is based on a formulation that incorporates errors in the measurements of NS mass estimates. Specifically, the model formulation:

$$\begin{equation} m_{i} = \mathcal {M}_{i}+w_{i}, \,\,\,\,\, i=1,\ldots ,n, \end{equation} \tag{ 12 }$$

where, for the ith NS, m_i is the estimate of the NS mass $\mathcal {M}_{i}$ and w_i is the associated error. We thus need a model for the NS mass distribution and the measurement error distribution. Evidently, the key focus of inference is the NS mass distribution, but a flexible specification for the error distribution is needed to ensure that this inference is not biased. At the same time, the model specification must take into account the limited amount of data. The proposed modeling approach achieves a balance between these considerations and, importantly, enables a relatively straightforward implementation of inference through posterior simulation computational methods.

Visual inspection of the pulsar mass estimates (see Tables 1 and 2 and Figure 1) suggests that skewness may be present in the NS mass distribution, at least for the NS–WD systems. It is therefore important to extend the normality assumption, which is implicit in the existing estimation methods. Furthermore, it is clear from the error bars of the pulsar mass estimates that an asymmetric measurement error distribution is needed for some of the observations, especially for the DNS systems. The statistical model developed below allows for skewness both in the NS mass distribution and the error distribution while encompassing the normal distribution for either as a special case.

The pulsars in less constrained systems (e.g., with only one PK parameter determined) typically have asymmetric measurement errors. The flexibility of the statistical modeling approach developed here allows us to take full advantage of all available mass measurements in Tables 1 and 2. It is noteworthy that the model is generic enough so it can be adopted to other similar astrophysical problems and serve as a useful reference.

Regarding the model for the NS mass distribution, we work with a skewed normal distribution with a density function given by

$$\begin{equation} {\rm SN}(\mathcal {M} \mid \mu , \sigma , \alpha) = \frac{2}{\sigma } \phi \left(\frac{\mathcal {M} - \mu }{\sigma } \right) \Phi \left(\frac{(\mathcal {M} - \mu) \alpha }{\sigma } \right), \end{equation} \tag{ 13 }$$

where ϕ(· ) and Φ(· ) denote the density function and cumulative distribution function (CDF), respectively, of the standard normal distribution. Here, $\mu \in \mathbb {R}$ is a location parameter, $\sigma \in \mathbb {R}^{+}$ is a scale parameter, and $\alpha \in \mathbb {R}$ is a skewness parameter. This model was studied by Azzalini (1985) and is one of the more commonly used skewed normal distributions. Note that α = 0 yields the normal distribution with mean μ and standard deviation σ as a special case of Equation (13), which highlights the role of α as a skewness parameter. In particular, positive/negative values of α result in right/left skewness for the density in Equation (13). Hence, an appealing feature of this model is that, within the context of Bayesian inference, we can make probabilistic assessments for skewness of the NS mass distribution relative to a normal distribution through, for instance, a posterior interval estimate for parameter α. As we will discuss in Section 4.2, we find some evidence for skewness in the NS mass distribution corresponding to the NS–WD systems, but not for the DNS systems.

Next, we describe the model for the error distribution, which is motivated by the process used to produce the pulsar mass estimates and the associated error bars. For each pulsar (either from a NS–WD or a DNS system), an empirical density curve for its mass is constructed based on how well the PK parameters of the system can be constrained. We generically denote the final constructed density for the ith pulsar as h_i(m) and note that, although it is unimodal, it may be asymmetric (especially for pulsars that are in a system for which only one PK parameter can be constrained), resulting in the asymmetric error bars reported for some of the systems in Tables 1 and 2. The pulsar mass estimate, m_i, is obtained as the mode of this density, whereas the error bars, +u_i / − ℓ_i, define the interval, (m_i − ℓ_i, m_i + u_i), of the smallest possible length with 68% probability coverage. Numerically, the interval is obtained by starting from the peak value of h_i(m) and slicing down the density until $\int _{m_{i} - \ell _{i}}^{m_{i} + u_{i}} h_{i}(m) \, {d}m = 0.68$, which, given the unimodal shape of density h_i(m), also implies that h_i(m_i − ℓ_i) = h_i(m_i + u_i).

Hence, in the context of model given by Equation (12), the errors w_i are realizations from a distribution with a mode at 0. We thus seek a measurement error distribution that is parameterized in terms of its mode, allows asymmetry around the mode, and yields the normal distribution in the special case of symmetry. A particularly suitable choice is the asymmetric normal distribution studied in Fernández & Steel (1998) with a density function given by

$$\begin{eqnarray} &&{\rm AN}(w \mid c,d) = \frac{2}{d\left(c+\frac{1}{c}\right)} \nonumber\\ &&\quad \times \left\lbrace \phi \left(\frac{w}{cd}\right) 1_{[0,\infty)}(w) + \phi \left(\frac{c w}{d}\right)1_{(-\infty ,0)}(w) \right\rbrace ,\qquad\qquad \end{eqnarray} \tag{ 14 }$$

where c > 0, d > 0, and 1_A(· ) denotes the indicator function of set A. The mode of this density is at 0, when c = 1 it reduces to the normal density with mean 0 and standard deviation d, and when c > 1 (c < 1) it is right skewed (left skewed). Therefore, d plays the role of a scale parameter, whereas c is the asymmetry parameter.

A critically important feature of the asymmetric normal density in Equation (14) is that it leads to a straightforward estimation of parameters c_i and d_i for the ith pulsar, using the values of the 68% central limits, +u_i / − ℓ_i, in Tables 1 and 2. First, from the condition AN(− ℓ_i∣c_i, d_i) = AN(u_i∣c_i, d_i) we obtain $\phi (- c_{i} d_{i}^{-1} \ell _{i})=$ $\phi (c_{i}^{-1} d_{i}^{-1} u_{i})$, which, in turn, yields c_i = (u_i/ℓ_i)^1/2. Next, with c_i determined, we specify d_i by solving numerically for its value that satisfies $\int _{-\ell _{i}}^{u_{i}} {\rm AN}(w \mid c_{i},d_{i}) {d}w = 0.68$. Computing this equation involves normal distribution function evaluations and it can be easily shown that there is a unique solution for d_i. Note that c_i = 1 when ℓ_i = u_i and thus the error distribution is normal for the pulsars with symmetric error bars. However, for pulsars with asymmetric error bars, asymmetry is introduced in the respective error distribution components. The extent of the asymmetry depends on the relative magnitude of ℓ_i and u_i; for instance, the maximum value for the asymmetry parameter c_i is 2.02 corresponding to PSR J1518+4904's companion and the minimum value is 0.38 corresponding to PSR J0514–4002A.

For each of the NS–WD and DNS systems, the data vector comprises data = {(m_i, c_i, d_i): i = 1, ..., n}, with (c_i, d_i) computed as detailed in Section 4.1. Then, combining the models for the NS mass and error distributions in Equations (13) and (14), respectively, the hierarchical model for the data can be written as

$$\begin{eqnarray} &&m_{i} \mid \mathcal {M}_{i} \, \sim \, {\rm AN}(m_{i} - \mathcal {M}_{i} \mid c_{i},d_{i}), \,\,\, i=1,\ldots ,n, \nonumber \\ &&\mathcal {M}_{i} \mid \mu ,\sigma ,\alpha \, \sim \, {\rm SN}(\mathcal {M}_{i} \mid \mu , \sigma , \alpha), \,\,\, i=1,\ldots ,n. \end{eqnarray} \tag{ 15 }$$

That is, given the respective NS masses $\mathcal {M}_{i}$, the pulsar mass estimates m_i arise conditionally independently from the asymmetric normal response distribution with mode at $\mathcal {M}_{i}$. Here, the $\mathcal {M}_{i}$ are modeled as conditionally independent realizations, given parameters (μ, σ, α), from the skewed normal NS mass distribution.

Now, the likelihood function for the NS mass distribution parameters arises from the hierarchical model in Equation (15) by marginalizing over $\mathcal {M}_{i}$, that is,

$$\begin{eqnarray*} \mathcal {L}(\mu ,\sigma ,\alpha ;{\rm data})&=& \prod _{i=1}^{n} \int {\rm AN}(m_{i} - \mathcal {M}_{i} \mid c_{i},d_{i}) \nonumber \\ &&\times\,{\rm SN}(\mathcal {M}_{i} \mid \mu , \sigma , \alpha) \, {d}\mathcal {M}_{i}. \nonumber \end{eqnarray*}$$

The integral is readily available analytically only in the special case of a normal distribution for both the errors and NS masses (c_i = 1 and α = 0, respectively). Therefore, in general, numerical maximization of the likelihood function to obtain the maximum likelihood estimates for parameters μ, σ, and α is not straightforward. Even more challenging is the uncertainty quantification for the point estimates and its subsequent effect on the NS mass density; this would require large-sample asymptotic results, the use of which is particularly problematic to justify given the small number of observations from both the DNS and NS–WD systems.

We thus employ the Bayesian approach to inference for the NS mass distribution using Markov Chain Monte Carlo (MCMC) methods for sampling from the posterior distribution of model parameters (Gelman et al. 2003). Under the Bayesian approach, the model given by Equation (15) is completed with priors for the NS mass distribution parameters. Details on the prior distributions are given in Appendix A.

In addition to providing a coherent probabilistic framework for inference, the Bayesian model formulation enables harnessing the full hierarchical structure in Equation (15). In particular, our posterior simulation method retains the individual NS masses $\mathcal {M}_{i}$ as part of the full parameter vector along with the NS mass distribution parameters (μ, σ, α). In fact, we work with a transformed version (μ, τ², ψ) of (μ, σ, α), such that σ² = τ² + ψ² and α = ψ/τ, with τ > 0 and $\psi \in \mathbb {R}$. This re-parameterization facilitates implementation of the computational method for posterior inference as an efficient Gibbs sampler algorithm. Key in this direction is also a stochastic representation of the skewed normal distribution in Equation (13) as a mixture of normal distributions (Henze 1986). The specific result is given in Appendix A, which includes also the technical details of the MCMC posterior simulation method.

Implementing the Gibbs sampler described in Appendix A, we obtain posterior samples for (μ, τ², ψ) and thus, through the transformation discussed above, for parameters (μ, σ, α). These samples can be used to explore a variety of inferences for the NS mass distribution.

First, the point estimate for the density of the NS mass distribution is given by the posterior predictive density, $\mathcal {P}(\mathcal {M}_{0} \mid {\rm data})$, where $\mathcal {M}_{0}$ denotes the (unknown) mass of a "new" unobserved pulsar that we seek to estimate (predict) given the observed data. The posterior predictive density is given by

$$\begin{eqnarray} \mathcal {P}(\mathcal {M}_{0} \mid {\rm data}) &=& \int {\rm SN}(\mathcal {M}_{0} \mid \mu , \sigma , \alpha) \nonumber\\ &&\times\,p(\mu ,\sigma ,\alpha \mid {\rm data}) \, {d}\mu \, {d}\sigma \, {d}\alpha , \end{eqnarray} \tag{ 16 }$$

where p(μ, σ, α∣data) denotes the posterior distribution of the model parameters. Using expression (16) along with the samples from p(μ, σ, α∣data), we can compute through straightforward Monte Carlo integration the NS mass density estimates over any grid of mass values of interest. Figure 2 plots the posterior predictive NS mass density estimates for the DNS and NS–WD systems, which have peaks at 1.33 M_☉ and 1.55 M_☉, respectively. We can also sample from the posterior predictive distribution by sampling from the ${\rm SN}(\mathcal {M}_{0} \mid \mu , \sigma , \alpha)$ distribution (using its normal mixture stochastic representation) for each posterior sample of (μ, σ, α). The resulting samples quantify the posterior predictive uncertainty around the NS mass density peaks. In particular, for the DNS systems, the 68% and 95% posterior predictive intervals are (1.21 M_☉, 1.43 M_☉) and (1.10 M_☉, 1.55 M_☉), whereas for the NS–WD systems the corresponding intervals are given by (1.35 M_☉, 1.81 M_☉) and (1.13 M_☉, 2.07 M_☉).

Zoom In Zoom Out Reset image size

It is noteworthy from Figure 2 that the NS–WD systems posterior predictive density suggests positive skewness in the NS mass distribution. This is also reflected in the posterior distribution for the skewness parameter α given the NS–WD systems data; specifically, the posterior mean for α is 0.90 and Pr(α > 0∣data) = 0.78. In contrast, the corresponding results for the DNS systems data are E(α∣data) = −0.03 and Pr(α > 0∣data) = 0.49 supporting symmetry (and normality) for the DNS systems mass distribution.

Next, we supplement the point estimates in Figure 2 with uncertainty bands for the NS mass density. To this end, using a grid of mass values in 0.5 M_☉ to 2.5 M_☉, we evaluate the skewed normal NS mass density in Equation (13) at each of the posterior samples for its parameters (μ, σ, α). This produces a sample of densities that can be averaged to obtain the posterior mean NS mass density estimate, given by the solid lines in Figure 3. (Formally, the posterior mean estimate is equivalent to the posterior predictive density and thus the solid lines in Figure 3 agree with the estimates in Figure 2.) However, we can now also depict the posterior uncertainty for the entire NS mass density through percentiles from the posterior sample of densities. In Figure 3, we use the 0.025 and 0.975 percentiles and thus the gray bands depicting the posterior uncertainty correspond to 95% interval estimates for the NS mass density.

Zoom In Zoom Out Reset image size

Finally, Figure 3 also plots the prior point and the 95% interval estimates for the NS mass density. These are produced as discussed above for the posterior inference results, but in this case using samples from the prior distribution assigned to the model parameters. Hence, Figure 3 shows the prior guess at the NS mass density (the prior mean estimate given by the dashed line), as well as the extent of variability in the prior (encapsulated by the dotted lines). This provides an effective means to summarize the extent of prior information for the NS mass density incorporated into the model through the specific priors for the model parameters. Moreover, the comparison with the corresponding estimates given the data illustrates the amount of prior-to-posterior learning, which is evidently significant for both the DNS and NS–WD systems.

The predictive performance of the statistical model developed in Sections 4.1 and 4.2 was evaluated using a well-founded technique for Bayesian model checking (see, for example, Chapter 6 of Gelman et al. 2003). Briefly, for each data point (m_i, c_i, d_i), we obtained the posterior predictive distribution, $\mathcal {P}(m_{i}^{{\rm rep}} \!\mid\! {\rm data})$, for replicated response $m_{i}^{{\rm rep}}$, that is, the pulsar mass estimate that we would observe if the experiment that produced the data was to be replicated. Details on sampling from these posterior predictive distributions are included in Appendix B. An indication of how well the model is performing predictively can be obtained by checking where the observed pulsar mass estimate m_i lies within the corresponding posterior predictive density. A relatively large number of observations falling in the tails of the respective predictive densities is indicative of a poor model fit. Under the proposed model, all pulsar mass estimates from both the DNS and NS–WD systems were effectively captured within their corresponding posterior predictive distributions; see the plots in Appendix B. These results provide a further illustration of the predictive power of the model as it replicates the appropriate type of asymmetry for the responses with asymmetric measurement errors.

The first article that reviewed pulsar mass measurements in order to deduce the range of masses NSs can attain was published by Joss & Rappaport (1976). They used mass measurements from five sources (PSR B1913+16, Her X-1, Cen X-3, SMC X-1, and 3U0900−40), which were predominantly X-ray sources, and found a marginally consistent range of 1.4–1.8 M_☉.

Finn (1994) used a Bayesian statistical approach for the first time to infer limits on the NS mass distribution. By using the mass measurements of four radio pulsars (PSRs B1913+16, B1534+12, B2127+11C, and B2303+46), he concluded that NS masses should fall mainly in the range between 1.3–1.6 M_☉.

A comprehensive paper on pulsar masses was published by Thorsett & Chakrabarty (1999). Their analysis based on 26 sources yielded a remarkably tight mass range at $1.38^{+0.06}_{-0.10}\, {M}_{\odot }$. The width of their NS mass inference was mainly driven by the narrow error bands of the DNS mass measurements.

The recent work by Schwab et al. (2010) analyzed masses of 14 sources with an approach based on comparing the CDF with an idealized Gaussian. It is well understood that the Kolmogorov–Smirnov (K-S) test should be used with caution in cases where deviations occur in the tails (Mason & Schuenemeyer 1983). Additionally, even in data samples where the number of outliers in the tails are considerably larger and the associated measurement errors are taken into account, a K-S approach will still remain limited in quantifying the significance of the outliers. Therefore, while the bimodal feature found for the initial mass function (i.e., M_birth) may be consistent with theoretical expectations for remnant masses produced by electron-capture versus Fe-core collapse SNe (Podsiadlowski et al. 2004), the evidence for a deviation of M_birth from a unimodal distribution is still tentative. In order to firmly establish a potential multi-modal feature for the NS birth mass distribution, a more diverse sample tested with statistics that allow sensitivity and performance checks are required.

Our model is flexible and sensitive enough for detecting signatures of a potential truncation in the underlying mass distribution. We pay particular attention to whether there are signatures of truncation in the implied mass distribution at the high-mass end for two reasons: (1) a high-mass limit set by the EOS of the NS matter, or by general relativity, should produce a relatively sharp cut off, which will manifest itself as a truncation. (2) The nature of potential skewness of the underlying mass distribution will provide insight into how NSs are produced and the evolutionary link between the two NS populations.

The masses of NSs in DNS systems imply a tight symmetric distribution while the distribution of masses of NSs in NS–WD systems show evidence of skewness with a heavy tail on the high-mass end. On the other hand, both populations show no evidence of a strong truncation limit on either end.

These have important implications: the stochastic nature of evolutionary processes such as long-term stable accretion naturally produces a wider distribution for NSs in NS–WD systems. This, along with the lack of a strong truncation, indicates that, in particular, the high-mass end of the NS mass distribution is driven by evolutionary constraints. As a result, this rules out the possibility that an upper mass limit is set by the EOS of matter or general relativity for NSs with masses M ≲ 2.1 M_☉. Therefore, the 2.1 M_☉ upper mass limit implied by NSs in NS–WD systems should be considered a minimum secure limit to the maximum NS mass rather than an absolute upper limit to NS masses. The heavy tail of the mass distribution of NSs in NS–WD systems favors the possibility that at least some of these pulsars are born as massive NSs.

All EOSs that require a maximum NS mass M_max ⩽ 2.1 M_☉ are ruled out. The implied stiffness of the EOS largely precludes the presence of meson condensates and hyperons at supranuclear densities. Consequently, lower central densities, larger radii, and thicker crusts for NSs are favored (Shapiro & Teukolsky 1983).

The energy density–radius relation implied by Tolman (1939), when combined with the causality limit, gives an analytical solution for an upper limit on the central density

$$\begin{eqnarray} &&\rho _{c}M^{2}=15.3 \times 10^{15} \, {M}_{\odot } ^{2} \,{\rm g}\,{\rm cm}^{-3}. \end{eqnarray} \tag{ 17 }$$

With a 2.1 M_☉ secure lower limit on the maximum NS mass, we set a 95% confidence upper limit to the central density of NSs, which is

$$\begin{eqnarray} &&\rho _{{\rm max}}<3.47\times 10^{15}\ {\rm g}\,{\rm cm}^{-3} \end{eqnarray} \tag{ 18 }$$

corresponding to ≈11ρ_s for a fiducial saturation threshold n_s ∼ 0.16 fm⁻³.

Exotic matter such as hyperons and Bose condensates significantly reduces the maximum mass of NSs. Therefore, a strict lower limit on the maximum NS mass M_max > 2.1 M_☉ rules out soft EOSs with extreme low-density softening and require the existence of exotic hadronic matter (see Lattimer & Prakash 2007 for a review). NSs with deconfined strange quark matter mostly have maximum predicted masses lower than 2.1 M_☉. Hence, EOSs with strange quark matter that predict maximum masses smaller than 2.1 M_☉ can also be ruled out as viable configurations for NS matter.

A 2.1 M_☉ upper limit to masses of NSs in NS–WD system poses a problem. If all millisecond pulsars were indeed NSs that are recycled from a first generation of normal pulsars, the implied distribution should be consistent with a recycled version of the initial mass distribution. While the peaks of the distributions for DNS and NS–WD systems appear to be consistent with the expectations of standard recycling (Section 2.2), the widths imply otherwise. As shown in Figure 1, Δm_acc = 0.22 M_☉ lies within the expected range. However, with typical accretion rates experienced during the LMXB phase ($\dot{m}_{{\rm acc}} \sim 10^{-3}\,\dot{{M}}_{{\rm Edd}}$), it is difficult to accumulate sufficient mass onto NSs that started their lives as ∼1.4 M_☉ NSs and produce NSs with masses ∼2.1 M_☉ such as PSR B1516+02B. Even with initial masses of ∼1.6 M_☉ these sources need to accrete Δm ≈ 0.4–0.5 M_☉ during their active accretion phase. This requires long-term, stable, active accretion at unusually high rates.

Based on the P–$\dot{{\rm P}}$ demographics of millisecond pulsars, Kiziltan & Thorsett (2009) argue that ≈ 30% of the millisecond pulsar population may be produced via a non-standard evolutionary channel. This prediction falls in line with a distribution that has a consistent recycled peak but has an unusual width (and some skewness with a heavy tail on the high-mass end), which extends up to 2.1 M_☉. While it is difficult to quantify the formation rate(s) of non-standard processes that may produce these NSs, it is clear that the standard scenario requires at least a revision. Such a revision should consistently account for the observed P–$\dot{{\rm P}}$ distribution of millisecond pulsars, along with the long term sustainability of unusually high accretion rates that would be required to produce a population of second generation massive NSs.

The only viable alternative to a major revision of the mass evolution implied by the standard recycling scenario, also corroborated by the lack of truncation of the underlying NS mass distribution and enforced by the heavier tail of the skewed distribution inferred from NS–WD systems, is then to form massive NSs.

The skewed normal distribution in Equation (13) admits a stochastic representation as a mixture of normal distributions with a truncated normal mixing distribution. Introducing truncated normal latent random variables z_i and reparameterizing SN(· ∣μ, σ, α) to SN(· ∣μ, τ², ψ), as described in Section 4.2, we have

$$\begin{eqnarray} \mathcal {M}_{i}\mid z_i,\mu ,\tau ^2,\psi &\stackrel{{\rm ind}}{\sim }&{\rm N}(\mathcal {M}_{i} \mid \mu +\psi z_i,\tau ^2); \nonumber\\ z_i &\stackrel{{\rm ind}}{\sim } &{\rm N}(0,1)1_{[0,\infty)}(z_i) \end{eqnarray} \tag{ A1 }$$

as a hierarchical mixture representation of ${\rm SN}(\mathcal {M}_{i} \mid \mu , \tau ^{2}, \psi)$ (that is, if we marginalize over the z_i in (A1), we obtain $\mathcal {M}_{i} \mid \mu ,\tau ^2,\psi \stackrel{{\rm ind}}{\sim }$$\,{\rm SN}(\mathcal {M}_{i} \mid \mu , \tau ^{2}, \psi)$). Here, N(0, 1)1_{[0, ∞)}(z_i) indicates a standard normal distribution restricted to $\mathbb {R}^+$.

To facilitate MCMC posterior simulation, we work with this formulation for the skewed normal distribution. The second line of the hierarchical model for the data given in Equation (15) is therefore changed as indicated in Equation (A1) to give the following hierarchical model for the data:

$$\begin{equation*} m_i \mid \mathcal {M}_i\stackrel{{\rm ind}}{\sim } {\rm AN}(m_i - \mathcal {M}_i \mid c_i,d_i),\quad i=1,\dots ,n \end{equation*}$$

$$\begin{equation*} \mathcal {M}_{i}\mid z_i,\mu ,\tau ^2,\psi \stackrel{{\rm ind}}{\sim }{\rm N}(\mathcal {M}_{i} \mid \mu +\psi z_i,\tau ^2),\quad i=1,\dots ,n \end{equation*}$$

$$\begin{equation*} z_i\stackrel{{\rm ind}}{\sim }{\rm N}(z_i\mid 0,1)1_{[0,\infty)}(z_i),\quad i=1,\dots ,n. \end{equation*}$$

The Bayesian model is completed with (conditionally conjugate) priors for the NS mass distribution parameters. Specifically, we place a bivariate normal prior, N((μ, ψ)^T∣θ, Σ), on (μ, ψ) with mean vector θ and covariance matrix Σ and an inverse-gamma prior, IG(τ²|a, b), on τ² with shape parameter a > 1 and scale parameter b, such that the prior mean of τ² is b/(a − 1).

In addition to the NS mass distribution parameters (μ, τ², ψ), the model parameters include the individual NS masses $\mathcal {M}_1,\dots ,\mathcal {M}_n$, as well as the auxiliary random variables z₁, ..., z_n. The likelihood function and the priors for (μ, ψ) and τ² combine to give the posterior distribution for all model parameters, which is proportional to:

$$\begin{eqnarray} && {\rm N}((\mu ,\psi)^T\mid \theta ,\Sigma) {\rm IG}(\tau ^2|a,b) \prod _{i=1}^{n} \lbrace {\rm AN}(m_i - \mathcal {M}_i \mid c_i,d_i)\nonumber\\ &&\quad \times\,{\rm N}(\mathcal {M}_i\mid \mu +\psi z_i,\tau ^2){\rm N}(z_i\mid 0,1)1_{[0,\infty)}(z_i)\rbrace. \end{eqnarray} \tag{ A2 }$$

We utilize an MCMC posterior simulation algorithm to sample from the posterior distribution. The MCMC algorithm dynamically updates the model parameters by simulating from their full conditional distributions in turn and, upon reaching convergence, the resulting samples are samples from the posterior distribution that may be used for inference. The key benefit of the augmented hierarchical model that includes the auxiliary variables z₁, ..., z_n is that it allows for the use of a Gibbs sampler, which is both efficient and straightforward to implement, since it involves sampling from standard distributions as described below.

For each i = 1, ..., n, the full posterior conditional distribution for $\mathcal {M}_i$ is proportional to ${\rm AN}(m_i - \mathcal {M}_i \mid c_i,d_i) {\rm N}(\mathcal {M}_i\mid \mu +\psi z_i,\tau ^2)$, which can be shown to result in a mixture of two truncated normals. The first truncated normal has mean $\lbrace (\mu +\psi z_i)c_i^2d_i^2+m_i\tau ^2\rbrace /(\tau ^2+c_i^2d_i^2)$ and variance $(\tau ^2c_i^2d_i^2)/(\tau ^2+c_i^2d_i^2)$ and is restricted to the interval (− ∞, m_i]. The second truncated normal has mean $\lbrace (\mu +\psi z_i)(d_i^2/c_i^2) + m_i\tau ^2\rbrace /(\tau ^2+(d_i^2/c_i^2))$ and variance $(\tau ^2d_i^2/c_i^2)/(\tau ^2+(d_i^2/c_i^2))$ and is restricted to the interval (m_i, ∞). The (unnormalized) weight of the first truncated normal is

$$\begin{eqnarray*} &&\exp \left\lbrace -\frac{(\mu +\psi z_i-m_i)^2}{2\left(\tau ^2+c_i^2d_i^2\right)}\right\rbrace \cdot \frac{\tau d_ic_i}{\left(\tau ^2+d_i^2c_i^2\right)^{1/2}}\nonumber\\ &&\quad \times\left\lbrace \Phi \left(\frac{m_i-\left((\mu +\psi z_i)c_i^2d_i^2+m_i\tau ^2\right)\left/\right. \left(\tau ^2+d_i^2c_i^2\right)}{\tau d_ic_i\left/\right. \left(\tau ^2+d_i^2c_i^2\right)^{1/2}}\right)\right\rbrace \end{eqnarray*}$$

and that of the second is

$$\begin{eqnarray*} &&{{ \displaystyle \exp \left\lbrace -\displaystyle\frac{(\mu +\psi z_i-m_i)^2}{2\left(\tau ^2+\left(d_i^2\left/\right. c_i^2\right)\right)}\right\rbrace \cdot \displaystyle\frac{\tau d_i/c_i}{\left(\tau ^2+\left(d_i^2\left/\right. c_i^2\right)\right)^{1/2}}}} \nonumber\\ && {{ \displaystyle\quad \times \left\lbrace 1\,{-}\,\Phi \displaystyle\left(\frac{m_i\,{-}\,\left((\mu +\psi z_i)\left(d_i^2\left/\right. c_i^2\right)+m_i\tau ^2\right)\left/\right. \left(\tau ^2+\left(d_i^2\left/\right. c_i^2\right)\right)}{(\tau d_i/c_i)\left/\right. \left(\tau ^2 + \left(d_i^2\left/\right. c_i^2\right)\right)^{1/2}}\right)\right\rbrace}}. \end{eqnarray*}$$

Each $\mathcal {M}_i$ is therefore updated at each MCMC iteration by drawing from a mixture of two truncated normals, with the parameters given above.

Next, the full conditional distribution for z_i is proportional to ${\rm N}(\mathcal {M}_{i} \mid \mu +\psi z_i,\tau ^2)N(z_i\mid 0,1)1_{[0,\infty)}(z_i)$, which results in a truncated normal distribution with mean $(\mathcal {M}_i-\mu)\psi /(\psi ^2+\tau ^2)$ and variance τ²/(ψ² + τ²), restricted to [0, ∞).

The NS mass distribution parameters can be sampled with standard updates. The posterior full conditional distribution for (μ, ψ) can be derived as a bivariate normal distribution,

$$\begin{eqnarray} &&{\rm N}((\mu ,\psi)^T\mid (\Sigma ^{-1}+Z^TZ/\tau ^2)^{-1} (\Sigma ^{-1}\theta +Z^T\boldsymbol{\mathcal {M}}/\tau ^2),\nonumber\\ &&\quad \times(\Sigma ^{-1}+Z^TZ/\tau ^2)^{-1}), \end{eqnarray} \tag{ A3 }$$

where $\boldsymbol{\mathcal {M}}=(\mathcal {M}_1,\dots ,\mathcal {M}_n)^T$ and Z is an n × 2 matrix with Z_i1 = 1 and Z_i2 = z_i, for i = 1, ..., n. Lastly, the full conditional for τ² is given by an inverse-gamma distribution with updated shape and scale parameters, ${\rm IG}(a+n/2,b+\sum _{i=1}^{n}(\mathcal {M}_i-\mu -\psi z_i)^2/2)$.

The Gibbs sampler is therefore implemented with the following general procedure:

• 1.  
  Begin with initial values $(\lbrace \mathcal {M}_{1}^{(0)}\!,\ldots ,\mathcal {M}_{n}^{(0)} \rbrace , \lbrace z_{1}^{(0)},\ldots ,z_{n}^{(0)} \rbrace {,} \mu ^{(0)},\psi ^{(0)},\tau ^{2(0)})$.
• 2.  
  If the current sample at iteration t is $(\lbrace \mathcal {M}_{1}^{(t)},\ldots ,\mathcal {M}_{n}^{(t)} \rbrace , \lbrace z_{1}^{(t)},\ldots ,z_{n}^{(t)} \rbrace ,\mu ^{(t)},\psi ^{(t)},\tau ^{2(t)})$, obtain the next sample by simulating from the following distributions:
• 3.  
  Repeat the previous step, for t = 1, ..., T.

The chain appeared to converge very rapidly; when two chains were run with different initial values, they converged within just a few iterations. There was some autocorrelation present for μ and ψ, but it fell below 0.2 after approximately a lag of 5 and samples were approximately uncorrelated by a lag of 10. In light of these results, all inferences are based on 10,000 posterior samples, which came from a longer chain that was thinned, retaining every 10 posterior samples.

Hyperparameters θ = (θ_μ, θ_ψ)^T, Σ, a, and b must be specified in this Bayesian model. First, note that $E(\mathcal {M})=$θ_μ + θ_ψE(z) and that the skewness parameter of the NS mass distribution is α = ψ/τ. To avoid favoring skewness in the prior distribution, we set E(ψ) = θ_ψ = 0 and then set θ_μ to a reasonable prior location for the NS mass distribution for each system. The set of priors used for the results in Section 4.2 had θ_μ = 1.3 for the DNS system and θ_μ = 1.45 for the NS–WD system. With θ_ψ = 0, we have that $Var(\mathcal {M})=$ b/(a − 1) + Σ_μ + Σ_ψ, where (Σ_μ, Σ_ψ) are the elements of matrix Σ, which is taken to be diagonal. We fix a = 3 in each set of priors and for the DNS system we set b/(a − 1), Σ_μ, and Σ_ψ all to 0.25²/3. For the NS–WD system, we set all three expressions to 0.3²/3. These result in a fairly dispersed (thus, non-informative) prior distribution for NS masses, as demonstrated by the prior interval estimates in Figure 3.