Measuring Mass and Radius of the Maximum-mass Nonrotating Neutron Star

The mass (M TOV) and radius (R TOV) of the maximum-mass nonrotating neutron star (NS) play a crucial role in constraining the elusive equation of state of cold dense matter and in predicting the fate of remnants from binary neutron star (BNS) mergers. In this study, we introduce a novel method to deduce these parameters by examining the mergers of second-generation (2G) black holes (BHs) with NSs. These 2G BHs are assumed to originate from supramassive neutron stars (SMNSs) formed in BNS mergers. Since the properties of the remnant BHs arising from the collapse of SMNSs follow a universal relation governed by M TOV and R TOV, we anticipate that by analyzing a series (∼100 detections) of mass and spin measurements of the 2G BHs using the third-generation ground-based gravitational-wave detectors, M TOV and R TOV can be determined with a precision of ∼0.01M ⊙ and ∼0.6 km, respectively.


INTRODUCTION
The equation of state (EOS) of cold dense matter holds significant importance and is currently undetermined.Observations of neutron stars (NSs), such as their masses and radii, can be utilized to constrain the unknown EOS.The theoretical maximum mass (M TOV ) that a nonrotating NS can sustain is a fundamental characteristic of an EOS.The observed maximum mass of NSs through pulsar timing (e.g., PSR J0740+6620, Fonseca et al. 2021) or the maximum cutoff mass inferred from NS population (Alsing et al. 2018;Shao et al. 2020a;Fan et al. 2023) serve as a stringent lower limit for M TOV , which can effectively eliminate various proposed EOS models.While the (multimessenger) analyses of the NS merger-driven events, in particular GW170817/GRB 170817A/AT2017gfo, have yielded a rough estimate of M TOV (Fan et al. 2013;Shao et al. 2020b) or an upper limit (e.g., Margalit & Metzger 2017;Shibata et al. 2017;Rezzolla et al. 2018;Ruiz et al. 2018) that rules out very stiff EOS.If one has a series of mass-radius (or tidal deformability) measurements of NSs, the EOS can be determined by solving the so-called inverse stellar problem (Lindblom & Indik 2014). 1 For some NSs within low-mass X-ray binary (LMXB) systems, their mass and radius can be measured through spectroscopic observations of thermonuclear (Type-I) X-ray bursts exhibiting photospheric radius expansion, or by observing the angular size of the NSs during periods of quiescent states (see, e.g., Özel & Freire 2016, for a comprehensive review).Additionally, mass and radius information can also be derived from X-ray timing observations using a novel pulse profile modeling of the hot-spot emission from the surfaces of rotating NSs (Bogdanov et al. 2019), a method considered more reliable than the former approach (Miller & Lamb 2016).The simultaneous measurements of mass and radius using this method have been achieved by the Neutron star Interior Composition Explorer (NICER) mission for two NSs: the nearby isolated NS PSR J0030+0451 (Riley et al. 2019;Miller et al. 2019) and the millisecond pulsar PSR J0740+6620 (Riley et al. 2021;Miller et al. 2021) with the highest mass known.In compact binary coalescing systems involving at least one NS, the imprints left by tidal effects (Flanagan & Hinderer 2008;Damour et al. 2012) in the emitted gravitational waves (GWs) provide a new avenue to investigate the EOS (Abbott et al. 2018;De et al. 2018).The first-ever measurements of mass and tidal deformability for a pair of coalescing NSs were obtained from the binary neutron star (BNS) merger event GW170817 (Abbott et al. 2019b), and subsequently other measurements were from GW190425 (Abbott et al. 2020).
The merger of a BNS system can lead to multiple outcomes, depending mainly on the total (gravitational) mass of the system (M tot ) and the maximum mass of a nonrotating NS.We briefly summarize the diverse outcomes as follows (see, e.g., Baiotti & Rezzolla 2017, for a comprehensive review): (1) If the baryonic mass of the remnant (M b,rem ) is less than that of the maximum-mass nonrotating NS (M b,TOV ), it will initially form a rapidly rotating NS and eventually evolve into a stable, slowly rotating NS. (2) When the gravitational mass of the remnant (M rem ) is below the critical maximum mass (M crit ) supported by uniform rotation, a supramassive neutron star (SMNS) will form.Subsequently, it will lose angular momentum (or accrete matter from the envelope; Margalit et al. 2022) and approach the stability line, leading to the formation of a rotating black hole (BH).( 3) If M rem is slightly larger than M crit , a temporary existence of a hypermassive neutron star (HMNS) that can sustain more mass than uniform rotation becomes highly probable.After dissipative effects remove the differential rotation, it will become unstable and collapse, transforming either into an SMNS or directly into a rotating BH. (4) In cases where M tot exceeds a certain threshold for prompt collapse (M thres ), a highly spinning BH (with a dimensionless angular momentum χ ≃ 0.7 − 0.8, Kastaun et al. 2013;Bernuzzi et al. 2014) can immediately form after the merger.In the second scenario, it has been found that, at the time of collapse, the gravitational mass M rem and the dimensionless angular momentum χ rem of the remnant satisfy a tight universal relation (Shao et al. 2020b) where M crit = M rem , c 2 = 0.0902(2), c 4 = 0.0193(2), and C TOV = GM TOV /R TOV c 2 (G and c are the gravitational constant and the speed of light, respectively) represents the compactness, which is associated with the mass and radius of a nonrotating NS at the maximum mass configuration.This relation also holds for the mass (M BH ) and spin (χ BH ) of the leftover BH, due to negligible GW radiation (Dhani et al. 2023) and angular momentum conservation during the collapse.Such a low-mass second-generation (2G) BH has the potential to encounter another NS, forming a binary system and possibly being detectable by GW detectors (Gupta et al. 2020;Liu & Lai 2021;Gayathri et al. 2023).
If the mass and spin of the 2G BH can be precisely/accurately measured through the GW signal, a series of such measurements could lead to the determination of M TOV and R TOV by using the universal relation, which is the main motivation of this work.

