Maximum mass cutoff in the neutron star mass distribution and the prospect of forming supramassive objects in the double neutron star mergers
Dong-Sheng Shao, Shao-Peng Tang, Jin-Liang Jiang, Yi-Zhong Fan
MMaximum mass cuto ff in the neutron star mass distribution and the prospect of formingsupramassive objects in the double neutron star mergers Dong-Sheng Shao,
1, 2
Shao-Peng Tang,
1, 2
Jin-Liang Jiang,
1, 2 and Yi-Zhong Fan
1, 2, ∗ Key Laboratory of Dark Matter and Space Astronomy,Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210033, China School of Astronomy and Space Science, University of Science and Technology of China, Hefei, Anhui 230026, China
The sample of neutron stars with a measured mass is growing quickly. With the latest sample, we adopt botha flexible Gaussian mixture model and a Gaussian plus Cauchy-Lorentz component model to infer the massdistribution of neutron stars and use the Bayesian model selection to explore evidence for multimodality anda sharp cuto ff in the mass distribution. The two models yield rather similar results. Consistent with previousstudies, we find evidence for a bimodal distribution together with a cuto ff at a mass of M max = . + . − . M (cid:12) (68%credible interval; for the Gaussian mixture model). If such a cuto ff is interpreted as the maximum gravitationalmass of nonrotating cold neutron stars, the prospect of forming supramassive remnants is found to be quitepromising for the double neutron star mergers with a total gravitational mass less than or equal to 2.7 M (cid:12) unlessthe thermal pions could substantially soften the equation of state for the very hot neutron star matter. Thesesupramassive remnants have a typical kinetic rotational energy of approximately 1 − × ergs. Togetherwith a high neutron star merger rate approximately 10 Gpc − yr − , the neutron star mergers are expected to besignificant sources of EeV(10 eV) cosmic-ray protons. PACS numbers:Keywords: Neutron stars − Compact binary stars − Gravitational waves
I. INTRODUCTION
The mass distribution of neutron stars (NSs) is very helpful in revealing the mechanism of supernova explosion, the accretiondynamics of binary neutron star, the equation of state of matter with ultrahigh density, the mechanism of cosmic-ray acceleration,etc. Since the discovery by Hewish et al. [1], more than 2500 NSs have been detected. The mass measurements, however, aremuch more challenging and such a goal has just been achieved for a small fraction of these extreme objects, usually in NS binarysystems. Nevertheless, benefiting from improvements on pulsar radio timing and x-ray observation, the growing population ofthe NSs with reliable mass measurements within decades makes it feasible to statistically infer the features of the distribution[2–6]. A bimodal distribution with two peaks at approximately 1 . M (cid:12) and approximately 1 . − . M (cid:12) was suggested in theabove literature, which can be well explained by di ff erent formation channels and evolution scenarios [7]. The other recentintriguing / remarkable finding is a significant cuto ff at the high end of the mass distribution ( M max = . + . − . M (cid:12) ), which is mostlikely the maximum gravitational mass ( M TOV ) of the nonrotating neutron star and plays an important role in bounding theequation of state (EoS) of NS matter [11]. It also directly sets a robust lower limit on the mass of stellar-origin black holes.Since the publication of Alsing et al. [11], rapid progress has been made on the mass measurements of NSs. In particular,several NSs are found to be very massive. For instance, PSR J1600 − . + . − . M (cid:12) [12], PSR J0740 + . + . − . M (cid:12) [13], PSR J1959 + . ± . M (cid:12) , PSR J2215 + . + . − . M (cid:12) [14](please note that the other group reported a mass of 2 . + . − . M (cid:12) [15]), and J1811 − . + . − . M (cid:12) [16]. In thiswork the uncertainties are for 68.3% confidence level unless specifically noticed. Therefore, it is necessary to update the analysisof Alsing et al. [11] with the latest sample of NSs with mass measurements / information and new fit functions. That is the mainmotivation of this study.Our work is structured as follows. In Sec. II we collect the mass measurements of NSs from the latest literature and thenanalyze the maximum mass cuto ff in the mass distribution. In Sec. III, as a direct application of the inferred maximum masscuto ff , we adopt an EoS-insensitive approach to examine the prospect of forming supramassive neutron stars (SMNSs) in thedouble neutron stars(DNS) mergers and discuss the possibility that the EoS softening e ff ect of thermal pions generated in thevery hot neutron star can be probed. Motivated by a promising formation prospect, we estimate the kinetic rotational energy ofthese SMNSs (again, in an EoS-insensitive way) and then their role in accelerating EeV cosmic rays. Finally we summarize ourresults with some discussions. ∗ Corresponding author. [email protected] Though some objects are millisecond pulsars, their rotations are not quick enough to enhance the maximum gravitational mass e ff ectively (see, e.g., Refs.[8–10]). a r X i v : . [ a s t r o - ph . H E ] S e p II. UPDATED ESTIMATE OF THE MAXIMUM MASS CUTOFF IN THE NEUTRON STAR MASS DISTRIBUTIONA. Neutron star mass measurements
Up to now, almost all the reliable mass measurements were carried out for neutron stars in binary systems. Benefiting fromKepler’s Third Law, the orbital parameters of those neutron stars and their companions, which can be measured by either radiotiming of pulsations or x-ray / optical observations, make it possible to determine the masses of NSs (see Refs. [17, 18] for recentreviews).Generally in the Newtonian frame, we can measure five Keplerian parameters for the orbital motions of binary: the binaryperiod P b , the eccentricity e , the component of the pulsar’s semimajor axis a p along the line of sight x p = a p sin i / c (where i isthe orbital inclination angle and c is the speed of light), and the time and longitude of periastron T and ω . Then, the so-calledmass function of the binary, defined as f ≡ ( m c sin i ) ( m p + m c ) = (cid:16) π P b (cid:17) x G [17], is dependent on P b and i , where m p and m c stand forthe masses of pulsar and its companion respectively, and G is the Newtonian gravitational constant. The degeneracies of theunknowns can be broken as long as the mass ratio q ≡ m p / m c = x c / x p (where x c is the projection of the companion’s semimajoraxis on the line of sight) and the mass of the companion are determined, as briefly summarized below (see Ref.[18] and thereferences therein). For some binary systems with a main-sequence star or bright white dwarf companion, q and m c can bemeasured by studying the spectrum of the companion (e.g., the Balmer line of hydrogen in the atmosphere via phase-resolvedoptical spectroscopy), leading to a reliable mass measurement. For double NS systems (or binaries with massive white dwarfcompanion), the components are compact enough to make the relativistic e ff ects on the orbital motion observation. These e ff ects,described by five post-Keplerian (PK) parameters, including the periastron precession ˙ ω , Einstein delay γ , the shape and rangeof the Shapiro delay s and r , and the orbital period decay ˙ P b , depend sensitively on the masses of components and the Keplerianparameters of their orbit (see Ref.[19] and references therein). Once some of them have been precisely measured by radio timingof pulsars, the individual mass can be determined with the least uncertainty, especially when both of the components happenedto be pulsars. For neutron stars with high or low stellar mass companion, the observations from x-ray / optical bands provide aviable approach to determine the orbital parameters and the masses. The measurement of eclipsing of the x ray from the high-mass x-ray binary (HMXB) yields some fundamental parameters of binary orbit, such as the period P b , the eccentricity e , thelongitude of periastron ω , and the semimajor axis of the neutron star’s orbit a x sin i . Together with the information of velocityand inclination of companion obtained from optical observations, the mass can be solved from the basic equations. Modeling thethermal emission of neutron star atmosphere can constrain both mass and radius of quiescent low-mass x-ray binary (qLMXB).The thermonuclear x-ray burst, a helium flash occurring in the surface layer of the accreting neutron star of the low-mass x-raybinary(LMXB), can also be used to measure mass and radius, by combining analysis of the apparent angular size, the Eddingtonflux, and the source distance. Most events included in our sample (TableI) are measured in the above ways.For the merging NSs in the deep universe, their masses can be measured with the gravitational wave data. Until January 2020,only two neutron star merger events had been reported. The GW170817 data strongly favored the double NS merger origin [20],while for GW190425, a NS-black hole binary system [21, 22] cannot be ruled out. In this work following the LIGO and Virgocollaborations [23], we attribute it to the merger of a pair of relatively massive NSs. The double neutron star merger sample,though still small right now, is expected to increase rapidly in the next decade. For completeness we include these objects inour analysis. As for the isolated neutron stars, the pulse-profile modeling can simultaneously yield the masses and radii of somenearby bright millisecond pulsars. With the NICER data, such a goal was achieved first for PSR J0030 + et al. [26] proposed a new method to infer the masses of a few isolatedneutron stars with the gravitational redshift measurements. Though interesting, such an approach is model dependent, and wedo not include these events in our sample.We updated the sample listed in Table 1 of Alsing et al. [11] in two aspects. First, the events with improved mass mea-surements have been updated. Second, the new events with mass measurements, dated from April 2018, including a few verymassive ones such as PSR J0740 + + + − et al. [11], thetotal number of NSs in our sample increased from 74 to 103 (see TableI in the Appendix for details). B. Evidence for a maximum mass cuto ff in the neutron star mass distribution: Updated analysis Benefiting from the accumulated NS mass measurements mentioned above, we can statistically investigate the propertiesof NS mass distribution. Previous works have proposed various models of NS population to fit the observation data, such asuniform, single Gaussian, bimodal Gaussian [27], multicomponent Gaussian [11], and skewed normal distribution [5]. In thiswork, we use two models to fit the latest sample. The first is a two-component Gaussian mixture model with a maximum masscuto ff M max , i.e., P ( m p | (cid:126)θ ) = r N ( m p | µ , σ ) / Φ + (1 − r ) N ( m p | µ , σ ) / Φ , for m p ∈ [ M min , M max ] , (1) FIG. 1: Distributions of parameters ( (cid:126)θ ) of the NS mass distribution obtained with some di ff erent priors and models. where (cid:126)θ = { µ , µ , σ , σ , r , M max } , and inside the brackets, µ ( µ ), σ ( σ ), r , and m p , respectively, denote the mean, standarddeviation, relative weight of the two components, and the pulsar mass. The second model is a mixture of Gaussian component N and a Cauchy-Lorentz component Ca , whose function form is given by replacing the second Gaussian component N ofEq.(1) with Ca ; thus, µ and σ represent the location and scale parameters of Cauchy-Lorentz distribution. In our analysis, weset M min = . M (cid:12) , which is a reasonable lower bound to contain most of the NS mass measurements to probe the distributionproperties. The normalization constants Φ i ( i = ,
2) are integrals over each component ( · ) i (over the allowed NS mass range)using Φ i = (cid:82) M max M min ( · ) i ( x | µ i , σ i )d x . To approximate well the non-Gaussian mass measurements in TableI, we use the asymmetricnormal distribution studied in Refs.[5, 28] to reproduce the error distribution, of which the density function is given byAN( w | c , d ) = d ( c + c ) (cid:26) φ (cid:18) wcd (cid:19) [0 , ∞ ) ( w ) + φ (cid:18) cwd (cid:19) ( −∞ , ( w ) (cid:27) , (2)where c > d > φ means normal distribution, and 1 A ( · ) denotes the indicator function of set A . Thus given the NS massmeasurements M i + u i − (cid:96) i ( + u i / − (cid:96) i are 68% central limits), parameters c i and d i for the i th NS can be estimated through c i = ( u i /(cid:96) i ) / and (cid:82) u i − (cid:96) i AN( w | c i , d i )d w = .
68. Then it is straightforward to calculate the probability for a specific pulsar mass m p via P ( D i | m p ) = AN( m p − M i | c i , d i ). For the data only having mass function f and mass ratio q (or total mass m T ) measurementsavailable in TableI, we adopt the Eqs.(3) and (4) of Alsing et al. [11] to evaluate the probability P ( D i | m p ). Therefore, thelikelihood constructed for our inference is given by L ( D | (cid:126)θ ) ∝ N (cid:89) i = (cid:90) P ( m p | (cid:126)θ ) P ( D i | m p )d m p , (3)with which we can use the nest sampling technique, e.g., P y M ulti N est sampler [29], to obtain the samples of the parameters( (cid:126)θ ) of NS mass distribution. The ranges of (cid:126)θ are chosen as follows: µ i ∈ [0 . , . M (cid:12) , σ i ∈ [0 . , M (cid:12) , r ∈ [0 . , . M max ∈ [1 . , . M (cid:12) . We take both uniform and uniform-in-log priors to perform nest sampling, and further request that µ < µ < M max and σ < σ . M / M P D F Best Fit 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 M / M P D F Best Fit
FIG. 2: Maximum aposteriori
NS mass distribution (blue) with 1000 independent posterior samples to give a visual guide for the uncertainties,under the considered model that is most preferred by the data. Left panel: the two-component Gaussian mixture with a sharp cuto ff of approx-imately 2 . M (cid:12) (median value). Right panel: the mixture of Gaussian and Cauchy-Lorentz components with a sharp cuto ff of approximately2 . M (cid:12) (median value). In Bayes’s statistic frame, the method to evaluate the model preference is the odds ratio of the probability of two hypotheseswhich is given by O = P (D | H ) P (D | H ) P ( H ) P ( H ) , (4)where P (D | H ) is the Bayesian evidence (or marginal likelihood) Z for a given model hypothesis H , P (D | H ) / P (D | H ) isthe Bayes factor K , and P ( H ) / P ( H ) is the prior odds ratio that defines our prior relative belief in these two models (here weset to unity as an apriori agnostic). To estimate the “evidence” of the maximum mass cuto ff , we follow the procedure of Alsing et al. [11] by comparing the evidence to the model of fixing M max = . M (cid:12) (i.e., without mass cuto ff ). For the two-componentGaussian model, the logarithm evidence of with (without) cuto ff is − .
47 ( − . − .
04 ( − . ff [30], consistent with the previouswork [11].Some results are reported in Fig.1, and for two-component Gaussian model, we have M max = . + . − . M (cid:12) (68% credibleinterval; the 95% credible interval is M max = . + . − . M (cid:12) ). We have also tested the removal of one or more very massiveneutron stars whose mass measurement methods are not so direct / widely accepted, e.g., PSR J2215 + + + M max = . + . − . M (cid:12) . If we further remove PSR J1959 + M max = . + . − . M (cid:12) , while for the mixture of Gaussian and Cauchy-Lorentz model, we have M max = . + . − . M (cid:12) , 2 . + . − . M (cid:12) ,and 2 . + . − . M (cid:12) for the full data and the sequential removal of PSR J2215 + + M max ; the inclusionof “heavier” NSs can e ff ectively shift the bounds on M max . In view of these facts and considering that the mass measurementsof Ref.[14] may su ff er from some systematic uncertainties, for a cross-check we instead adopt the mass of 2 . + . − . M (cid:12) for PSRJ2215 + et al. [15] in our modeling and infer the M max to be 2 . + . − . M (cid:12) (two-Gaussian-componentmodel) and 2 . + . − . M (cid:12) (Gaussian plus Cauchy-Lorentz component model), respectively. Since they are very consistent withthe result using the data of Ref.[14], we still take our results based on the “latest” mass measurement sample as the fiducialones. We show in Fig.2 the maximum aposteriori mass distribution for the two models with free M max with 1000 independentposterior samples plotted over the top to give a visual impression of the uncertainties on the shape of the distribution. We takethe best-fit result of the two-component Gaussian mixture (shown in the blue line) as our fiducial model.Note that our sample, though significantly extended in comparison to Alsing et al. [11], still su ff ers from some selectione ff ects. In particular, almost all the mass data available come from NSs in binaries, and all the most precise mass measurementscome from double NS systems; there is still no strong evidence yet that these events represent well the whole population ofneutron stars. Anyhow, a very recent study shows that the neutron stars, even born in the death of the very massive stars with azero metallicity, have a gravitational mass below 1 . M (cid:12) [31]. Therefore, it is unlikely to find some isolated NSs as massive as M max ∼ . M (cid:12) . (Indeed, the current measurement / estimates of the isolated NSs find a low mass of approximately 1 . − . M (cid:12) [24, 26], though the sample is still quite small.) Moreover, we have shown that the heaviest objects rather than the preciselymeasured neutron star masses in the double neutron star binaries play the key role in governing M max . The more relevantselection e ff ects may arise from the observations / identifications of the neutron star binary systems. To examine this possibilitywe have collected some measured / derived properties of our pulsar sample from the ATNF (Australia Telescope National Facility)catalog [32]. These properties include the spin period, magnetic field strength on the surface, age, and distance available for 46objects with mass measurements. For the radio luminosity and flux, the sample is a bit smaller. Nevertheless, these most massiveevents (except Vela X1 which was detected in optical and x-rays) have been included. Below, we focus on the objects that arepossibly more massive than 2 M (cid:12) (see the colored points in Fig.3). All these very massive objects are rather old (greater than orequal to 10 yrs) and have very low surface magnetic fields (approximately a few 10 G; the “exception” is PSR J0348 + . × G). Their rotation periods are short (less than 10 ms, exceptPSR J0348 + . + + + / indication for a sizeable nondetecion / misidentification probability of thevery massive NSs in radio. Considering the above facts, in agreement with Alsing et al. [11], we suggest that some selectione ff ects might leave some imprint on the inferred mass distribution, while it seems unlikely that they are responsible for theinferred hard cuto ff at M max . In order to further check whether the radio emission plays an important role in shaping the massdistribution, we have further carried out a Chi-square test of the independence between the masses and radio luminosities ofthe current sample and got a p-value as low as 6 . × − , suggesting no linear correlation among the masses and luminosities.To check potential nonlinear correlation between the two variables, we have also fitted the data with a multivariate adaptiveregression code called pyearth [33], which is e ff ective at identifying a linear or nonlinear relation. The best fit is a constantfunction, indicating the absence of any linear or nonlinear correlation between the mass and luminosity. This is also evident inFig.4 which displays the masses and radio luminosities of 40 pulsars. In the right part we show the distribution of luminosities oftwo groups of NSs separated by M = . M (cid:12) , namely the high-mass group and the low-mass group. The two sided Kolmogorov-Smirnov test has been adopted to examine whether the two groups of luminosities share the same intrinsic distribution or not andwe have got a p-values of 0 .
96, which favors the same intrinsic luminosity distribution hypothesis (a highly relevant / consistentconclusion was also drawn in Ref.[34]). We there ore conclude that there is no evidence for a correlation between the mass andradio luminosity of the neutron stars. As a result of the lack of identification of the mass-dependent selection e ff ects (see, e.g.,Refs. [5, 35] for previous investigations), we do not expect that these selection e ff ects will introduce serious bias to the observedmass distribution (see also Ref.[5]). III. SUPRAMASSIVE NEUTRON STARS: FORMATION PROSPECT IN DOUBLE NS MERGERS AND THEACCELERATION OF EEV COSMIC RAYS
The inferred M max can be used to tighten the bounds on the equation of state of the extreme dense matter of the neutron stars([e.g., Ref. 11]). The other direct / important application is to estimate the fate of the remnants formed in the double neutron starmergers. As shown below, we find a promising prospect of forming supramassive neutron stars (SMNSs) which has a typicalkinetic rotational energy approximately 1 − × erg. With a reasonably high surface magnetic field (i.e., greater than orequal to 10 G), the release of such a huge amount of energy into the surrounding material is quick and the driven energeticblast wave can accelerate the EeV cosmic rays e ff ectively. http: // / research / pulsar / psrcat / , version 1.63. https: // contrib.scikit-learn.org / py-earth / content.html. Mass ( M ) P e r i o d ( s ) J1600-3053J2215+5135J1959+2048J0740+6620J0348+0432J1811-2405
Mass ( M ) B s u r f ( G a u ss ) J1600-3053J2215+5135J1959+2048J0740+6620J0348+0432J1811-2405
Mass ( M ) A g e ( y e a r ) J1600-3053J2215+5135J1959+2048J0740+6620J0348+0432J1811-2405
Mass ( M ) D i s t a n c e ( k p c ) J1600-3053J2215+5135J1959+2048J0740+6620J0348+0432J1811-2405
FIG. 3: Statistical characteristics of pulsar mass vs spin period, magnetic field strength on surface, age and distance. The most massive neutronstars (which may be heavier than 2.0 M (cid:12) ) are displayed in colored circles, triangles, and squares with 1 σ error bar. A. Prospect of forming supramassive neutron star in double neutron star mergers
The fates of the remnants formed in the double NS mergers have been extensively examined before. Since the EoS of NSmatter is essentially unknown, the relevant examinations were based on some representative models (e.g., Refs.[39–43]). Nowwe reexamine the general prospect of forming SMNSs in the double NS mergers with some EoS-insensitive relationships.The remnants formed in the mergers of the heaviest double NS binary systems may collapse promptly into black holes.But in most cases the initial remnant should be a di ff erentially rotating very massive NS. The di ff erential rotation may beterminated quickly (in a timescale of approximately 0 . − et al. [44]derived an empirical relation among the critical total gravitational mass of DNSs ( M tot , c ), the mass and compactness of NSin the nonrotation maximum equilibrium configuration (denoted by M TOV and ζ TOV , respectively), the dimensionless angularmomentum of remnant at the onset of collapse ( j c ), and the total mass lost apart from the remnant core ( m loss , including thekilonova / macronova ejecta and the accretion torus / disk), which reads M tot , c ≈ M TOV (1 + . ζ − j + . ζ − j )(0 . + . ζ TOV )(1 − . M − (cid:12) m loss ) + m loss , (5)where j kep is the dimensionless angular momentum of NS rotating at Keplerian velocity.With the empirical relation of j kep ≈ . ζ . , Eq.(5) reduces to M tot , c ≈ M TOV + . (cid:32) j c j kep (cid:33) + . (cid:32) j c j kep (cid:33) (0 . + . ζ TOV )(1 − . M − (cid:12) m loss ) + m loss . (6)The merger remnant, if not supported by the thermal pressure and if it has entered the phase of uniform rotation, will collapse toa black hole as M tot > M TOV (0 . + . ζ TOV )(1 − . M − (cid:12) m loss ) + m loss , (7) M / M l o g ( L / ( m J y k p c )) low-mass grouphigh-mass groupJ0740+6620J1600-3053J1811-2405J1959+2048J2215+5135J0348+0432 low-mass grouphigh-mass group FIG. 4: Pulsar mass vs radio luminosity. The massive neutron stars (which may be heavier than 1.6 M (cid:12) ) are displayed in magenta crosses with1 σ error bar. The radio luminosities were measured at 1.4 (most events) or 0.35GHz (a few objects) [32] except J0348 + − . where M tot is the total gravitational mass of the double NS system. A stable massive NS will be the output if instead we have M tot < M TOV (0 . + . ζ TOV )(1 − . M − (cid:12) m loss ) + m loss . (8)For the rest, a SMNS remains as long as the decreasing angular momentum meets the condition of j ≥ j c , where j c = (cid:113) − . + (cid:112) . + . C − ( M tot − m loss ) − . j kep (9)and C = M TOV (0 . + . ζ TOV )(1 − . M − (cid:12) m loss ).In addition to j c / j kep and m loss , M TOV and R TOV (note that ζ TOV ≡ GM TOV / R TOV c ) play some roles in shaping M tot , c . Thesetwo parameters are still unknown yet. Anyhow reasonable constraints have been set by the jointed data analysis of gravitationalwave, NS observations, and some nuclear experiments. We take the very recent constraints obtained in the analysis of PSRJ0030 + M TOV is restricted in the range of (2 . , . M (cid:12) . The two-dimensional probability distribution is shown in Fig.5, with which it is straightforward to predict M tot , c (more exactly, the range)for the given M TOV , m loss , and j / j kep . The results are shown in Fig.6. Clearly, for some light double neutron star binary systems(note that about half of the Galactic binary NS systems have a mass less than or equal to 2 . M (cid:12) , as listed in TableI), the mergerremnant should be SMNSs as long as j / j kep ≥ . M max ≥ . M (cid:12) (note that j ≥ . j kep unless the postmerger gravitationalwave radiation is more e ffi cient than that found in the literature [44, 46–48]), in agreement with previous research (e.g., Refs.[39, 41–43]).There is, however, one caution that Eqs.(5) and (6) are derived with some empirical relationships that are established for agroup of representative EoSs of NSs at zero temperature. For the nascent NSs formed in the mergers, the typical temperature isfound to be approximately 30 −
50 MeV (e.g., Ref.[40]). At such high temperatures, strong interactions may sizably enhance thenumber density of negatively charged pions. Very recently, Fore and Reddy [50] have calculated such an e ff ect for a range ofdensity (below 1.4 times of the nuclear saturation density) and temperature using the virial expansion and found that the thermalpions increase the proton fraction and soften the EoS. Therefore a new parameter R ≡ M TOV ( T ) / M TOV ≤ potential reduction of the “e ff ective” maximum gravitational mass at the high temperature ( T is approximately tensMeV; note that if it finally turns out that R > T O V Best fit:
TOV = 0.081 M TOV + 0.109Data points2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 M TOV / M f i t T O V / T O V [ % ] FIG. 5: M TOV − ζ TOV correlation inferred from the joint constraints set by the data of PSR J0030 + σ , 2 σ , and 3 σ contours, respectively. promising than estimated in this work). In reality, R may be a complicated function. For illustration we simply assume a constantand then discuss the range of R that can be probed in the near future. As demonstrated in Fig.7, in principle, the mergers ofsome light double NS binary systems can e ff ectively probe R in a wide range supposing the remnant nature (hypermassive NS,SMNS, or even stable NS) can be reasonably inferred from the follow-up electromagnetic observations. Therefore, the possiblee ff ect of high temperature on softening the EoS can be unambiguously clarified. B. Energy reservoir of SMNSs and the acceleration of EeV cosmic rays
The SMNSs, if formed, have a huge amount of energy that could give rise to some very interesting phenomena. In thissubsection, we estimate the kinetic rotational energy of these objects and discuss the possible role in accelerating EeV cosmicrays. In the Appendix of [44], we have introduced the details of deriving the empirical relation among M crit , M TOV , and j / j kep .Our current approach is basically the same as [44] except that the EoS sample has been further expanded to include also STOS[51], H4 [52], MS1B, MS1, MS1B PP, MS1 PP [53], AP3, APR [54–56], SKI5 [57–59], HS TM1, HS TMA, HS NL3 [60–64].We aim to get an empirical relationship among E rot , M TOV , R TOV and j / j kep . For our purpose, first, we examine the rotationalkinetic energy of the SMNS at the mass shedding limit. Motivated by the facts of E rot , kep ∝ I kep Ω , I kep ∝ ζ kep M kep R [65], M kep ∝ M TOV , R kep ∝ R TOV and Ω kep ∝ (cid:113) M TOV / R [10], we have E rot , kep ∝ ζ M TOV . The polynomial fit to the numericalresults of a set of widely discussed EOSs yields (see the left panel of Fig.8) E rot , kep ≈ erg [13 . ζ M TOV / M (cid:12) ) − . . (10)In the right panel of Fig.8, we show the results of E rot / E rot , kep as a function of j / j kep , which reads E rot ( j ) = E rot , kep (cid:20) . (cid:16) j / j kep (cid:17) + . (cid:16) j / j kep (cid:17) (cid:21) , (11)with which we can estimate the total amount of kinetic energy lost before the collapse of the SMNS if its initial (i.e., at the birth, j int ) and final (i.e., at the collapse, j c ) values of the dimensionless angular momentum j are known. The fit for j / j kep < . M TOV / M M t o t , c / M j / j kep = 0.0 j / j kep = 0.4 j / j kep = 0.6 j / j kep = 0.8 j / j kep = 1.0 FIG. 6: Critical total gravitational mass of DNSs vs M TOV , supposing M loss = . M (cid:12) . The vertical cyan region represents our evaluated masscuto ff in the NS mass distribution (see Fig.1). The horizontal gray regions, from the top to the bottom, represent the total gravitational mass ofGW190425, GW170817, and PSR J0514 − relatively poor, but it is still acceptable since the kinetic energy is already small in such a range. It is challenging to reliablyinfer j int with the current information from the observations and the numerical simulations (see Refs.[44, 66] for the extendeddiscussion). Anyhow, the dominant contribution of the early time (in the first approximately 100ms or so) angular momentumloss is likely due to the postmerger gravitational wave radiation. If such a process just carries away the energy of a few percentof solar mass (i.e., in the low end found in Ref.[67]), one would have j int ≈ j kep (see e.g., Refs.[42, 46]).After the formation of SMNS (or stable NS), the quadrupole gravitational wave radiation is approximately10 ( (cid:15)/ − ) ( P int / − erg s − (note that the upper limit on the ellipticity set by LIGO / Virgio observations for the Galac-tic pulsars is (cid:15) ≤ − [68]. Then, usually the energy release would be dominated by magnetic dipole radiation and the spindownluminosity can be estimated [69] by L dip ≈ erg s − (cid:32) R s cm (cid:33) (cid:32) B ⊥ G (cid:33) (cid:32) P int (cid:33) − , (12)where R s is the radius of SMNS, B ⊥ = B s sin α , P int is the initial spin period of the SMNS, and B s is the magnetic field onthe surface and α is the angle between the spin and dipole axes. Then the dipole radiation energy injected into the blast waveand ejecta is dE inj / dt ≈ L dip (1 + t /τ ) − , where τ = E rot , int / L dip is the so-called spindown timescale. For quickly rotating SMNS,supposing R s ∼ P int ∼ . τ ≈ × s ( E rot , int / erg)( R s /
14 km) − ( B ⊥ / G) − ( P int / . . Pleasebear in mind that at t ≈ τ , the total energy released is E inj ≈ E rot , int /
4, before which j may have already dropped below j c andthe SMNS has collapsed. In this work, for illustration, we ignore such a possibility and simply estimate the possible EeV cosmicray accelerated in the SMNS wind-driving blast wave (see also Ref.[70], in which the discussions, however, are mixed with otherpossible models of the fast radio bursts and are hence not systematic). If the frequency ω wind = / P int ≈ P int / . − of the electromagnetic wind of the SMNS is lower than the plasma frequency ω plasma ≈ . × Hz n / , the wind would beabsorbed by the surrounding material (either the shocked circumburst medium or the merger-driving subrelativistic outflow).Such a condition can be rewritten as n e , c ≥ × − cm − .If the absorption is dominated by the merger-driving subrelativistic outflow (i.e., the circumburst medium has a lownumber density n m (cid:28) n e , c ) with a mass of M ej ∼ . M (cid:12) , we would have R c ≤ . × cm ( M ej / . M (cid:12) ) / andthe bulk Lorentz factor of the outflow will peak with Γ τ ≈ + ( E inj / erg)( M ej / . M (cid:12) ) − , supposing τ < R c / c ∼ . × s ( M ej / . M (cid:12) ) / (i.e., B ⊥ ≥ G). The deceleration radius can be estimated to be R dec ∼ ( M ej / Γ τ n m m p ) / ∼ M TOV ( M ) = M T O V ( T ) / M T O V M tot = 2.7MM tot = 2.5MM tot = 2.34M FIG. 7: The “hypothesized” R to collapse for the given M tot (in particular for the lightest NSs or double NS binaries) and the zero temperature M TOV . The dashed and dotted lines are for j / j kep = .
8, respectively. A minimum M tot = . M (cid:12) is adopted because the lightest neutronstar detected so far has an accurately measured mass of 1 . M (cid:12) [49]. For simplicity, the best-fit ζ TOV − M TOV relation presented in Fig.5 isadopted. E r o t / e r g AP3AP4BSK20BSK21ENGH4MPA1MS1BMS1SLYWFF1WFF2SLY4APR APR4_EPPMS1B_PPMS1_PPRSSK255SK272SKASKBSKI2SKI3SKI4SKI5SKI6SKMP SLY2SLY230ASLY9HQC18SFHXHS_NL3HS_DD2STOSLS220HS_TMAHS_TM1DD2YSFHo M TOV / M | E f i t r o t / E r o t | [ % ] E : = E r o t / E m a x AP3AP4BSK20BSK21ENGH4MPA1MS1BMS1SLYWFF1WFF2SLY4APR APR4_EPPMS1B_PPMS1_PPRSSK255SK272SKASKBSKI2SKI3SKI4SKI5SKI6SKMP SLY2SLY230ASLY9HQC18SFHXHS_NL3HS_DD2STOSLS220HS_TMAHS_TM1DD2YSFHo j / j kep | E f i t / E | [ % ] FIG. 8: Left panel: kinetic rotational energy of NS with mass shedding configuration. Right panel: kinetic rotational energy of NS at thecritical point of certain angular momentum. . × cm ( M ej / . M (cid:12) ) / ( Γ τ / − / ( n m / − cm − ) − / , where m p is the rest mass of proton. The maximum energy of theaccelerated particles can be estimated as ε CR , max ∼ β ZeB dec R dec ∼ × eV Z β Γ (cid:18) n m − cm − (cid:19) (cid:32) M ej . M (cid:12) (cid:33) (cid:18) (cid:15) B − (cid:19) , (13)where B dec ∼ . × − G β Γ ( n m / − cm − ) / ( (cid:15) B / . / is the strength of shock amplified magnetic field at R dec . In thecase of the binary NS mergers taking place in the relatively dense circumburst medium (i.e., n m ∼ . − ), the deceleration of1the blast wave starts at the radius of R dec ∼ ( M ej / Γ τ n m m p ) / ∼ . × cm ( M ej / . M (cid:12) ) / ( Γ τ / − / ( n m / . − ) − / , andmost of the kinetic energy of the SMNS would have been injected into the blast wave supposing B ⊥ ≥ × G. In this case,e ffi cient EeV protons are still accelerated. For B ⊥ (cid:28) × G, the significant energy injection lasts much longer time and theaccelerated protons can only reach sub-EeV energy region. To widely explore the various possibilities, we present the numericalcalculation results in Fig.9. Clearly, for B ⊥ > G, EeV cosmic-ray protons are indeed the plausible outputs. t/ C R , m a x ( E e V ) M ej = 0.05 M , V ej = 0.2 c , R s = 14 km , B = 0.01 B = 10 G B = 10 G B = 10 G B = 10 G FIG. 9: Capability of the shocks driven by the kilonova ejecta and the possible energy injection from the central SMNS. The solid, dotted,dashed, and dot-dashed lines are for B s = , , , G, respectively, while the colors of the lines (black, blue, red) correspond tothe cases of n m = (10 − , − , − ) cm − respectively. For n m = − cm − , the red solid line terminates at t /τ ∼ .
03 because at later timeswe have R > . × cm, for which the plasma frequency is lower than the frequency of pulsar wind and the pulsar wind will escape freelyrather than be absorbed by the outmoving ejecta. The EeV cosmic-ray proton flux accelerated by the double NS merger formed SMNSs at approximately 10 eV is estimatedas F EeV − CR ∼ − η (cid:18) ω . (cid:19) (cid:32) E inj erg (cid:33) × (cid:32) R DNS , SMNS yr − Gpc − (cid:33) m − s − sr − eV − , (14)which is consistent with the observed cosmic-ray flux F obs ( ε CR ) = C ( ε CR / . × eV) − . ± . with C = (9 . ± . × − m − s − sr − eV − [71], for the EeV cosmic-ray acceleration e ffi ciency η ∼ .
03 and ω ∼ .
1, the fraction of total cosmic-ray energy at each energy decade. Here we normalize the SMNS formation rate to that of the double NS mergers [20] becauseSMNSs may be produced in a good fraction of such a kind of gravitational wave events (see Fig.6). A good fraction of doubleNS mergers is expected to take place in the elliptical galaxies that are short of the star formation now [72]. If some EeV cosmicrays or PeV(10 eV) neutrinos can be found in the directions of the elliptical galaxies, our arguments will be strongly favored. IV. CONCLUSIONS
The mass distribution of the neutron stars is known to be important in shedding light on the supernova explosion and theaccretion dynamics of binary neutron star. Moreover, the inferred cuto ff mass may represent the maximum mass of the nonro-tating cold neutron stars, which is essential for revealing the equation of state of matter with ultrahigh density. However, themass measurements are challenging and the sample just consists of 74 objects in the work by Alsing et al. [11]. Thanks to thededicated worldwide joint e ff orts, recently, the sample of neutron stars with a measured mass has been growing very quickly.In this work, we collect the mass measurements of NSs from the latest literature (including the updates of the masses of some2“old” objects) and increase the total number of neutron stars in the sample to 103. With this new sample, we adopt a flexibletwo-component mixture model and a Gaussian plus Cauchy-Lorentz component model to infer the mass distribution of neutronstars and use Bayesian model selection to explore evidence for multimodality and a sharp cuto ff in the mass distribution. Theresults for these two models are consistent with each other (see Fig.1). In agreement with previous studies, we find evidence for abimodal distribution together with a cuto ff at a mass of M max = . + . − . M (cid:12) (68% credible interval; for the 95% credible interval, M max = . + . − . M (cid:12) ). Our M max is larger than that reported by Alsing et al. [11] by approximately 0 . M (cid:12) mainly because of theinclusion of the recent measurements of two massive objects PSR J1959 + + + . + . − . M (cid:12) [15], with which we have M max = . + . − . M (cid:12) ). Compared toprevious works, our method faithfully reproduces the asymmetric mass measurement errors and avoids the bias caused by theapproximation of Gaussian measurement errors. Our resulting distributions of µ , σ , and M max are less a ff ected by the priorsand NS population models, indicating the robustness of the approaches. We have discussed some possible selection e ff orts, forexample that our sample mainly consists of the neutron stars found in binary systems and the most accurately measured objectsare from the double neutron star systems. However, the most massive neutron stars are heavily recycled and the isolated neutronstars are not expected to have a mass close to M max . We also show that both PSR J0348 + + M max , are among the luminous radio pulsars. Hence at least for the current sample, thereis no evidence for a sizeable nondetecion / misidentification probability of the very massive NSs in radio (these recycled pulsarsrotate very quickly, which generate strong radio radiation though the surface magnetic fields are low). In agreement with Alsing et al. [11], we suggest that the selection e ff ects are unlikely to account for the inferred hard cuto ff at M max . We further examinedthe possible dependence between radio luminosity and mass but did not find any evidence. Together with the previous findings[5, 35], we suggest that the selection e ff ects will not introduce serious bias to the observed mass distribution (see also Ref.[5]).If the inferred cuto ff mass M max represents the maximum gravitational mass of nonrotating cold neutron stars ( M TOV ), theprospect of forming supramassive remnants is found to be promising for the double neutron star mergers with a total gravitationalmass less than or equal to 2 . M (cid:12) (see Fig.6) unless the thermal pions could substantially soften the equation of state for the veryhot neutron star matter or alternatively the postmerger gravitational wave radiation has carried away the kinetic rotational energyof approximately 0 . M (cid:12) (for which j / j kep < . M tot = . M (cid:12) ) and PSR J1946 + M tot = . M (cid:12) ), may e ff ectively probe the potential e ff ectof EoS softening by the thermal pions generated at high temperatures, supposing the remnant nature (hypermassive NS, SMNSor even stable NS) can be reliably inferred from the follow-up electromagnetic observations. The SMNSs are expected to havea typical kinetic rotational energy of approximately 1 − × ergs. If not radiated mainly in a gravitational wave, thanksto a high neutron star merger rate of approximately 10 Gpc − yr − and the plausible high chance of forming SMNSs in suchmergers, the neutron star mergers are likely the significant sources of EeV cosmic-ray protons.In the O3 run of advanced LIGO / Virgo, almost all double neutron star event / candidates were just detected by one ofthe two LIGO detectors (https: // gracedb.ligo.org / superevents / public / O3 / ). These events were poorly localized and no elec-tromagnetic counterparts were reliably identified. The situation will change considerably in the late observing runs (seehttps: // dcc.ligo.org / public / / P1900218 / / SummaryForObservers.pdf for the latest schedule for the future plans of the sec-ond generation gravitational wave detectors). The Kamioka Gravitational Wave Detector (KAGRA) had already joined the O3run in March 2020. Moreover, the sensitivities of Virgo and KAGRA detectors will be enhanced by a factor of a few in the O4run that is expected to start in January 2022. LIGO-India is anticipated to join the O5 run in 2025. So the future gravitationalwave events will be significantly better localized and their electromagnetic counterparts, in particular the kilonovae / macronovae,will be much more frequently detected. These events are expected to be able to test some speculations of this work, in particularthe possibilities of R <
1, and the formation of SMNSs in a good fraction of neutron star mergers.
Acknowledgments
We thank the referees and Professor C. M. Zhang for very helpful suggestions. This work was supported in part by NSFCunder Grants No. 11525313 (i.e., Funds for Distinguished Young Scholars) and No. 11921003.
Appendix A: Neutron star data
Here, we summarize the measurements of the neutron star masses, totally 103 items classified into six types. This sample isupdated to February 2020.3
TABLE I: Catalog of neutron stars with mass measurementsName Type m p ( M (cid:12) ) f ( M (cid:12) ) m T ( M (cid:12) ) q ReferenceB1534 +
12 NS-NS 1.3330 ± et al. [73]B1534 +
12 comp. NS-NS 1.3455 ± et al. [73]B1913 +
16 NS-NS 1.4398 ± et al. [74]B1913 +
16 comp. NS-NS 1.3886 ± et al. [74]B2127 +
11C NS-NS 1.358 ± et al. [75]B2127 +
11C comp. NS-NS 1.354 ± et al. [75]J0453 + ± et al. [49]J0453 + ± et al. [49]J0509 + ± et al. [76]J0509 + ± et al. [76]J0514-4002A NS-NS 1 . + . − . Rifolfi et al. [77]J0514-4002A comp. NS-NS 1 . + . − . Rifolfi et al. [77]J0737-3039A NS-NS 1.3381 ± et al. [78]J0737-3039B NS-NS 1.2489 ± et al. [78]J1756-2251 NS-NS 1.341 ± et al. [79]J1756-2251 comp. NS-NS 1.230 ± et al. [79]J1757-1854 NS-NS 1.3384 ± et al. [80]J1757-1854 comp. NS-NS 1.3946 ± et al. [80]J1807-2500B NS-NS 1.3655 ± et al. [37]J1807-2500B comp. NS-NS 1.2064 ± et al. [37]J1906 + ± et al. [81]J1906 + ± et al. [81]GW170817A NS-NS 1 . + . − . Abbott et al. [20]GW170817B NS-NS 1 . + . − . Abbott et al. [20]GW190425A NS-NS 1 . + . − . Abbott et al. [23]GW190425B NS-NS 1 . + . − . Abbott et al. [23]J1411 + ± et al. [82]J1518 + ± et al. [83]J1811-1736 NS-NS 0.128121 2.57 ± et al. [84]J1829 + ± et al. [85]J1913 + ± et al. [86]J1930-1852 NS-NS 0.34690765 2.54 ± et al. [87]J1946 + ± et al. [88]B1855 +
09 NS-WD 1.37 ± et al. [12]J0337 + ± et al. [89]J0348 + ± et al. [90]J0437-4715 NS-WD 1.44 ± et al. [91]J0621 + . + . − . Kasian L. E. [92]J0740 + . + . − . Cromartie et al. [13]J0751 + ± et al. [93]J1012 + ± et al. [6]J1141-6545 NS-WD 1.27 ± et al. [94]J1600-3053 NS-WD 2 . + . − . Arzoumanian et al. [12]J1614-2230 NS-WD 1.908 ± et al. [12]J1713 + ± et al. [95]J1738 + ± et al. [96]J1741 + . + . − . Arzoumanian et al. [12]J1748-2446am NS-WD 1 . + . − . Andersen and Ransom [97]J1802-2124 NS-WD 1.24 ± et al. [98]J1811-2405 NS-WD 2 . + . − . Ng et al. [16]J1909-3744 NS-WD 1.48 ± et al. [12]J1911-5958A NS-WD 1.34 ± et al. [99]J1918-0642 NS-WD 1.29 ± et al. [12]J1946 + ± et al. , [100]J1949 + . + . − . Zhu et al. [101]J1950 + ± et al. [101]J1959 + ± + ± et al. [12]J2045 + ± et al. [102]Table I ( Continued on next page ) Table I (
Continued )Name Type m p [ M (cid:12) ] f [ M (cid:12) ] m T [ M (cid:12) ] q ReferenceJ2053 + ± et al. [102]J2215 + . + . − . Kandel and Romani [14]J2222-0137 NS-WD 1.76 ± et al. [103]J2234 + . + . − . Stovall et al. [104]B1516 +
02B NS-WD 0.000646723 2.29 ± et al. [105]B1802-07 NS-WD 0.00945034 1.62 ± +
46 NS-WD 0.246261924525 2.64 ± ± et al. [107]J1748-2446I NS-WD 0.003658 2.17 ± et al. [108]J1748-2446J NS-WD 0.013066 2.20 ± et al. [108]J1750-37A NS-WD 0.0518649 1.97 ± et al. [109]J1824-2452C NS-WD 0.006553 1.616 ± et al. [105]NGC6440B NS-WD 0.0002266235 2.69 ± ff ord and Ransom [110]J1311-3430 NS-WD(?) 3 × − ± et al. [111]J1723-2837 NS-WD(?) 0.005221 3.45 ± ± et al. [113]J1816 + ± et al. [114]J0045-7319 NS-MS 1.58 ± + ± et al. [116]J1903 + ± et al. [12]J0030 + . + . − . Riley et al. [24](NICER)4U1538-522 HMXB 1.02 ± et al. [117]4U1700-377 HMXB 1.96 ± et al. [117]Cen X-3 HMXB 1.57 ± et al. [117]EXO 1722-363 HMXB 1.91 ± et al. [117]Her X-1 HMXB 1.07 ± et al. [118]J013236.7 + ± et al. [119]LMC X-4 HMXB 1.57 ± et al. [117]OAO 1657-415 HMXB 1.74 ± et al. [117]SAX J1802.7-2017 HMXB 1.57 ± et al. [117]SMC X-1 HMXB 1.21 ± et al. [117]Vela X-1 HMXB 2.12 ± et al. [117]XTE J1855-026 HMXB 1.41 ± et al. [117]2S 0921-630 LMXB 1.44 ± . + . − . ¨Ozel et al. [121]4U1702-429 LMXB 1.9 ± et al. [122]4U 1724-207 LMXB 1 . + . − . ¨Ozel et al. [121]4U 1820-30 LMXB 1 . + . − . ¨Ozel et al. [121]Cyg X-2 LMXB 1.71 ± et al. [123]KS 1731-260 LMXB 1 . + . − . ¨Ozel et al. [121]EXO 1745-248 LMXB 1 . + . − . ¨Ozel et al. [121]SAX J1748.9-2021 LMXB 1 . + . − . ¨Ozel et al. [121]X 1822-371 LMXB 1.96 ± et al. [124]XTE J2123-058 LMXB 1.53 ± et al. [125]Notes: NS-NS, double neutron star system; NS-WD, neutron star-white dwarf binary; NS-MS, neutron star-main sequence star system;HMXB, high mass x-ray binary; LMXB, low mass x-ray binary; INS, isolated neutron star.The question mark means the nature of the companion is uncertain.[1] A. Hewish, S. J. Bell, J. D. H. Pilkington et al. , Nature , 709 (1968).[2] C. M. Zhang, J. Wang, Y. H. Zhao, H. X. Yin, L. M. Song, D. P. Menezes, D. T. Wickramasinghe, L. Ferrario, and P. Chardonnet,Astron. Astrophys. , A83 (2011).[3] R. Valentim, E. Rangel, and J. E. Horvath, Mon. Not. R. Astron. Soc. , 1427 (2011).[4] F. ¨Ozel, D. Psaltis, R. Narayan, and A. S. Villarreal, Astrophys. J. , 55 (2012).[5] B. Kiziltan, A. Kottas, M. De Yoreo, and S. E. Thorsett, Astrophys. J. , 66 (2013).[6] J. Antoniadis, T. M. Tauris, F. ¨Ozel, E. Barr, D. J. Champion, and P. C. C. Freire, arXiv:1605.01665.[7] J. Horvath and R. Valentim, Handbook of Supernovae (Springer-Verlag, Berlin, 2017). [8] J. L. Friedman, L. Parker, and J. R. Ipser, Astrophys. J. , 115 (1986).[9] Y. Z.Fan, X. F. Wu, and D. M. Wei, Phys. Rev. D , 067304 (2013).[10] C. Breu and L. Rezzolla, Mon. Not. R. Astron. Soc. , 646 (2016).[11] J. Alsing, H. O. Silva, and E. Berti, Mon. Not. R. Astron. Soc. , 1377 (2018).[12] Z. Arzoumanian, A. Brazier, S. Burke-Spolaor et al., Astrophys. J. Suppl. Ser. , 37 (2018).[13] H. T. Cromartie, E. Fonseca, S. M. Ransom et al., Nat. Astron. , 72 (2020).[14] D. Kandel and R. W. Romani, Astrophys. J. , 101 (2020).[15] M. Linares, T. Shahbaz, and J. Casares, Astrophys. J. , 54 (2018).[16] C. Ng, L. Guillemot, P. C. C. Freire, M. Kramer, D. J. Champion, I. Cognard, G. Theureau, and E. D. Barr, Mon. Not. R. Astron.Soc. , 1261 (2020).[17] R. A. Remillard and J. E. McClintock, Ann. Rev. Astron. Astrophys. , 49 (2006).[18] F. ¨Ozel and P. Freire, Ann. Rev. Astron. Astrophys. , 401 (2016).[19] I. H. Stairs, Living Rev. Relativity , 5 (2003).[20] B. P. Abbott, et al. , Phys. Rev. Lett , 161101 (2017).[21] M. Z. Han, S. P. Tang, Y. M. Hu, Y. J. Li, J. L. Jiang, Z. P. Jin, Y. Z. Fan, and D. M. Wei, Astrophys. J. Lett. , L5 (2020).[22] R. J. Foley, D. A. Coulter, C. D. Kilpatrick, A. L. Piro, E. Ramirez-Ruiz, and J. Schwab, Mon. Not. R. Astron. Soc. , 190 (2020).[23] B. P. Abbott, et al. , Astrophys. J. Lett. , L3 (2020).[24] T.E. Riley, A. L. Watts, S. Bogdanov et al., Astrophys. J. Lett. , L21 (2019).[25] M. C. Mille, F. K. Lamb, A. J. Dittmann et al., Astrophys. J. Lett. , L24 (2019).[26] S. P. Tang, J. L. Jiang, W. H. Gao, Y. Z. Fan, and Da-Ming Wei, Astrophys. J. , 45 (2020).[27] N. Farrow, X. J. Zhu, and E. Thrane, Astrophys. J. , 18 (2019).[28] C. Fern´andez and F. Steel, J. Am. Stat. Assoc. , 359 (1998).[29] J. Buchner, PyMultiNest: Python interface for MultiNest, ascl:1606.005, 2016.[30] R. E. Kass and A. E. Raftery, J. Am. Stat. Assoc. , 773 (1995).[31] K. Ebinger, S. Curtis, S. Ghosh, C. Fr¨ohlich, M. Hempel, A. Perego, M. Liebend¨orfer, and F. Thielemann, Astrophys. J. , 91 (2020).[32] R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, Astrophys. J. , 1993 (2005)[33] J. H. Friedman, Multivariate adaptive regression splines. Ann. Stat. , 1 (1991).[34] A. Szary, B. Zhang, G. I. Melikidze, J. Gil, and R.-X. Xu, Astrophys. J. , 59 (2014).[35] D. R. Lorimer, A. J. Faulkner, A. G. Lyne, R. N. Manchester, M. Kramer, M. A. McLaughlin, G. Hobbs, A. Possenti, I. H. Stairs, and F.Camilo, Mon. Not. R. Astron. Soc. , 777 (2006).[36] R. S. Lynch, J. Boyles, S. M. Ransom, I. H. Stairs, and D. R. Lorimer et al., Astrophys. J. , 81 (2013).[37] R. S. Lynch, P. C. C. Freire, S. M. Ransom, and B. A. Jacoby, Astrophys. J. , 109 (2012).[38] S. D. Bates, D. R. Lorimer, and J. P. W. Verbiest, Mon. Not. R. Astron. Soc. , , 1352 (2013).[39] I. A. Morrison, T. W. Baumgarte, and S. L. Shapiro, Astrophys. J. , 941 (2004).[40] K. Hotokezaka, K. Kiuchi, K. Kyutoku, H. Okawa, Y.-i. Sekiguchi, M. Shibata, and K. Taniguchi, Phys. Rev. D , 024001 (2013).[41] S. Lawrence, J.G. Tervala, P.F. Bedaque, and M.C. Miller, Astrophys. J. , 186 (2015).[42] A. L. Piro, B. Giacomazzo, and R. Perna, Astrophys. J. Lett. , L19 (2017).[43] P.X. Ma, J.L. Jiang, H. Wang, Z.P. Jin, Y.Z. Fan, and D.M. Wei, Astrophys. J. , 74 (2018).[44] D. S. Shao, S. P. Tang, S. Xin, J. L. Jiang, Y. Z. Wang, Z. P. Jin, Y. Z. Fan, and D. M. Wei, Phys. Rev. D , 063029 (2020).[45] J. L. Jiang, S. P. Tang, Y. Z. Wang, Y. Z. Fan, and D. M. Wei, Astrophys. J. , 55 (2020).[46] F. Zappa, S. Bernuzzi, D. Radice, A. Perego, and T. Dietrich, Phys. Rev. Lett , 111101 (2018).[47] D. Radice, A. Perego, S. Bernuzzi, and B. Zhang, Mon. Not. R. Astron. Soc. , 3670 (2018).[48] Y. Z. Fan, J. L. Jiang, S. P. Tang, Z. P. Jin, and D. M. Wei, arXiv:2005.10482.[49] J. G. Martinez, K. Stovall, P. C. C. Freire, J. S. Deneva, F. A. Jenet, M. A. McLaughlin, M. Bagchi, S. D. Bates, and A. Ridolfi, Astrophys.J. , 143 (2015).[50] B. Fore and S. Reddy, Phys. Rev. C , 035809 (2020).[51] H. Shen, H. Toki, K. Oyamatsu, and K Sumiyoshi, Nucl. Phys. A637 , 435 (1998).[52] B. D. Lackey, M. Nayyar, and B. J. Owen, Phys. Rev. D , 024021 (2006).[53] H. M¨uller and B. D. Serot, Nucl. Phys. A606 , 508 (1996).[54] G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. , 299 (1971).[55] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C , 1804 (1998).[56] F. Douchin and P. Haensel, Astron. Astrophys. , 151 (2001).[57] P. G. Reinhard and H. Flocard, Nucl. Phys. A584 , 467 (1995).[58] P. Danielewicz and J. Lee, Nucl. Phys.
A818 , 36 (2009).[59] F. Gulminelli and A. R. Raduta, Phys. Rev. C , 055803 (2015).[60] Y. Sugahara and H. Toki, Nucl. Phys. A579 , 557 (1994).[61] H. Toki, D. Hirata, Y. Sugahara, K Sumiyoshi, and I Tanihata, Nucl. Phys.
A588 , 357 (1995).[62] G. A. Lalazissis, J. K¨onig, and P. Ring, Phys. Rev. C , 540 (1997).[63] L. Geng, H. Toki, and J. Meng, Prog. Theor. Phys. , 785 (2005).[64] M. Hempel and J. Scha ff ner-Bielich, Nucl. Phys. A837 , 210 (2010).[65] P. S. Koliogiannis and C. C. Moustakidis, Phys. Rev. C , 015805 (2020).[66] M. Shibata, E. Zhou, K. Kiuchi, and S. Fujibayashi, Phys. Rev. D , 023015 (2019).[67] S. Bernuzzi, D. Radice, C. D. Ott, L. F. Roberts, P. M¨osta, and F. Galeazzi, Phys. Rev. D , 024023 (2016).[68] B. P. Abbott, et al. , Phys. Rev. D , 12 (2019). [69] S. L. Shapiro and S. A. Teukolsky, Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects (Wiley, New York,1983).[70] X. Li, B. Zhou, H.-N. He, Y.-Z. Fan, and D.-M. Wei, Astrophys. J. , 33 (2014).[71] M. Nagano and A. A. Watson, Rev. Mod. Phys. , 689 (2000).[72] E. Berger, Annu. Rev. Astron. Astrophys. , 43 (2014).[73] E. Fonseca, I. H. Stairs, and S. E. Thorsett, Astrophys. J. , 82 (2014).[74] J. M. Weisberg, D. J. Nice, and J. H. Taylor, Astrophys. J. , 1030 (2010).[75] B. A. Jacoby, P. B. Cameron, F. A. Jenet, S. B. Anderson, R. N. Murty, and S. R. Kulkarni, Astrophys. J. , L113 (2006).[76] R. S. Lynch, J. K. Swiggum, V. I. Kondratiev et al., Astrophys. J. , 93 (2018).[77] A. Ridolfi, P. C. C. Freire, Y. Gupta and S. M. Ransom, Mon. Not. R. Astron. Soc. , 3860 (2019).[78] M. Kramer, I. H. Stairs, R. N. Manchester et al., Science , 97 (2006).[79] R. D. Ferdman, I. H. Stairs, M. Kramer et al., Mon. Not. R. Astron. Soc. , 2183 (2014).[80] A. D. Cameron, D. J. Champion, M. Kramer et al., Mon. Not. R. Astron. Soc. , L57 (2018).[81] J. van Leeuwen, L. Kasian, I. H. Stairs et al., Astrophys. J. , 118 (2015).[82] J. G. Martinez, K. Stovall, P.C.C. Freire, J.S. Deneva, T.M. Tauris, A. Ridolfi, N. Wex, F.A. Jenet, M.A. McLaughlin, and M. Bagchi,Astrophys. J. Lett. , L29 (2017).[83] G. H. Janssen, B. W. Stappers, M. Kramer, D. J. Nice, A. Jessner, I. Cognard, and M. B. Purver, Astron. Astrophys. , 753 (2008).[84] A. Corongiu, M. Kramer, B. W. Stappers, A. G. Lyne, A. Jessner, A. Possenti, N. D’Amico, and O. L¨ohmer, Astron. Astrophys. , 703(2007).[85] D. J. Champion, D. R. Lorimer, M. A. McLaughlin, K. M. Xilouris, Z. Arzoumanian, P. C. C. Freire, A. N. Lommen, J. M. Cordes, and F.Camilo, Mon. Not. R. Astron. Soc. , 929 (2005).[86] P. Lazarus, P. C. C. Freire, B. Allen, C. Aulbert, O. Bock, S. Bogdanov, A. Brazier, F. Camilo, F. Cardoso, and S. Chatterjee, Astrophys.J. , 150 (2016).[87] J. K. Swiggum, R. Rosen, M. A. McLaughlin, D. R. Lorimer, S. Heatherly, R. Lynch, S. Scoles, T. Hockett, E. Filik, and J. A. Marlowe,Astrophys. J. , 156 (2015).[88] K. Stovall, P. C. C. Freire, S. Chatterjee, P. B. Demorest, D. R. Lorimer, M. A. McLaughlin, N. Pol, J. van Leeuwen, R. S. Wharton, andB. Allen, Astrophys. J. Lett. , L22 (2018).[89] A. M. Archibald, N. V. Gusinskaia, J. W. T. Hessels, A. T. Deller, D. L. Kaplan, D. R. Lorimer, R. S. Lynch, S. M. Ransom, and I. H.Stairs, Nature , 73 (2018).[90] J. Antoniadis, P. C. C. Freire, N. Wex et al., Science , 448 (2013).[91] D. J. Reardon, G. Hobbs, W. Coles et al., Mon. Not. R. Astron. Soc. , 1751 (2016).[92] L. E. Kasian, Ph.D. thesis, University of British Columbia, 2012.[93] G. Desvignes, R. N. Caballero, L. Lentati et al., Mon. Not. R. Astron. Soc. , 3341 (2016).[94] V. V. Krishnan, M. Bailes, W. van Straten et al., Science , 577 (2020).[95] W. W. Zhu, G Desvignes, N Wex et al., Mon. Not. R. Astron. Soc. , 3249 (2019).[96] J. Antoniadis, M. H. van Kerkwijk, D. Koester, P. C. C. Freire, N. Wex, T. M. Tauris, M. Kramer, and C. G. Bassa, Mon. Not. R. Astron.Soc. , 3316 (2012).[97] B. C. Andersen and S. M. Ransom, Astrophys. J. , L13 (2018).[98] R. D. Ferdman, I. H. Stairs, M. Kramer et al., Astrophys. J. , 764 (2010).[99] C. G. Bassa, M. H. van Kerkwijk, D. Koester, and F. Verbunt, Astron. Astrophys. , 295 (2006).[100] E. D. Barr, P. C. C. Freire, M. Kramer, D. J. Champion, M. Berezina, C. G. Bassa, A. G. Lyne, and B. W. Stappers, Mon. Not. R. Astron.Soc. , 1711 (2017).[101] W. W. Zhu, P. C. C. Freire, B. Knispe et al., Astrophys. J. , 165 (2019).[102] M. Berezina, D. J. Champion, P. C. C. Freire et al., Mon. Not. R. Astron. Soc. , 4421 (2017).[103] I. Cognard, P. C. C. Freire, L. Guillemot et al., Astrophys. J. , 128 (2017).[104] K. Stovall, P. C. C. Freire, J. Antoniadis et al., Astrophys. J. , 74 (2019).[105] P. C. C. Freire, A. Wolszczan, M. van den Berg, and J. W. T. Hessels, Astrophys. J. , 1433 (2008).[106] S. E. Thorsett and D. Chakrabarty, Astrophys. J. , 288 (1999).[107] P. C. C. Freire, A. Ridolfi, M. Kramer et al. , Mon. Not. R. Astron. Soc. , 857 (2017).[108] S. M. Ransom, J. W. T. Hessels, I. H. Stairs, P. C. C. Freire, F. Camilo, V. M. Kaspi, and D. L. Kaplan, Science , 892 (2005).[109] P. C. C. Freire, S. M. Ransom, S. Begin, I. H. Stairs, J. W. T. Hessels, L. H. Frey, and F. Camilo, Astrophys. J. , 670 (2008).[110] N. Cli ff ord and S. M. Ransom, B.S. thesis, University of Virginia, 2019.[111] R. W. Romani, A. V. Filippenko, J. M. Silverman, S. B. Cenko, J. Greiner, A. Rau, J. Elliott, and J. P. Holger, Astrophys. J. , L36(2012).[112] A. D. van Staden and J. Antoniadis, Astrophys. J. Lett. , L12 (2016).[113] F. R. Ferraro, E. Sabbi, R. Gratton, A. Possenti, N. D’Amico, A. Bragaglia, and F. Camilo, Astrophys. J. Lett. , L13 (2003).[114] D. L. Kaplan, V. B. Bhalerao, M. H. van Kerkwijk, D. Koester, S. R. Kulkarni, and K. Stovall, Astrophys. J. , 158 (2013).[115] D. J. Nice, in IAU Proceedingss (Cambridge University Press, Cambridge, 2003).[116] A. T. Deller, A. M. Archibald, W. F. Brisken et al., Astrophys. J. , L25 (2012).[117] M. Falanga, E. Bozzo, A. A. Lutovinov, J. M. Bonnetbidaud, Y. Fetisova, and J. Puls, Astron. Astrophys. , A130 (2015).[118] M. L. Rawls, J. A. Orosz, J. E. McClintock, M. A. P. Torres, C. D. Bailyn, and M. M. Buxton, Astrophys. J. , 25 (2011).[119] V. Bhalerao, M. H. van Kerkwijk, and F. Harrison, Astrophys. J. , 10 (2012).[120] D. Steeghs and P.G. Jonker, Astrophys. J. , L85 (2007).[121] F. ¨Ozel, D. Psaltis, G. Tolga, G. Baym, C. Heinke, and S. Guillot, Astrophys. J. , 28 (2016). [122] J. N¨attil¨a, M. C. Miller, A. W. Steiner, J. J. E. Kajava, V. F. Suleimanov, and J. Poutanen, Astron. Astrophys. , A31 (2017).[123] J. Casares, J. I. Gonz´alez Hern´andez, G. Israelian, and R. Rebolo, Mon. Not. R. Astron. Soc. , 2517 (2010).[124] T. Munoz-Darias, J. Casares, and I. G. Martinez-Pais, Astrophys. J. , 502 (2005).[125] D. M. Gelino, J. A. Tomsick, and W. A. Heindl, Bull. Am. Astron. Soc.34