Radius and equation of state constraints from massive neutron stars and GW190814
Yeunhwan Lim, Anirban Bhattacharya, Jeremy W. Holt, Debdeep Pati
RRevisiting constraints on the maximum neutron star mass in light of GW190814
Yeunhwan Lim,
1, 2, 3, ∗ Anirban Bhattacharya, † Jeremy W. Holt,
5, 6, ‡ and Debdeep Pati § Max-Planck-Institut f¨ur Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany Institut f¨ur Kernphysik, Technische Universit¨at Darmstadt, 64289 Darmstadt, Germany ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum f¨ur Schwerionenforschung GmbH, 64291 Darmstadt, Germany Department of Statistics, Texas A&M University, College Station, TX 77843, USA Cyclotron Institute, Texas A&M University, College Station, TX 77843, USA Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA (Dated: July 16, 2020)We investigate the maximum neutron star mass based on constraints from low-energy nuclearphysics, neutron star tidal deformabilities from GW170817, and simultaneous mass-radius mea-surements of PSR J0030+045 from NICER. Our prior distribution is based on a combination ofnuclear modeling valid in the vicinity of normal nuclear densities together with the assumption ofa maximally stiff equation of state at high densities. The transition density is treated as a modelparameter with uniform prior. Bayesian likelihood functions involving measured neutron star tidaldeformabilities and radii are subsequently used to generate equation of state posteriors. We demon-strate that a modification of the highly uncertain supra-saturation density equation of state allowsfor the support of 2 . − . M (cid:12) neutron stars without strongly modifying the properties (radius,tidal deformability, and moment of inertia) of ∼ . M (cid:12) neutron stars. In our analysis, only thesoftest equations of state are eliminated under this scenario. However, the properties of neutronstars with masses ∼ . M (cid:12) are significantly different under the two competing assumptions thatthe GW190814 secondary was a black hole or a neutron star. PACS numbers: 21.30.-x, 21.65.Ef,
Introduction—
Recently, the LIGO/Virgo Collabora-tion (LVC) has reported measurements [1] of gravita-tional waves resulting from a 2 . − . M (cid:12) “mass-gap”object [2] in binary coalescence with a heavy (22 . − . M (cid:12) ) companion black hole. Not only are the massratio of q = 0 . +0 . − . and inferred merger rate of1 −
23 Gpc − yr − challenging to describe [1, 3, 4] withtraditional binary evolutionary models, but taken at facevalue, the mass-gap secondary component in the observa-tion represents the discovery of either the heaviest knownneutron star (NS) or the lightest known black hole (BH),though see Ref. [5] for an alternative scenario in whichthe source of GW190814 is conjectured to be a normalNSBH merger amplified via gravitational lensing. Nei-ther the absence of a measurable tidal deformation sig-nature in the gravitational waveform nor the absence ofan electromagnetic counterpart would be unexpected [6]for a NSBH merger at the extreme mass ratio reportedin GW190814. However, equation of state inferences [7]based on GW170817 and properties of its electromagneticcounterpart [8–12] suggest that such heavy neutron starswould be challenging to describe with traditional neutronstar equations of state founded in nuclear physics modelswell constrained up to one or two times normal nucleardensities.Given the highly uncertain nature of matter at densi-ties exceeding two to three times normal nuclear matter ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] density ( n = 0 .
16 fm − = 2 . × g/cm ), where thereexists no theoretical description of the strong interactionwith controlled uncertainties, in this work we explore theextreme scenario in which the high-density equation ofstate is maximally stiff and therefore can support theheaviest neutron stars. Our low-density equation of stateis constrained by nuclear theory and experiment as wellas recent radius and tidal deformability measurements of ∼ . M (cid:12) neutron stars, while the transition region to themaximally stiff equation of state is varied between 2 − n .We explore the minimum transition density required tosupport 2 . − . M (cid:12) neutron stars and find that it lies inthe region n ∼ . n , which is below the central densityof neutron stars with masses M ∼ . M (cid:12) . Nevertheless,we find that the existence of massive 2 . − . M (cid:12) neutronstars does not strongly constrain the bulk properties oftypical lighter neutron stars, and only the softest equa-tions of state with small radii and tidal deformabilitiesare excluded. In contrast, the radii and tidal deformabil-ities of heavy neutron stars with M ∼ . M (cid:12) differ moresignificantly under the two competing scenarios that thesecondary component of GW190814 is a black hole or aneutron star. Bayesian modeling of the neutron star equation ofstate—
Experimentally measured nuclear binding ener-gies and bulk oscillation modes constrain [13, 14] the nu-clear equation of state around normal nuclear density n for matter consisting of nearly equal numbers of neu-trons and protons. Neutron-rich matter, on the otherhand, is challenging to produce and study in the labo-ratory, and therefore the principal nuclear physics con-straints on the neutron star equation of state rely in oneway or another on nuclear theory models, which nowa- a r X i v : . [ nu c l - t h ] J u l days have a firm foundation in chiral effective field the-ory [15–17], the low-energy realization of quantum chro-modynamics. Previously, we have constructed [18, 19]Bayesian posterior probability distributions for the neu-tron star equation of state that incorporate constraintsfrom chiral effective field theory [20, 21] and experiment[22, 23]. When these models were extrapolated to highdensities, the maximum neutron star mass was found tobe M (cid:39) . M (cid:12) . Numerous other works have employedchiral effective field theory to study the neutron matterequation of state [24–32], neutron star radii [33–35], tidaldeformabilities [36–38], and moments of inertia [39, 40].In describing the properties of the heaviest neutron stars,whose central densities can reach up to n = 5 − n , all ofthese models perform extrapolations into regions wherethe composition and dynamics are poorly understood.We take as a model for the low-density equation ofstate a Taylor series expansion in the Fermi momentum k F ∼ n / . This is justified since the bulk nuclear matterequation of state is parametrized by baryon number den-sities with polytropic index [41, 42] even in the presenceof quark matter or phase transitions. The energy den-sity functional for neutron star matter is built from boththe pure neutron matter and isosopin-symmetric nuclearmatter equations of state, interpolated to beta equilib-rium conditions ( µ n = µ p + µ e ) enforcing a well-justified[43] quadratic dependence on the proton-neutron asym-metry parameter δ = ( n n − n p ) / ( n n + n p ): E ( n, δ ) = 12 m τ n + 12 m τ p + [1 − δ ] f s ( n ) + δ f n ( n ) , (1)where τ n ( τ p ) is the neutron (proton) kinetic energy den-sities and f s ( f n ) refers to the isospin-symmetric nuclearmatter (pure neutron matter) potential energy densityexpanded as follows: f s ( n ) = (cid:88) i =0 a i n (2+ i/ , f n ( n ) = (cid:88) i =0 b i n (2+ i/ . (2)In Eq. (2) the isospin-symmetric nuclear matter co-efficients (cid:126)a = { a , a , a , a } are obtained by fitting to10 equation of state calculations in chiral effective fieldtheory [20] up to the density 2 n . We have shown inprevious works [18] that lowering the maximum fittingdensity to 1 . n does not qualitatively modify our priordistributions. We then implemented experimental likeli-hood functions involving the { a , a , a , a } parametersfrom empirical nuclear matter properties, such as the sat-uration energy, saturation density, incompressibility, andskewness averaged over 205 realistic mean field modelsfitted to the binding energies and bulk properties of fi-nite nuclei [13]. For the parameters (cid:126)b = { b , b , b , b } entering in the pure neutron matter energy density func-tional f n ( n ), we first fit to a set of 10 chiral effective fieldtheory neutron matter calculations [20] up to the den-sity 2 n . The resulting multivariate distribution is thenrefined by imposing nuclear experimental constraints onthe isospin-asymmetry energy at saturation density and its higher-order derivatives in the density [21, 44]. Inall of our neutron star structure models, we construct arealistic outer and inner crust using the same parame-ters ( (cid:126)a,(cid:126)b ) in a unified way implementing the liquid dropmodel as explained in more detail in Ref. [45].To explore the widest range of maximum neutron starmasses, we extend this previous model for the equationof state probability distribution to include a transitionto the maximally-stiff equation of state consistent withrelativity, defined when the speed of sound is equal tothe speed of light. The transition density n t is takento have a uniform prior in the range 2 n < n t < n .A critical density beyond 4 n is of course possible, butwe find that it gives no significant modification to theequation of state prior. Formally, we employ a second-order phase transition where the phase transition startsat E = E and ends at E = E . Beyond E , the speed ofsound is assumed to be equal to the speed of light, andthus the pressure and energy density have a linear rela-tion. Between E and E , the speed of sound is assumedto increase linearly as a function of energy density: c s ( E ) = c + (1 − c ) E − E ∆ E , ∆ E = E − E . (3)The pressure between the phase transition density is thenobtained from the integration of the speed of sound: P = P + c ( E − E ) + (1 − c )2∆ E ( E − E ) , (4)where ∆ E = E .The approach described above defines the prior dis-tribution π ( · ) associated with the neutron star equationof state parameters θ = ( (cid:126)a ; (cid:126)b ) T . We construct Bayesianposterior probability distributions as follows. Having ob-served neutron star tidal deformabilities associated withGW170817 [7, 46–48] and simultaneous mass-radius mea-surements [49, 50] of PSR J0030+045 from the NICERmission, the posterior distribution of θ is proportional to L ( θ ) π ( θ ), where L ( θ ) = (cid:26) (cid:89) i =1 L MRi ( θ ) (cid:27) (cid:26) (cid:89) i =1 L M Λ i ( θ ) (cid:27) , (5)is the likelihood function of θ . In Eq. (5), L MRi ( θ ) for i = 1 , L M Λ i ( θ ) for i = 1 , L MR likelihood(the L M Λ terms are similarly derived and we omit thedetails here). We first introduce some notation to con-nect the parameter θ to the NICER mass-radius measure-ments. Let R θ ( M ) denote the (unique) radius-versus-mass curve corresponding to the set of parameters in θ .Each R θ ( · ) curve has its own maximum mass M max θ abovewhich the neutron star would collapse to a black hole, andhence the domain of R θ ( · ) is ( M min θ , M max θ ), where in allcases we take M min θ = 1 . M (cid:12) . The exact choice of M min θ is not particularly crucial, since neither NICER nor LIGOmeasurements have significant statistical weight around1 M (cid:12) . The main purpose of this additional notation isto provide a prescription to randomly generate a ( M, R )pair a priori , which proceeds by (i) sampling θ ∼ π , (ii)given θ , sampling M uniformly between ( M min θ , M max θ ),and (iii) setting R = R θ ( M ). The uniform generationof M is justified since the equation of state is agnosticabout the location on an R ( M ) curve.Although the correlated uncertainty corresponding toeither NICER measurement resembles a tilted ellipse, acloser inspection of the contour plots reveal departuresfrom normality. As a result, we refrain from using a para-metric Gaussian likelihood function, and instead builda non-parametric likelihood using a kernel density esti-mator (kde). Specifically, we separately fit kernel den-sity estimators (cid:98) f and (cid:98) f to the ( M, R ) posterior sam-ples corresponding to Fig. 7b of Ref. [49] and Fig. 19(“ST+PST”) of Ref. [50]. We used the R package ks to fitthe kde, employing a Gaussian kernel and the bivariatesmoothed cross-validation estimator for the bandwidthmatrix. Then, we consider an “average” of these fitteddensities over an R ( M ) curve as the corresponding like-lihood, i.e., for i = 1 , L MRi ( θ ) = (cid:90) M max θ M min θ (cid:98) f i (cid:0) M, R θ ( M ) (cid:1) dMM max θ − M min θ . (6)Next, we describe how we incorporate the secondary“mass-gap” object into our likelihood function. That itsdistribution is constrained in the interval 2 . − . M (cid:12) makes it a candidate for either the lightest BH or theheaviest NS ever observed. From Fig. 4 of [7], the distri-bution of the secondary mass M s resembles a Gaussiandistribution N( µ s , σ ) with mean and 90% intervals givenby 2 . +0 . − . . This leads to µ s = 2 .
59 and σ s = 0 . θ from the equation of state,the secondary object is realizable as a neutron star withprobability given by P ( M s ≤ M max θ ) = Φ (cid:18) M max θ − µ s σ s (cid:19) , (7)where Φ( · ) denotes the standard Gaussian cumulativedistribution function, Φ( x ) = (cid:82) x −∞ (2 π ) − / e − x / dx .Eq. (7) then defines the likelihood of the object assum-ing it to be a neutron star, denoted L NSs ( θ ), which whenmultiplied with L ( θ ) defined in Eq. (5) gives the overalllikelihood for θ . Similarly, if we assume the object is ablack hole, then the likelihood involves the probabilitygiven by L BHs ( θ ) : = 1 − Φ (cid:8) ( M max θ − µ s ) /σ s } . Results—
In Fig. 1 we show the mass and radius prob-ability distributions based on the Bayesian analysis de-scribed above. In all subpanels of Fig. 1 the green andblue contours represent the 68% (solid lines) and 95%(dashed lines) credibility bands obtained from our kernel . . . . . . . M ( M (cid:12) ) Smooth high-density prior MillerRiley Modified high-density prior R (km) . . . . . . . M ( M (cid:12) ) Posterior : M max > . M (cid:12) R (km) Posterior : M max < . M (cid:12) − − − FIG. 1. (Color online) Mass and radius probability dis-tributions for the (top-left) prior without high-density ex-trapolation, (top-right) prior with high-density extrapolation,(bottom-left) posterior supporting ∼ . M (cid:12) neutron stars,and (bottom-right) posterior not supporting ∼ . M (cid:12) neu-tron stars. The green [50] and blue [49] contours representthe NICER 68% (solid) and 95% (dashed) credibility bands. density estimators associated with the Riley et al. [50]and Miller et al. [49] analyses of NICER x-ray waveformdata from PSR J0030+045. The top-left figure is our pre-vious prior [18] without a high-density extrapolation, thetop-right panel is our new prior with uniformly varyingtransition density 2 n < n t < n . In order to sup-port ∼ . M (cid:12) neutron stars, we find that the transi-tion density must satisfy n t < . n , indicating that therelatively soft neutron star equations of state predictedby chiral effective field theory must become fairly stiffsoon after their natural breakdown scale in the range1 − n . We see that the inclusion of the maximallystiff equation of state at high densities naturally leadsto much larger maximum neutron star masses, up to M max = 2 . M (cid:12) for the lowest value of the transitiondensity considered n t = 2 n . We note that this newmaximum neutron star mass of M max = 2 . M (cid:12) is al-most certainly unphysical since it lies above the totalmass M tot (cid:39) . M (cid:12) of the GW170817 remnant, whichis expected [51] to have collapsed to a black hole after be-ing supported initially through differential rotation. Thebottom-left and bottom-right panels of Fig. 1 representthe posterior mass-radius probability distributions underthe assumption that the secondary in GW190814 wasa neutron star or a black hole, respectively. Interest-ingly, we see that for typical neutron stars with masses M ∼ . M (cid:12) , the distribution of radii is not stronglydifferent under the two interpretations of the GW190814 R (km) . . . P D F , P ( R ) P ( R | M max > . M (cid:12) ) P ( R | M max < . M (cid:12) ) M = 2 . M (cid:12) FIG. 2. (Color online) Radius distribution for a heavy neutronstar with mass 2 . M (cid:12) (shaded) under the two assumptionsthat the GW190814 secondary was a black hole (blue) or aneutron star (red). The dashed lines correspond to varyingneutron star masses in the range 2 . − . M (cid:12) with spacing∆ M = 0 . M (cid:12) . secondary. This is due to the fact that the bulk proper-ties of the average neutron star are strongly correlated[18, 52, 53] with the pressure of beta-equilibrium matterat the density n = 2 n , which is close to the regime wherenuclear physics places strong constraints on the equationof state. However, we do observe that the existence ofmassive (2 . − . M (cid:12) ) neutron stars would rule out thesoftest equations of state.Our previous finding [18] for the radius of a 1 . M (cid:12) neutron star at the 95% credibility level was 10 . ≤ R . ≤ . . . ≤ R . ≤ . L BHs ( θ ) and 11 . ≤ R . ≤ . L NSs ( θ ). Our finding is consistent with the determi-nation of mass and radius from the cooling tail method[54], where the source for the analysis is different. Wesee from Fig. 1 that heavy neutron stars, such as PSRJ0740+6620 with mass 2 . +0 . − . M (cid:12) [55], have signifi-cantly different radius probability distributions under ourtwo assumptions for the GW190814 likelihood, L BHs ( θ )and L NSs ( θ ). In Fig. 2 we show in the shaded regionsthe posterior probability distributions for the radius of a2 . M (cid:12) neutron star under the two assumptions that theGW190814 secondary was a black hole (blue) or a neu-tron star (red). The notation P ( R | M max > . M (cid:12) ) cor-responds to the likelihood assumption L NSs ( θ ) and like-wise P ( R | M max < . M (cid:12) ) corresponds to L BHs ( θ ). Thestiff equations of state needed to support the heaviestneutron stars produce a narrow and large neutron starradius at this mass, while softer equations of state lead tostatistically significant smaller radii. In Fig. 2 the dashedlines correspond to different heavy neutron star massesranging from 2 . − . M (cid:12) at spacing ∆ M = 0 . M (cid:12) .The radius distributions for the lightest neutron stars ex- . . . . M ( M (cid:12) ) Λ P (Λ | M max > . M (cid:12) ) . . . . M ( M (cid:12) ) P (Λ | M max < . M (cid:12) )10 − − − FIG. 3. (Color online) Probability distributions for the tidaldeformability versus mass under the two assumptions L NSs ( θ )(left) and L BHs ( θ ) (right). The blue and green contours showthe 68% (solid) and 95% (dashed) credibility bands associatedwith the primary and secondary in GW170817, respectively. tend to the smallest radii for both posteriors. We find atthe 95% credibility level that the radius of a 2 . M (cid:12) neutron star is 10 . < R . < . L BHs ( θ ) and 11 . < R . < . L NSs ( θ ).In Fig. 3 we show the two posterior probability distri-butions for the tidal deformability as a function of massunder the two assumption for the likelihood L BHs ( θ ) and L NSs ( θ ). In both subpanels the blue and green contoursdenote the 68% (solid lines) and 95% (dashed lines) cred-ibility bands from our kde associated with the primaryand secondary components, respectively, of GW170817.In Fig. 4 we show the tidal deformability of a typi-cal 1 . M (cid:12) neutron star under the assumptions that theGW190814 secondary was a neutron star (red) or blackhole (blue). Our previous 95% credibility interval in Ref.[18] was found to be 136 < Λ . < L BHs ( θ ), we find170 < Λ . < L NSs ( θ ),we likewise find 313 < Λ . <
575 at the 95% credibilitylevel. There remains a significant overlap between thetwo distributions, but we observe a broad low tidal de-formability region possible only in the absence of heavyneutron stars with masses 2 . − . M (cid:12) .Following the discovery of the double pulsar systemJ0737-3039, it was suggested [56, 57] that precise radiotiming measurements could enable the extraction of spin-orbit coupling effects on the system’s periastron advanceand hence the moment of inertia of PSR J0737-3039A.The mass of PSR J0737-3039A is precisely known to be1 . M (cid:12) , and in Fig. 5 we plot the associated predictionsfor its moment of inertia assuming that the equation ofstate can support 2 . − . M (cid:12) neutron stars (red) or not(blue). For such a relatively light neutron star, there isan even smaller difference between the moment of iner-tia probability distributions under the two assumptions Λ . . . . . P D F , P ( Λ . ) P (Λ | M max > . M (cid:12) ) P (Λ | M max < . M (cid:12) ) M = 1 . M (cid:12) FIG. 4. (Color online) Probability distributions for the tidaldeformability of a 1 . M (cid:12) neutron star under the assumptionthat the GW190814 secondary was a black hole (blue) or aneutron star (red). . . . . . . I (10 g cm ) . . . . P D F , P ( I ) P ( I | M max > . M (cid:12) ) P ( I | M max < . M (cid:12) ) M = 1 . M (cid:12) FIG. 5. (Color online) Probability distributions for the mo-ment of inertia of PSR J0737-3039A with mass 1 . M (cid:12) un-der the two assumptions that the GW190814 secondary wasa black hole (blue) or a neutron star (red). L BHs ( θ ) and L NSs ( θ ). Previously, we found [39] that atthe 95% credibility level the moment of inertia of J0737-3039A should lie in the range 1 . × < I . < . × g cm . In our revised modeling, includingNICER and GW170817 data, we now find a new sta-tistical range 1 . × < I . < . × g cm under the assumption L BHs ( θ ) and 1 . × < I . < . × g cm under the assumption L NSs ( θ ). Afteraccounting for the NICER likelihood functions, we pre-dict a somewhat larger moment of inertia for 1 . M (cid:12) as well as a reduced statistical uncertainty. Summary—
The existence of heavy neutron stars withmasses 2 . − . M (cid:12) are a challenge to explain with equa-tions of state smoothly extrapolated from the low-densityregime (1 − n ) constrained by nuclear physics to thehighest density regime (5 − n ) encountered in neutronstar cores. We have demonstrated that a modificationof the highly uncertain supra-saturation density equa-tion of state allows for the support of 2 . − . M (cid:12) neu-tron stars consistent with state-of-the-art nuclear theorymodeling within the framework of chiral effective fieldtheory, nuclear experiments involving medium-mass andheavy isotopes, as well as current observations of neutronstar radii and tidal deformabilities, all integrated withina consistent Bayesian statistical framework. While thenature of the secondary in GW190814 cannot be deter-mined within our present modeling (see also Refs. [58–61]), we note that we have observed strong correlationsbetween the maximum neutron star mass and the radiiof heavy neutron stars. We suggest that measurementsof very massive ( ∼ . M (cid:12) ) neutron star radii (or tidaldeformabilities), such as a NICER measurement of thePSR J0740+6620 radius, may provide a useful and strongconstraint on the nuclear equation of state at supra-saturation density. ACKNOWLEDGMENTS
Y. Lim was supported by the Max Planck Societyand the Deutsche Forschungsgemeinschaft (DFG, Ger-man Research Foundation) – Project ID 279384907 –SFB 1245. Dr. Pati and Dr. Bhattacharya acknowledgesupport from NSF DMS (1854731, 1916371) and NSFCCF 1934904 (HDR-TRIPODS). In addition, Dr. Bhat-tacharya acknowledges NSF CAREER 1653404 for sup-porting this project. The work of J. W. Holt is sup-ported by the National Science Foundation under GrantNo. PHY1652199 and by the U. S. Department of EnergyNational Nuclear Security Administration under GrantNo. de-na0003841. Portions of this research were con-ducted with the advanced computing resources providedby Texas A&M High Performance Research Computing. [1] B. P. Abbott et al. (LIGO Scientific Collaborationand Virgo Collaboration), Astrophys. J. Lett. , L44(2020).[2] F. Ozel, D. Psaltis, R. Narayan, and J. E. McClintock,Astrophys. J. , 1918 (2010).[3] M. Zevin, M. Spera, C. P. Berry, and V. Kalogera,arXiv:2006.14573 (2020). [4] K. Vattis, I. S. Goldstein, and S. M. Koushiappas,arXiv:2006.15675 (2020).[5] T. Broadhurst, J. M. Diego, and G. F. Smoot,arXiv:2006.13219 (2020).[6] F. Foucart, L. Buchman, M. D. Duez, M. Grudich, L. E.Kidder, I. MacDonald, A. Mroue, H. P. Pfeiffer, M. A.Scheel, and B. Szilagyi, Phys. Rev. D , 064017 (2013). [7] B. P. Abbott et al. (The LIGO Scientific Collabora-tion and the Virgo Collaboration), Phys. Rev. Lett. ,161101 (2018).[8] A. Bauswein, O. Just, H.-T. Janka, and N. Stergioulas,Astrophys. J. Lett. , L34 (2017).[9] B. Margalit and B. D. Metzger, Astrophys. J. Lett. ,L19 (2017).[10] M. Ruiz, S. L. Shapiro, and A. Tsokaros, Phys. Rev. D , 021501 (2018).[11] L. Rezzolla, E. R. Most, and L. R. Weih, Astrophys. J.Lett. , L25 (2018).[12] D. Radice, A. Perego, F. Zappa, and S. Bernuzzi, Astro-phys. J. Lett. , L29 (2018).[13] M. Dutra, O. Louren¸co, J. S. S´a Martins, A. Delfino,J. R. Stone, and P. D. Stevenson, Phys. Rev. C ,035201 (2012).[14] M. Dutra, O. Louren¸co, S. S. Avancini, B. V. Carlson,A. Delfino, D. P. Menezes, C. Providˆencia, S. Typel, andJ. R. Stone, Phys. Rev. C , 055203 (2014).[15] S. Weinberg, Physica A , 327 (1979).[16] E. Epelbaum, H.-W. Hammer, and U.-G. Meißner, Rev.Mod. Phys. , 1773 (2009).[17] R. Machleidt and D. R. Entem, Phys. Rept. , 1(2011).[18] Y. Lim and J. W. Holt, Phys. Rev. Lett. , 062701(2018).[19] Y. Lim and J. W. Holt, Eur. Phys. J. A , 209 (2019).[20] J. W. Holt and N. Kaiser, Phys. Rev. C , 034326(2017).[21] J. W. Holt and Y. Lim, Phys. Lett. B , 77 (2018).[22] M. Dutra, O. Lourenco, J. S. S´a Martins, A. Delfino,J. R. Stone, and P. D. Stevenson, Phys. Rev. C ,035201 (2012).[23] I. Tews, J. M. Lattimer, A. Ohnishi, and E. E. Kolomeit-sev, Astrophys. J. , 105 (2017).[24] K. Hebeler and A. Schwenk, Phys. Rev. C , 014314(2010).[25] L. Coraggio, J. W. Holt, N. Itaco, R. Machleidt, andF. Sammarruca, Phys. Rev. C , 014322 (2013).[26] A. Gezerlis, I. Tews, E. Epelbaum, S. Gandolfi,K. Hebeler, A. Nogga, and A. Schwenk, Phys. Rev. Lett. , 032501 (2013).[27] A. Carbone, A. Rios, and A. Polls, Physical Review C , 054322 (2014).[28] C. Drischler, A. Carbone, K. Hebeler, and A. Schwenk,Phys. Rev. C , 054307 (2016).[29] I. Tews, S. Gandolfi, A. Gezerlis, and A. Schwenk, Phys-ical Review C , 024305 (2016).[30] M. Piarulli, I. Bombaci, D. Logoteta, A. Lovato,and R. Wiringa, Phys. Rev. C , 045801 (2020),arXiv:1908.04426 [nucl-th].[31] C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev.Lett. , 042501 (2019).[32] C. Drischler, R. Furnstahl, J. Melendez, and D. Phillips,arXiv:2004.07232 (2020).[33] K. Hebeler, J. M. Lattimer, C. J. Pethick, andA. Schwenk, Phys. Rev. Lett. , 161102 (2010). [34] G. Raaijmakers, S. K. Greif, T. E. Riley, T. Hinderer,K. Hebeler, A. Schwenk, A. L. Watts, S. Nissanke,S. Guillot, J. M. Lattimer, and R. M. Ludlam, Astro-phys. J. Lett. , L21 (2020).[35] C. D. Capano, I. Tews, S. M. Brown, B. Margalit, S. De,S. Kumar, D. A. Brown, B. Krishnan, and S. Reddy,Nature Astron. , 625 (2020).[36] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen,Phys. Rev. Lett. , 172703 (2018).[37] E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, arXiv:1803.00549 (2018).[38] I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C ,045804 (2018).[39] Y. Lim, J. W. Holt, and R. J. Stahulak, Phys. Rev. C , 035802 (2019).[40] S. Greif, K. Hebeler, J. Lattimer, C. Pethick, andA. Schwenk, arXiv:2005.14164 (2020).[41] K. Hebeler, J. M. Lattimer, C. J. Pethick, andA. Schwenk, Astrophys. J. , 11 (2013).[42] S. Greif, G. Raaijmakers, K. Hebeler, A. Schwenk, andA. Watts, Mon. Not. Roy. Astron. Soc. , 5363 (2019).[43] I. Lagaris and V. Pandharipande, Nuclear Physics A ,470 (1981).[44] J. Margueron and F. Gulminelli, Phys. Rev. C , 025806(2019).[45] Y. Lim and J. W. Holt, Phys. Rev. C , 065805 (2017).[46] B. P. Abbott et al. (LIGO Scientific Collaboration andVirgo Collaboration), Phys. Rev. Lett. , 161101(2017).[47] S. De, D. Finstad, J. M. Lattimer, D. A. Brown,E. Berger, and C. M. Biwer, Phys. Rev. Lett. ,091102 (2018).[48] https://dcc.ligo.org/LIGO-P1800115/public .[49] M. Miller et al. , Astrophys. J. , L24 (2019).[50] T. E. Riley et al. , Astrophys. J. , L21 (2019).[51] D. Kasen, B. Metzger, J. Barnes, E. Quataert, andE. Ramirez-Ruiz, Nature (London) , 80 (2017).[52] J. M. Lattimer and M. Prakash, Astrophys. J. , 426(2001).[53] M. Tsang, W. Lynch, P. Danielewicz, and C. Tsang,Phys. Lett. B , 533 (2019).[54] V. F. Suleimanov et al. , Mon. Not. Roy. Astron. Soc. , 906 (2016).[55] H. T. Cromartie et al. , Nature Astronomy (2019),arXiv:1904.06759.[56] A. G. Lyne, M. Burgay, M. Kramer, A. Possenti,R. Manchester, F. Camilo, M. A. McLaughlin, D. R.Lorimer, N. D’Amico, B. C. Joshi, J. Reynolds, andP. C. C. Freire, Science , 1153 (2004).[57] J. M. Lattimer and B. F. Schutz, Astrophys. J.629