METHODS
To illustrate our approach for constraining M TOV and R TOV , we first select BSk21 (Goriely et al. 2010;Potekhin et al. 2013) as the underlying EOS due to its good agreement with current observational constraints (see, e.g., Tang et al. 2021).Then, we carry out the calculation of the stability line (M crit − χ) for the given EOS using the Rapidly Rotating Neutron Star (RNS) code (Cook et al. 1994;Stergioulas & Friedman 1995;Stergioulas 1996), which involves the following two steps: (1) Compute the static and Keplerian sequences to determine their respective maximum (critical) masses (including, e.g., M TOV , M b,TOV , and M b,kep ), along with other associated quantities, such as the central energy densities ϵ c,TOV and ϵ c,kep , as well as the angular momentum J kep .(2) Slice the identified J kep and apply the 'fminbound' algorithm from 'scipy.optimize' with a boundary of (ϵ c,kep , ϵ c,TOV ) to precisely locate the turning point within the J-constant model.Through analyses of NS mass distribution and EOS-insensitive relations, Shao et al. (2020a) found that the prospect of forming SMNS during BNS mergers is quite promising (see also Fan et al. 2023, for the latest validation).Here we perform calculations to determine the properties of the final outcome resulting from the merger within the presently identified BNS population (see Table 1 for details), while considering the specified EOS.The brief steps are summarized in the flowchart depicted in Fig. 1.First, we calculate the threshold mass for prompt collapse (M thres ) using the fitting formula provided in Tootle et al. (2021) (see also Köppel et al. 2019 andBauswein et al. 2021).If the total gravitational mass of the binary is below M thres , we then employ the empirical fitting formula P 2 2 (q, Λ) (with recommended coefficients) from Nedora et al. (2022) to estimate the masses of both the dynamical ejecta (M eje ) and the remnant disk (M disk ).Utilizing the conservation of baryonic mass within the system, we subsequently estimate the baryonic mass of the remnant However, in reality, the measured values of M BH (associated with M crit ) and χ BH in the merger of a 2G BH with an NS are subject to uncertainties.To investigate how well the mass and spin parameters can be measured for such BH-NS systems by future GW detectors, e.g., the Einstein Telescope (ET; Punturo et al. 2010) and the Cosmic Explorer (CE; Reitze et al. 2019), we estimate the measurement uncertainties for a simulated population of these systems using a novel Fisher-matrix code GWFAST (Iacovelli et al. 2022a,b) with IMRPhenomXP NRTidalv2 waveform template (Dietrich et al. 2019;Pratten et al. 2021).We assume that χ BH is uniformly distributed in the range of (0.1, 0.7), and M BH is determined using Eq. ( 1) with the true values of M TOV and R TOV .The spin of the NS is assumed to have a low magnitude, distributed uniformly between 0 and 0.05.The NS mass follows a bimodal Gaussian distribution as defined in Eq. ( 1) of Shao et al. (2020a), and is characterized by 57, and a maximum-mass cutoff at M TOV = 2.28M ⊙ .Besides, it is assumed that both the spin orientations of the BH and the NS are isotropically distributed, and the azimuthal position of orbital angular momentum on its cone of precession about the total angular momentum ϕ JL and the azimuthal angle separating the spin vectors ϕ 12 are uniform in (0,2π).We employ the redshift distribution similar to that of Iacovelli et al. (2022a), and with a longer minimum time delay t d,min = 100 Myr.Additionally, we maintain the assumption of no redshift evolution in the mass distribution.The luminosity distance is computed from the redshift assuming Planck18 flat ΛCDM (Planck Collaboration et al. 2020).We then use the complete GW parameters (for a list and description of these parameters, one can refer, for example, to Table 1 of Tang et al. (2020)) to calculate the Fisher Information Matrix (FIM), except for fixing the tidal deformability of the BH to 0. We discard all events where the inversion error of the FIM exceeds 0.05 and where the optimal signal-to-noise ratio (S/N) is below 12.After obtaining the distributions of ∆M BH (source frame mass) and ∆χ BH , we then randomly assign measurement errors to a set of mock events.For each event, we introduce random shifts drawn from normal distributions with standard deviations corresponding to the measurement errors.These shifts are added to the actual values of M BH and χ BH , yielding the corresponding 'measured' mean values.This process is intended to simulate the potential effects encountered in real data analysis, where slight biases in parameter inference may arise.Finally, we utilize the orthogonal distance regression (ODR) method implemented in the 'scipy' library to perform data fitting.This is specifically applied to the mock data M BH and χ BH , taking into consideration their respective measurement errors.This approach allows us to determine the best-fit parameters for M TOV and R TOV , along with their associated uncertainties.

RESULTS
As illustrated in Fig. 2, based on the BSk21 EOS, the majority of Galactic double NS systems are expected to form SMNS if they were to merge.The heaviest BNS system currently observed, GW190425, is associated with a high-spin value of the remnant, as indicated by numerical-relativity simulations (Dudi et al. 2022).In contrast, some lighter double NS systems are presumed to produce low-spin stable NSs.The universal relation, characterized by the true M TOV and R TOV parameters of BSk21 (represented by the red curve), fits the 'observed' remnant masses and spins well.By varying these parameters, we notice that M TOV has a strong influence on the overall data fit, while R TOV predominantly impacts the fitting of high-spin data.

M(M )
Figure 2. Stability line (Mcrit − χ, depicted as a gray dotted line) of SMNS for the BSk21 EOS, along with the anticipated remnant masses and spins for BNS systems.These remnant characteristics can be divided into three categories, symbolized by blue, orange, and black markers.These markers correspond to SMNS remnants, stable NSs, and promptly formed BHs, respectively.The cyan, red, pink, and magenta lines illustrate curves based on universal relation, each with varying MTOV and RTOV (corresponding to values shown in the parentheses of the legend).The vertical black dotted lines mark the χ = 0.1 and χ = 0.7 boundaries, which roughly separate the three categories.
For (2G)BH-NS systems that are detectable (with S/N≥ 12) by future GW detectors, we showcase the distributions of measurement uncertainties for the source frame mass, M BH , and spin magnitude, χ BH , of the BH in Fig. 3.The detector configuration is consistent with the 'ET+2CE' settings as used in Iacovelli et al. (2022a).We determine the measurement uncertainty of the source frame mass of BHs by where ∆M z BH is the measurement uncertainty of the detector frame mass of BHs.The redshift uncertainty ∆z is calculated through where H(z) is the Hubble parameter as a function of redshift.As illustrated in the top panels of Fig. 3, it is evident that the proposed third-generation GW detectors exhibit significant capability in constraining the component masses.
For BH-NS mergers, there is the possibility of producing electromagnetic (EM) radiations, such as gamma-ray bursts 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 (GRBs) and kilonovae.This can help determine the redshift (we assume ∆z ∼ 0 in this case), consequently reducing the measurement uncertainty of M BH .The ∆M BH derived solely from GW observations (hereafter referred to as 'GW only') is approximately 0.093 +0.060 −0.049 M ⊙ .This capability is further enhanced with the assistance of EM counterparts (hereafter referred to as 'GW+EM'), with ∆M BH = 0.026 +0.026 −0.015 M ⊙ .While the ability to constrain the spin magnitude is not as pronounced as for the mass, as shown in the bottom right panel of Fig. 3.
With the distribution of measurement uncertainties in hand, we can proceed to generate mock measurements and employ ODR fitting to determine the values of M TOV and R TOV .To illustrate this method, we initially utilize a set of 20 'golden' events, each with precise measurements of M BH and χ BH (specifically, both have measurement error values ≲ 0.05).As shown in the top left panel of Fig. 4, it is feasible to determine the unknown parameters (i.e., M TOV and R TOV ), and their true values can be accurately recovered.In principle, with precise measurements of just two samples of M BH and χ BH , we can solve for these two parameters.Increasing the number of data points will effectively reduce the uncertainties of M TOV and R TOV .To simulate a realistic scenario, we generate samples based on the capabilities of third-generation GW detectors for constraining mass and spin parameters.As shown in the top right panel of Fig. 4, we notice that 100 'GW only' detections can produce a precision of 0.014 M ⊙ for M TOV and 0.9 km for R TOV .Furthermore, when we include additional EM counterpart information (i.e., the 'GW+EM' case), the uncertainties further decrease to 0.01 M ⊙ for M TOV and 0.6 km for R TOV .However, if there are potential contaminants from additional BH populations (e.g., low-mass 1G BHs with high spins produced by accretion), M TOV and R TOV may be biased as shown in the bottom right panel of Fig. 4. Considering that there are fluctuations in the random process used to generate the mock events, we then simulate 5000 instances of event samples.The comprehensive results are illustrated in Fig. 5, showcasing the distributions of recovered best-fit values along with their corresponding uncertainties from the ODR fitting.These results demonstrate that as the number of detected events increases, both M TOV and R TOV can be measured with good accuracy and precision.In the optimistic scenario, the precision of determining M TOV and R TOV can reach approximately 0.01 M ⊙ and 0.6 km, respectively.1).The top left panel serves as an illustration of the method, showcasing the fitting results with 20 golden events.The top right and bottom left panels exhibit the fitting results for 100 samples randomly generated from the distribution of ∆MBH for the 'GW only' and 'GW+EM' cases, respectively.The bottom right panel evaluates the effect of outlier contamination on the fitting procedure.In each panel, the red solid curve represents the universal relation with the best-fit parameter values, and the black dashed line represents that with the true values.The light blue band indicates the 2σ interval, and the gray error bars denote the mock measurements.The numbers shown above each panel represent the best-fit values and their corresponding 1σ uncertainties.

SUMMARY AND DISCUSSION
In this work, we present a method for determining the mass and radius of a maximum-mass nonrotating neutron star, which exploits the tight universal relation between the critical mass of a uniformly rotating SMNS and its dimensionless angular momentum at the point of collapse.Such a universal relation is associated with the mass (M TOV ) and radius (R TOV ) of the maximum-mass nonrotating neutron star.The characteristics of the BHs formed after the collapse of SMNSs also largely conform to this universal relation.Assuming that the remnant BH has the opportunity to merge with another NS, and that this event is detectable by GW detectors, it is theoretically possible to determine the parameters within the universal relation by utilizing measurements of the mass and spin of the remnant BH.To evaluate the feasibility of this approach, we initially assess the capability of the proposed third-generation GW detectors to measure the mass and spin parameters.This is done using the Fisher-matrix approach, as implemented in the GWFAST code.Our findings suggest that it is feasible to accurately measure the source frame BH mass, while the measurement of the spin parameter is satisfactory, though not excellent.Following these findings, we use ODR fitting to analyze mock measurement samples, considering different numbers of detected events and the inclusion or exclusion of EM counterpart information.We find that it is highly probable to accurately and precisely measure both the M TOV and R TOV .In scenarios with 100 detections that include both GW and EM counterpart information, the M TOV and R TOV can be measured with a precision of approximately 0.01 M ⊙ and 0.6 km, respectively.This precise determination of the M TOV and R TOV holds significant potential for constraining the unknown EOS of dense matter and for enhancing our understanding of high-energy phenomena in astrophysics, including predicting the outcome of the BNS merger and determining the central engine of GRBs.Such independent measurements can also directly test the M TOV and R TOV inferred in other ways (see Fan et al. 2020, 2023, andreferences therein).Additionally, since the M crit − χ universal relation differs between NSs and quark stars (Bozzola et al. 2018;Zhou et al. 2019), our method could potentially be further extended to differentiate the nature of compact stars.Finally, we would like to discuss some potential difficulties/challenges when people apply our method to real observations.Firstly, distinguishing between high-mass NSs and low-mass BHs based solely on GW signals is challenging.However, there are promising methods proposed for third-generation GW detectors to mitigate this issue (Chen et al. 2020;Brown et al. 2022).Moreover, observation data from EM counterparts, such as kilonova observations, can also assist in distinguishing between the two types (Kawaguchi et al. 2020) and in constraining mass and spin parameters (Barbieri et al. 2019) when available.Secondly, confirming the presence of 2G BHs presents its own unique challenges.Current mass-spin population studies for BBHs suggest that first-generation (1G) BHs typically exhibit low spin (Abbott et al. 2023;Wang et al. 2022;Li et al. 2023).And binary evolution simulation also shows that the BHs in ordinary BH-NS systems have dimensionless spin parameter less than about 0.2 (Xing et al. 2023).If we attribute high-spin values of low-mass BHs mainly to BNS mergers, our sample remains relatively uncontaminated.However, if we consider accretion onto low-mass 1G BHs as a contributing factor, then these additional BH populations could potentially contaminate our sample and cause biased measurements of M TOV and R TOV as discussed in the Results section.In such instances, these populations should be carefully modeled and subsequently excluded.Thirdly, the detection rate of (2G)BH-NS mergers carries a significant degree of uncertainty.Nevertheless, some studies have indicated that certain detected GW events support the hypothesis of a (2G)BH-NS merger origin (Gupta et al. 2020;Liu & Lai 2021;Gayathri et al. 2023), and the estimated merger rate appears promising (Kritos & Cholis 2021).Fourth, the 2G BHs may not remain in the same state as they were initially formed from the collapse of SMNSs.This is due to the fact that BNS mergers might produce GRBs, and the Blandford-Znajek mechanism (Blandford & Znajek 1977) could cause the newly formed BHs to spin down.In addition, the mass that falls into these BHs may contribute to their growth in mass.As indicated by the results from numerical simulations (Fujibayashi et al. 2020), the shift in χ BH appears insignificant compared to its measurement uncertainty.However, the mass that falls into the BH is a critical factor, and the masses of 2G BHs (M BH ) need to be accurately calibrated.
We would like to thank Francesco Iacovelli for the kind help in using GWFAST.This work is supported in part by the Project for Special Research Assistant and the Project for Young Scientists in Basic Research (No. YSBR-088) of the Chinese Academy of Sciences, the National Natural Science Foundation of China under grant Nos. 12233011, 11921003, 11933010, 12073080, and 12303056, and by the General Fund (No. 2023M733736) of the China Postdoctoral Science Foundation.

Figure 1 .
Figure 1.Flowchart for determining the properties of the final outcome resulting from the merger of a BNS system.

Figure 3 .
Figure 3. Distributions of measurement uncertainties for the redshift z and the source frame mass MBH and spin magnitude χBH of a BH under the 'ET+2CE' detector configuration.Top panels present the ∆MBH results with (left) and without (right) electromagnetic counterparts information.Bottom panels show the results of ∆z (left) derived from the uncertainty of luminosity distance and the results of ∆χBH (right).The smooth, black solid lines represent the kernel density estimates (KDEs) of the distributions.The numbers reported above each panel are the median and 68.3% symmetric credible interval.

Figure 4 .
Figure4.Results of ODR fitting on randomly generated measurement samples using the universal relation described by Eq. (1).The top left panel serves as an illustration of the method, showcasing the fitting results with 20 golden events.The top right and bottom left panels exhibit the fitting results for 100 samples randomly generated from the distribution of ∆MBH for the 'GW only' and 'GW+EM' cases, respectively.The bottom right panel evaluates the effect of outlier contamination on the fitting procedure.In each panel, the red solid curve represents the universal relation with the best-fit parameter values, and the black dashed line represents that with the true values.The light blue band indicates the 2σ interval, and the gray error bars denote the mock measurements.The numbers shown above each panel represent the best-fit values and their corresponding 1σ uncertainties.

Figure 5 .
Figure 5. Distributions of the best-fit values and corresponding 1σ uncertainties of MTOV and RTOV.The blue and green violin plots represent the 'GW only' and 'GW+EM' results, respectively.The orange dashed/dotted lines denote the median and the 68.3% interval.The red star markers represent the true values.

Table 1 .
Mass Measurements of BNS Systems and Predicted Remnant's Properties (based on EOS BSk21) All the masses are measured in units of solar mass.The M b,TOV , M b,kep , and χ kep for EOS BSk21 are 2.72M⊙, 3.22M⊙, and 0.7, respectively.