Limiting masses and radii of neutron stars and their implications
Christian Drischler, Sophia Han, James M. Lattimer, Madappa Prakash, Sanjay Reddy, Tianqi Zhao
IINT-PUB-20-035
Limiting masses and radii of neutron stars and their implications
Christian Drischler,
1, 2, ∗ Sophia Han,
1, 3, † James M. Lattimer, ‡ Madappa Prakash, § Sanjay Reddy, ¶ and Tianqi Zhao ∗∗ Department of Physics, University of California, Berkeley, CA 94720, USA Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA Department of Physics and Astronomy, Ohio University, Athens, OH 45701, USA Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA Institute for Nuclear Theory, University of Washington, Seattle, WA 98195, USA (Dated: September 12, 2020)We combine equation of state of dense matter up to twice nuclear saturation density ( n sat =0 .
16 fm − ) obtained using chiral effective field theory ( χ EFT), and recent observations of neutronstars to gain insights about the high-density matter encountered in their cores. A key element inour study is the recent Bayesian analysis of correlated EFT truncation errors based on order-by-order calculations up to next-to-next-to-next-to-leading order in the χ EFT expansion. We refine thebounds on the maximum mass imposed by causality at high densities, and provide stringent limitson the maximum and minimum radii of ∼ . (cid:12) and ∼ . (cid:12) stars. Including χ EFT predictionsfrom n sat to 2 n sat reduces the permitted ranges of the radius of a 1 . (cid:12) star, R . , by ∼ . R . < . c s > / n sat , or that χ EFT breaks down below 2 n sat . We also commenton the nature of the secondary compact object in GW190814 with mass (cid:39) . (cid:12) , and discussthe implications of massive neutron stars > . (cid:12) (2 . (cid:12) ) in future radio and gravitational-wavesearches. Some form of strongly interacting matter with c s > .
35 (0 .
55) must be realized in thecores of such massive neutron stars. In the absence of phase transitions below 2 n sat , the small tidaldeformability inferred from GW170817 lends support for the relatively small pressure predicted by χ EFT for the baryon density n B in the range 1 − n sat . Together they imply that the rapid stiffeningrequired to support a high maximum mass should occur only when n B (cid:38) . − . n sat . I. INTRODUCTION
The maximum mass, M max , and radii of neutron stars(NSs) are related to each other by the equation of state(EOS) of dense matter and both can be accessed by ob-servations. Primary constraints on M max come from ob-servations and have a number of astronomical and phys-ical implications. M max is predominately determined bythe EOS at densities higher than three times nuclear sat-uration density, n sat (cid:39) .
16 fm − [1], and is therefore aprobe of the nature of high-density matter. Pinning down M max enables the exploration of the phases of cold anddense matter in the strongly coupled region of quantumchromodynamics (QCD) as well as the determination ofthe pressure vs energy density relation (or the EOS) ofsuch phases. The radii of canonical NSs with masses (cid:39) . (cid:12) , on the other hand, are largely determined bythe EOS at densities less than 3 n sat [2]. M max also fixes the minimum mass of a stellar mass O (M (cid:12) ) black hole (BH). It is therefore a crucial factorin determining the final fate of core-collapse supernovaeand binary neutron star (BNS) mergers. In core-collapsesupernovae, the formation of a BH will depend on the ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] ∗∗ [email protected] amount of fall-back matter and will be sensitive to the na-ture of the progenitor and neutrino emission after the ini-tial formation of a proto-neutron star. In BNS mergers,the formation of a BH depends on the total inspirallingmass, mass ejection, and the extent of rotational andmagnetohydrodynamic support [3, 4]. Now that at leasta few mergers involving NS have been detected throughgravitational-wave (GW) radiation, and many more areanticipated in the near future, improved constraints on M max will become available. As the high-frequency ca-pabilities of GW detectors are improved, the detectionof post-merger radiation will profoundly influence ourknowledge of M max . Already, knowledge of M max woulddetermine the nature of the components of the recentlyobserved mergers GW190425 and GW190814, both ofwhich show indications of having a component with amass larger than 2 M (cid:12) that either could be a heavy NSor a light BH. If concomitant electromagnetic (EM) sig-nals are also detected from future GW events, as theywere in the BNS merger GW170817 [5–7], additional in-formation about M max becomes available [3, 4].On the theoretical front, M max plays a crucial rolein determining both the minimum and maximum ra-dius as a function of the NS mass M . Therefore, be-sides the important contributions from radio and X-raybinary pulsar observations that have accurately mea-sured several NS masses and provided a lower bound M max (cid:38) (cid:12) [8–12], GW and X-ray data that can simul-taneously determine NS masses and radii offer importantconstraints. So far, the radii inferred from X-ray observa- a r X i v : . [ nu c l - t h ] S e p tions (see Ref. [13] for a review) of quiescent low-mass X-ray binaries (QLMXBs) [14], photospheric radius expan-sion bursts (PREs) [15], and pulse-profiles from rotation-powered millisecond pulsars [16], together with the firstGW detection of the BNS merger GW170817 [5, 6],have mostly been of NSs with canonical masses around1 . (cid:12) . Consequently, the Neutron Star Interior Com-position ExploreR (NICER) proposal [17] to measurethe radii of relatively massive NSs such as PSR J1614-2230 ( M (cid:39) .
91 M (cid:12) [8, 10, 11]) and PSR J0740+6620( M (cid:39) .
14 M (cid:12) [12]) is of considerable interest. The sameis true of future radio observations using the
Square Kilo-metre Array (SKA) telescope [18], etc. from binary pul-sars that could reveal even more massive NSs.The purpose of this paper is to explore the interplaybetween M max and NS radii and to confront theoreti-cal expectations with currently available observationalconstraints. An earlier study [2] showed that the radiiof (cid:39) . (cid:12) NSs are strongly correlated with the pres-sure of matter in the density range 1 − n sat . In theimportant density range (cid:46) n sat , chiral effective fieldtheory ( χ EFT) with pion and nucleon degrees of free-dom has become the dominant microscopic approach todescribing nuclear interactions [19–22]. Applying χ EFTto the EOS of infinite nuclear matter and the structureof neutron stars [23–37] has enabled significant progress(see also Refs. [38–40] for recent reviews). A notable de-velopment is the quantification of theoretical uncertain-ties in a statistically robust way. Recently, the BUQ-EYE collaboration [41] has introduced a novel frame-work for quantifying correlated EFT truncation errorsin infinite-matter calculations. In a first application,they have conducted a statistical analysis of the zero-temperature EOS with χ EFT nucleon-nucleon (NN) andthree-nucleon (3N) interactions up to next-to-next-to-next-to-leading order (N LO) [42, 43] using GaussianProcesses (GPs). This study was motivated by recent ad-vances in many-body perturbation theory (MBPT) [35]that have enabled improved χ EFT predictions of the pureneutron matter (PNM) EOS at this order, and, for thefirst time, order-by-order calculations in symmetric nu-clear matter (SNM) based on χ EFT NN and 3N forcesup to N LO [35, 42, 44].In this paper, we use BUQEYE’s analysis of the EOSin the limits of PNM and SNM up to 2 n sat to constructthe EOS of charge neutral and beta-stable neutron-starmatter (NSM). This is coupled to a standard NS crustfor densities (cid:46) . n sat and extrapolations for densities (cid:38) . n sat to assess the overall impact on NS structure.One goal of this study is to address quantitatively theextent to which EOS knowledge at ∼ . n sat can in-form us about the NS maximum mass, and how it can becombined with observations of massive NSs to constrainthe properties of matter encountered at the highest den-sities in their cores. Another goal is to derive model-independent bounds on the radii of NSs with masses inthe range 1 − (cid:12) . As the squared speed of sound c s reflects the stiffnessof the EOS, we probe both maximum and minimum ra-dius bounds by matching the N LO results, includingpossible extrapolations up to 3 . n sat , with a constantsound speed beyond a matching density n m . The exis-tence of nuclei, observations of accreting NSs that impli-cate the presence of neutron-rich nuclei in the NS crust,and heavy ion collisions (HICs) at intermediate energiestogether provide compelling circumstantial evidence toindicate that n m > n sat , and in this work we consider1 . ≤ n m /n sat ≤ .
0. The use of the maximally stiffEOS with c s = 1 (the causal limit) for n B > n m estab-lishes firm upper bounds both on M max and the radius asa function of mass. In addition, we also consider energydensity discontinuities at n m to refine minimum boundson radii as functions of mass for specified values of M max .We also explore models with smaller c s at high densityto ascertain maximum possible sound speeds from valuesof M max and mass-radius ( M – R ) observations.The discovery of a massive secondary compact objectwith mass ∼ . (cid:12) through GW observations of thebinary merger GW190814 generated a flurry of articlesaddressing if this object can be a NS, and, if so, its possi-ble implications [45–52]. Our results complement earlierstudies, but go beyond in several aspects. Most signifi-cantly,(i) we consistently include statistically meaningfulEFT truncation errors in the EOS of NSM up toN LO, and determine its range of applicability, toprovide a framework for constraining M max and NSradii,(ii) we identify correlations of NS radii and tidal de-formabilities with M max , together with their possi-ble implications for the EOS at densities (cid:38) . n sat ,and(iii) we show how these correlations and future obser-vations can tighten current bounds on NS massesand radii.This paper is organized as follows. In Sec. II we discussthe bounds imposed by causality and the scaling relationsfor the masses and radii of NSs. Section III contains de-tails of the various EOSs used along with the rationalefor their choice. Our results and their discussion in lightof the current observational constraints and possible fu-ture findings are presented in Sec. IV. An overall dis-cussion and comparison with pertinent recent works arecontained in Sec. V. Our concluding remarks are given inSec. VI. We use natural units in which (cid:126) = c = 1 unlessexplicitly specified. II. BOUNDS IMPOSED BY CAUSALITY
The assumption of causality, i.e., that the maximumsound speed c s = (cid:112) dP/dε is unity in units of c , can es-tablish relations limiting both minimum and maximumradii, as functions of mass, for NS. These limits will ex-plicitly depend on assumptions concerning the NS maxi-mum mass M max . These causal bounds can be improvedwith the consideration of nuclear physics inputs as willbe discussed in Sec. IV. The causality limit is imposedby using the EOS P ( ε ) = P + ( ε − ε ) (1)for the pressure P > P and the energy density ε > ε .The minimum radius as a function of mass R min ( M )for any EOS is conjectured [53] to result from usingEq. (1) with P = 0, P = 0 for ε ≤ ε (i.e., a self-bound star). In this case, the EOS has a single parameter( ε ) and solutions of the Tolman-Oppenheimer–Volkoff(TOV) equation [54, 55] scale with it. Letting m be themass enclosed within the radius r , one can define r = x c √ Gε , m = y c √ G ε , and P = zε , (2)where y ( x ) and z ( x ) are dimensionless functions, with theboundary conditions y c = y ( x = 0) = 0 and z c = z ( x =0) > y s = y ( x = x s ) and z ( x = x s ) = 0 at the stellar surface x s . The quantities y s and x s depend on z c . For small x s , y s ∝ x s , as expected.It should also be noted that the EOS Eq. (1) implies thatthe baryon number density is n B = n (cid:114) ε + Pε + P , (3)with n = ( ε + P ) /µ and µ being the baryon chemicalpotential at ε .In the case that P = 0, the central baryon densityis n cent = n √ z c . Also, the maximum mass con-figuration occurs for dy s /dx s = 0, or when x max , s =0.2405, y max , s = 0 . z max , c = 2 .
023 (and there-fore n max , c /n = 2 . M max = y max , s c √ G ε (cid:39) . (cid:114) ε sat ε M (cid:12) , (4)and the radius of the maximum mass configuration is R M max = x max , s c √ Gε (cid:39) . (cid:114) ε sat ε km . (5)The central energy density for the maximum mass con-figuration is ε max , c = ( z max , c + 1) ε , or using Eq. (4) toeliminate ε , ε max , c (cid:39) . (cid:18) M (cid:12) M max (cid:19) ε sat , (6)where ε sat (cid:39)
150 MeV fm − is the energy density at n sat .This must be the largest energy density found in any NS and it scales with M − . The maximum baryon densityis n max , c (cid:39) . m B µ (cid:18) M max M (cid:12) (cid:19) n sat , (7)where µ ∼ m B , the baryon mass. As an example, ifone assumes that M max = 2 . (cid:12) and µ = m B , itis found that ε = 2 . ε sat , ε max , c = 7 . ε sat and n max , c = 5 . n sat . FIG. 1. The mass as a function of radius for the EOS Eq. (9)with P = 0, P = 0 for ε < ε , and various values for c s withfixed M max = 2 . (cid:12) , are shown as five black curves (seelegend). These curves correspond to the minimum possibleradius R min ( M ), for different maximum values of the soundspeed. The four red curves correspond to ε = ε sat , either c s = 1 (and M max (cid:39) .
09 M (cid:12) ) or c s = 1 / M max (cid:39) .
48 M (cid:12) ) for
P > P , and either P = 0 (self-bound) or P =0 . ε (cid:39) − and a normal crust EOS for P < P (maximum possible radii R max ( M )); the configuration where ε c = ε sat is indicated by a diamond. The dimensionless M – R curve for the causal self-bound configuration is thus defined by y s ( x s ). Its dimen-sionful radius, as a function of mass, is conjectured to bethe minimum radius for any configuration, R min ( M ). Itscales with ε and therefore with the assumed value ofthe maximum mass: R min = GM max x s y max , s c = GM max y max , s c y − s (cid:18) y max , s MM max (cid:19) , (8)and increases as M max increases. The cases with M max =2 . (cid:12) and M max = 4 .
09 M (cid:12) for which ε is 2 . ε sat and ε sat , respectively, are shown in Fig. 1. For the casethat M max = 2 . (cid:12) for which ε = 4 . ε sat , we obtain R min (1 . (cid:12) ) = 8 . R M max = 8 . TABLE I. Maximum mass solutions for the EOS Eq. (9) with P = 0. The last two columns give the minimum radii in kmfor 1 . (cid:12) and 2 . (cid:12) stars, respectively, assuming M max =2 . (cid:12) . c s x max ,s y max ,s z max ,c R min , . R min , . If the assumed maximum sound speed is less than c , R min ( M ) will increase. Assuming the sound speed neverexceeds a given value of c s , R min ( M ) can be found using P = P + c s ( ε − ε ) , (9)with P = 0 and P = 0 for ε < ε . Once again, theTOV equation can be rendered into dimensionless formusing Eq. (2). Now, however, the baryon number densitybecomes n B = n (cid:18) P + εP + ε (cid:19) / (1+ c s ) (10)and n cent = n (cid:2) z c (cid:0) c − s (cid:1)(cid:3) / (1+ c s ) . (11)The dimensionless M – R curve y s ( x s ) changes, as do theproperties of the maximum mass configuration x max , s , y max , s and z max , c . Figure 1 shows M – R solutions for c s = 1 , / , / , /
2, and 1 /
3, all scaled so that M max =2 . (cid:12) . In the case c s = 1, one finds R min , . = 9 .
75 km and R min , . = 10 . . (12)Approximately, the minimum radii for smaller values of c s scale as c − / s [56], and for c s = 1 /
3, one finds that R min , . (cid:39) . R min , . (cid:39) . . (13) ε max , c /ε sat is proportional to z max , c + 1, which for c s < c / s . Relevant properties ofthese solutions are given in Table I.Stars with P = 0 are often referred to as self-bound stars. In contrast, normal NSs have a low-density crustwith P >
0. For normal stars, R min ( M ) will be largerthan those shown in Fig. 1. Generally, the radius willincrease with the assumed values of ε and P for a givenvalue of c s , and, to a lesser degree, will also depend onthe crust EOS for P < P . Most importantly, since M max and ε remain closely related, R min ( M ) will be very sensi-tive to the lower limit to M max . Details and implicationsare discussed in Sec. IV C.Ironically, the maximum radius as a function of mass R max ( M ) can also be found by appending the same EOS Eq. (9) at a matching density n m or ε m onto an as-sumed lower-density (crust) EOS. This is because Eq. (9)is the stiffest possible EOS for an assumed maximumvalue of the sound speed c s . Although the same EOS isused, the R min ( M ) bound involves a finite surface energydensity ε = ε m , while the R max ( M ) bound is assumed tolack a discontinuity in ε when appending the crust . Theresulting R max ( M ) trajectory, and M max , will depend onthe matching density ε m and pressure P m , the crust EOS,and assumed maximum sound speed c s , and both roughlyscale as ε − / . Since there is no evidence that a transi-tion to a non-hadronic EOS occurs for densities smallerthan ε sat , a limiting set of R max ( M ) curves is found as-suming ε m = ε sat . As the matching pressure P m is notnegligibly small, P m (cid:39) . ε sat for ε m (cid:39) ε sat , the M − R curve is considerably altered, and forms a maximum ra-dius trajectory R max ( M ) which lies at a larger radius foreach mass than R min ( M ), as can be seen by comparingthe two solid red curves for c s = 1 in Fig. 1. R max for c s = 1 can be safely assumed to give, approximately, thelargest possible radii for normal NS (it varies with theassumed EOS below ε m ). It is interesting that the max-imum masses with ε m = ε sat for a self-bound star (leftred solid curve) and for a normal star with a crust (rightred solid curve) are nearly identical and are substantiallylarger than 2 . (cid:12) , for example. Lower maximum massesare obtained if the matching density is increased, whichdecreases R max ( M ) as well. An observed upper limit on M max below 4 .
09 M (cid:12) will automatically alter the R max boundary, however, because in this case either ε m wouldhave to increase or c s would have to decrease to corre-spondingly reduce M max .The situation is similar if a lower fixed sound speed isassumed. Figure 1 also displays R min ( M ) and R max ( M )trajectories for c s = 1 / M max val-ues than for c s = 1. Note that R max ( M ) for c s = 1 / R min ( M ) for c s = 1 and P = 0 (left red solid curve) for M (cid:38) . (cid:12) , suggesting that c s = 1 / M max = 2 . (cid:12) ; the maxi-mum value of c s must be larger than 1 / . (cid:12) star, or P > M max < . (cid:12) .A more realistic maximum radius boundary will de-pend on both the matching density and the EOS belowthat density. In the next section we discuss realistic con-straints on this portion of the EOS stemming from theo-retical studies of NSM. Note that if a discontinuity in ε is assumed at ε m , a smaller R max trajectory is obtained, but one with a correspondingly smaller M max as well. This situation is briefly discussed in Sec. IV B. III. CONSTRUCTION OF THE EOS
The EOS up to the core-crust boundary of NSs at ≈ . n sat is generally considered to be well-understood [57,58]. Because nucleons contribute (cid:46)
10% to the crustpressure, uncertainties in the NN potential only weaklypropagate into the crust EOS. The proton fraction inmatter at densities higher than the core-crust bound-ary is relatively small, so that the EOS in the vicinityof the saturation density is dominated by that of PNM.The admixture of protons and leptons produce small cor-rections, which are effectively minimized because of therequirement that NSM be in beta equilibrium; that is,the total energy is minimized with respect to the protonfraction.During the past few years, there have been impor-tant theoretical advances in the understanding of densenuclear matter in the density regime (cid:46) . n sat from χ EFT studies. χ EFT [19–22] exploits the separation ofscales that arises due to the gap between the massesof the (pseudo-)Goldstone bosons of chiral symmetrybreaking—the pions—and energy scales associated withinteractions at shorter length scales. While the long-range pion exchanges are explicitly resolved in χ EFT,short-range interactions are given by contact interactionswhose low-energy couplings need to be fit to experimentaldata. The momentum scale associated with the interac-tions at short-distances is denoted by Λ b , and is called thebreakdown scale of the EFT. When the relative momen-tum between nucleons is small compared to Λ b , χ EFTaims to provide a model-independent and systematicallyimprovable description of nuclear interactions and ob-servables.Weinberg power counting in χ EFT organizes the mostgeneral nuclear Hamiltonian that is consistent with allsymmetries of low-energy QCD in powers of Q =max { m π , p } / Λ b , where p is a typical momentum associ-ated with the interacting nucleons (soft scale) and Λ b themomentum at which the χ EFT expansion breaks down(hard scale) because the details of the short-distance in-teraction become relevant [19–22]. If p < Λ b , χ EFTcalculations can (in principle) be improved up to thedesired accuracy by working to higher orders; and theresidual uncertainties due to truncation of the χ EFTexpansion at finite order can be quantified. The con-vergence of the χ EFT has been studied in detail in thecontext of free-space NN scattering using Bayesian meth-ods [59, 60]. These studies indicate that Λ b ≈
600 MeV—consistent with the phenomenological expectation thatvector mesons such as the ρ − meson with mass m ρ (cid:39)
770 MeV contribute to nuclear forces at short distances.In dense matter, the typical momentum scale associ-ated with nuclear interactions will in general depend onthe nature of the ground state. When interactions arenot too strong, the ground state is expected to resemblea Fermi liquid, so the typical momentum associated withinteractions will be proportional to the Fermi momentum k F = (6 π n B /g ) / , where n B is the baryon number den- sity and g the spin-isospin multiplicity (see Refs. [42, 43]for details). Thus, with increasing density we should ex-pect the efficacy of the χ EFT to diminish as the trunca-tion errors associated with finite order calculations grow.An important recent development is the consistentquantification and propagation of these EFT truncationerrors in the energy per particle and derived quantitiesin dense matter using Bayesian methods [42, 43]. Themany-body calculations used in these studies were con-ducted order-by-order in the χ EFT expansion, which al-lowed for Bayesian inference of the a priori unknownΛ b in infinite matter by making verifiable assumptionsabout the χ EFT expansion that can be tested using di-agnostic tools for model checking [59, 60]. This in turnprovides the systematic means to estimate and propagateEFT truncation errors. As noted in the previous section,realistic determinations of both R min ( M ) and R max ( M )depend on the effective boundary ε as well as the pres-sure of nucleonic matter at P < P down to the crust.The χ EFT uncertainties will therefore contribute to theuncertainties in R min ( M ) and R max ( M ) for an assumedvalue of ε .In what follows, we describe the results of the many-body calculations, the associated EFT truncation erroranalysis, and the propagation of theoretical uncertaintiesused in constructing the EOS of beta-equilibrated NSM. A. Many-body results for PNM and SNM
The BUQEYE collaboration [41] has recently intro-duced a novel framework for uncertainty quantification ofinfinite-matter observables calculated from EFT [42, 43].In this approach, Gaussian Processes with physics-basedhyper-parameters are trained on order-by-order calcula-tions of the EOS as a function of the baryon density andproton fraction. This allows for quantification and prop-agation of statistically robust uncertainties of the EOSto quantities involving derivatives, while self-consistentlyaccounting for correlations in density and across observ-ables (see Ref. [43] for details). Nowadays, EFT trunca-tion errors dominate these theoretical uncertainties, butmany-body approximations, etc., also contribute.As the framework’s first application [42, 43], BUQEYEstudied the χ EFT convergence of the EOS in the limitsof PNM and SNM up to N LO in the density region n B (cid:46) . n sat . They also inferred posterior distributionsfor nuclear saturation properties and key quantities forneutron stars: the pressure and speed of sound in PNM,and the nuclear symmetry energy as well as its slope indensity.The statistical analysis indicates that the EOS isstrongly correlated, which leads to correlation lengthscomparable to the k F associated with n sat in PNM andSNM, respectively. In other words, perturbing the EOSat one point in the density (or the proton fraction) per-turbs neighboring points (with respect to the correlationlength) as well. Without including these correlations, un- .
10 0 .
15 0 .
20 0 .
25 0 . n B [fm − ]051015202530 P r e ss u r e P N S M [ M e V f m − ] (a) N LON LO 0 .
10 0 .
15 0 .
20 0 .
25 0 . n B [fm − ] − . − . − . − . . . . D i ff e r e n c e i n P r e ss u r e ∆ P [ M e V f m − ] (b)∆ P ( n B ) = P PNM ( n B ) − P NSM ( n B ) FIG. 2. Panel (a): pressure of neutron-star matter (NSM) as a function of the baryon number density at N LO (orange-shadedband) and N LO (blue-shaded band) in the chiral expansion; panel (b): differences of the PNM and NSM pressures as annotatedin the panel. The uncertainty bands depict the 1 σ confidence region. certainties in derived quantities of the EOS, such as thenuclear symmetry energy, can be overestimated.Assuming p = k F , the inferred χ EFT breakdown scalein PNM and SNM, Λ b ≈ ±
50 MeV, was found tobe consistent with free-space NN scattering. If no newdegrees of freedom emerge at high densities, this could beassociated with n B ( k F = Λ b ) = Λ b / (3 π ) ≈ . ± . n sat in PNM, although the χ EFT expansion will start losingits predictive power at much lower densities. This should,however, only be considered as a rough estimate (withlarge uncertainties).Specifically, BUQEYE’s analysis is based on recentorder-by-order MBPT calculations in PNM and SNMwith chiral NN and 3N interactions up to N LO [35,42, 44]. The range in density covers n B = 0 . − .
34 fm − . These calculations significantly improved pre-vious MBPT studies in PNM at N LO [23, 33, 61], andassessed, for the first time, the SNM EOS with NN and3N interactions order-by-order up to N LO. The high-order MBPT calculations were performed by the novelMonte Carlo framework introduced in Ref. [35], whichenables MBPT calculations of the EOS with controlledmany-body uncertainties for these χ EFT interactions.The underlying nuclear interactions were constrainedin Ref. [35] as follows: NN potentials by Entem, Mach-leidt, and Nosyk [62] up to N LO were combined with3N forces at the same order and momentum cutoff so asto construct a set of order-by-order NN and 3N interac-tions. The two 3N low-energy couplings c D and c E , whichgovern the intermediate- and short-range 3N contribu-tions, respectively, at N LO were constrained by the tri-ton binding energy and the empirical saturation point ofSNM. Several combinations of c D and c E with reasonablesaturation properties could be obtained at N LO and N LO for the momentum cutoffs Λ = 450 and 500 MeV.A momentum cutoff is a typical scale in the regulatorfunction that is applied to χ EFT interactions to sup-press contributions from high-momentum modes. Notethat Λ b is a physical scale inherent to the EFT, whereasthe results should not be sensitive to the artificial scaleΛ; in practice, however, this has not yet been achieved in χ EFT for infinite matter. The BUQEYE collaborationfound that their results do not significantly dependenton which c D and c E combination is chosen for a givenmomentum cutoff. Furthermore, the 3N contributionsproportional to c D and c E vanish in PNM for nonlocalregulator functions [63]. Consequently, they consideredonly one combination for each cutoff, and focused theiranalysis on the Hamiltonian with Λ = 500 MeV, while theresults for the Λ = 450 MeV interaction were provided inthe Supplemental Material there.We follow this strategy here, and note that theresidual cutoff dependence is well within the EFTtruncation-error estimates at the 1 σ level; i.e., forΛ = 450 MeV, P PNM (2 . n sat ) = 17 . ± .
56 MeV fm − and E PNM (2 . n sat ) = 42 . ± .
01 MeV, whereas forΛ = 500 MeV, P PNM (2 . n sat ) = 18 . ± .
14 MeV and E PNM (2 . n sat ) = 41 . ± .
77 MeV fm − .BUQEYE’s EOSs are given as GPs and publicly avail-able as Jupyter notebooks [41]. A GP is an infinite di-mensional generalization of a multi-variate normal dis-tribution, where each point in density is associated witha random variable. Using the Jupyter notebooks we ex-tract the mean values, standard deviations, and correla-tion information of the energy per particle, pressure, andspeed of sound in PNM and SNM, and also the symmetryenergy. R [km]0 . . . . . . . M a ss M [ M (cid:12) ] n m = 2 . n sat χ EFT(a) χ EFT + matched linear EOS at c s, match = 1 9 10 11 12 13 14 15 16Radius R [km]0 . . . . . . . M a ss M [ M (cid:12) ] n m = 1 . n sat χ EFT(b) χ EFT + matched linear EOS at c s, match = 1N LON LOself-bound, c s = 1 FIG. 3. Panel (a): M – R diagram for NSM based on MBPT calculations shown in Fig. 2 (a), including N LO-NSM (orange-shaded band) and N LO-NSM (blue-shaded band) at low densities (cid:46) . n sat (both with ± σ uncertainties), matched to alinear causal ( c s, match = 1 .
0) EOS at n m = 2 . n sat , and the maximally compact EOS for self-bound stars with the same valueof M max (black solid). Horizontal lines indicate M = 1 .
4, 2.0, 2 . (cid:12) . The colored bands above n m = 2 . n sat represent upperbounds on the NS radius for a given mass, as the high-density matter is assumed maximally stiff without discontinuities in theoverall EOS (see detailed discussions in Sec. IV C). Panel (b): similar to (a) but with a lower matching density, n m = 1 . n sat . We consider correlations between the EOSs in PNMand SNM explicitly, neglecting correlations in density.These observables are then given by independent normaldistributions sampled on a fine grid in density using theGPs; e.g., E PNM ∼ N (cid:0) µ PNM , σ (cid:1) , (14) E SNM ∼ N (cid:0) µ SNM , σ (cid:1) . (15)The nuclear symmetry energy is defined as E sym = E PNM − E SNM ∼ N (cid:0) µ NSM , σ (cid:1) , (16)and, hence, has mean and variance (see, e.g. Ref. [64]): µ sym = µ PNM − µ SNM , (17) σ = σ + σ − ρσ PNM σ SNM , (18)where ρ is the correlation coefficient between the ener-gies per particle in PNM and SNM. For subsequent dis-cussion, we introduce here also the usual parameters S v and L in the density expansion of the nuclear symmetryenergy Eq. (16), E sym = S v + L (cid:18) n B − n sat n sat (cid:19) + . . . . (19)The correlation between the coefficients in the χ EFTexpansions for the PNM and SNM energy per particlewas quantified to be ρ ∗ = 0 . verystrong correlations [65, 66]. A detailed discussion can befound in Sec. IV A of Ref. [43]. We have checked that ρ (cid:39) ρ ∗ by comparing E sym against the values obtained in Ref. [43]: the maximum deviation between the twoapproaches is 37 keV (340 keV for its ± σ bounds) atthe highest density, n B = 0 .
34 fm − , which is negligiblecompared to the overall EFT truncation error at thatdensity.We also found that numerical integration of the pres-sure of PNM and SNM agreed well with the energy foundin the GP approach, the maximum deviation being 3 keVand 1 keV for PNM and SNM, respectively (290 keV and500 keV for their respective ± σ bounds) at the high-est density. There are mainly two related reasons whyfinite differencing for the pressure, discrete integrationfor the energy, and subtraction for the symmetry energy,works so well. First, the correlation length of the EOSis much longer than the length scale used for finite dif-ferencing. That means numerical differentiation followsclosely the curves µ ± σ , which are two realizations ofthe underlying GP. Secondly, the raw EOS data has al-ready been preprocessed by BUQEYE’s truncation errormodel. Numerical noise from the many-body method hasbeen smoothed out, and the EOS has been sampled on afine grid in density using the GP interpolant. This under-lines that GP interpolants are efficient tools for analyzing χ EFT calculations of the EOS.
B. Uncertainty propagation to NSM
To construct the EOS of charge-neutral, beta-equilibrated NSM between the crust and n B (cid:46) . n sat ,we use the standard approximation of keeping only thequadratic term in the nuclear energy expanded in the FIG. 4. Panel (a): M – R diagram for matched linear EOSs that give rise to M max = 2 . (cid:12) with N LO-NSM ( ± σ ) appliedfor low densities ≤ . n sat . Corresponding values of n m and c s, match are indicated. Panel (b): sound speed profiles c s ( n B )for N LO-NSM only (black-solid for the central value and black-dashed for ± σ uncertainties), and matched linear EOSs withdifferent values of c s, match associated with M max = 2 . (cid:12) in panel (a) (colored horizontal lines). The open triangles mark thecentral densities of the maximum-mass stars M max = 2 . (cid:12) . isospin asymmetry parameter β = 1 − x , where x = n p /n B is the proton fraction, and n p the proton density.The total energy of NSM is then E NSM = E PNM (1 − x ) + E SNM x (1 − x )+ E e + E µ , (20)where E e and E µ are the energies per baryon of elec-trons and muons, respectively. Microscopic calculationsof asymmetric matter based on chiral NN and 3N inter-actions at n B (cid:46) n sat have confirmed that the quadraticexpansion Eq. (20) is a reasonable approximation of thefull isospin dependence of the EOS [26, 67–70].Beta equilibrium is the condition that the total charge-neutral energy be minimized with respect to x , or ∂E NSM ∂x = 0 , (21)which is equivalent to µ n − µ p = 4 E sym (1 − x ) = µ e = µ µ . (22)Here, µ denotes a chemical potential. Propagating theEFT uncertainties to E NSM associated with Eq. (20) isstraightforward because of the condition (21). We obtain σ E NSM = (cid:18) ∂E NSM ∂E PNM (cid:19) σ E PNM + (cid:18) ∂E NSM ∂E SNM (cid:19) σ E SNM + 2 ρ ∂E
NSM ∂E PNM ∂E NSM ∂E SNM σ E PNM σ E SNM , (23)with the derivatives ∂E NSM /∂E
PNM = (1 − x ) and ∂E NSM /∂E
SNM = 4 x (1 − x ). Figure 2 (a) shows the pressure of NSM (including con-tributions from the leptons) P NSM = n ( dE NSM /dn B ).The blue (orange) uncertainty band corresponds to theN LO (N LO) results at the 1 σ level. Panel (b) displaysthe difference in pressures between PNM and NSM. Thezero crossings indicate where the pressure of NSM equalsthat of PNM. Depending on the chiral order, these cross-ings occur at n ≈ . − . n sat . They are due to a soften-ing of E sym at the higher densities; nevertheless, E NSM isalways less than that of E PNM . In no case does x exceedabout 0 .
055 for n B ≤ .
34 fm − . IV. RESULTSA. NSM and the role of χ EFT inputs
In order to compute bounds on NS masses and radii, weassume a typical crust EOS [57, 58] below 0 . n sat and theEOS for NSM based on MBPT- χ EFT calculations [43] upto n m = 1 . − . n sat , and then match to a linear EOSEq. (9) characterized by c s at higher densities.Figure 3 (a) shows the M – R relation for N LO-NSM,N LO-NSM EOSs in Fig. 2 (a) matched at n m = 2 . n sat to the stiffest linear EOS that with c s = 1; the solidcolored curves refer to the central values and the color-shaded bands refer to ± σ uncertainties in the MBPT- χ EFT calculations. Results for matching at a lower den-sity n m = 1 . n sat are shown in Fig. 3 (b). As expected,for a given value of n m , the largest radii result from the FIG. 5. ˜Λ– M and Λ– M relations confronted with constraints from GW170817 [6, 7] (vertical lines with arrows), with fixed M max = 2 . (cid:12) as an example. Parameters for matched EOSs are the same as in Fig. 4, except for the special case with n m = 1 . n sat (with c s, match = 0 . M = 1 .
186 M (cid:12) ) ≤ largest matching pressure P m , and thus N LO +1 σ , whileN LO-1 σ shows little difference compared to N LO-1 σ ;here − σ (+ σ ) refers to the lower (upper bound) on P ( n B )in Fig. 2. In general, the lower the matching density n m ,the larger the maximum radii. Any discontinuities in theenergy density, such as from a phase transition, wouldserve to decrease R ( M ), emphasizing the results in thesefigures as being upper bounds (see Sec. II). It can beseen that upper bounds on R . (where the bands in-tersect with the M = 1 . (cid:12) horizontal line) of about12 . . n m = 2 . n sat ( n m = 1 . n sat ).We defer a detailed discussion on the radius bounds toSec. IV C. Self-bound (crustless) stars with c s = 1 fora given M max represent the “maximally compact” con-figurations that exhibit the smallest possible radii at allmasses, and for comparison their mass-radius relationsare also displayed (black solid curves).The overall matched EOSs ε ( P ) are continuous buttheir sound speeds are not, and the causal limit c s = 1leads to maximum masses as high as ≈ .
93 M (cid:12) for n m = 2 . n sat and ≈ .
36 M (cid:12) for n m = 1 . n sat (differ-ences at low densities, e.g., between N LO and N LOhave negligible effects on M max , as already noted inSec. II). For a given value of c s, match , M max is essen-tially determined by n m and is relatively insensitive to P m . With smaller (and constant) values of c s, match , themaximum mass decreases.It is therefore of interest to examine what matchingconditions relating n m and c s, match ensue from a restric-tion such as M max = 2 . (cid:12) . Fig. 4 (a) depicts the M – R relations for n m and c s, match that lead to M max = 2 . (cid:12) , and the corresponding c s profiles are explicitly shown inpanel (b). The required values of c s, match are indicated inthe plot (solid horizontal lines for the N LO-central (de-noted as N LO-cen) and dashed for ± σ uncertainties),which increase with the matching density n m . At fixedmatching density indicated by the vertical dotted lines,the variation in c s, match above n m is consistent with theuncertainties in c s from χ EFT calculations at n m , and asofter EOS (smaller c s ) at low densities is compensatedby a stiffer EOS (larger c s, match ) at higher densities.The simple linear parametrization of high-density EOSused here can be viewed as a guide to assess the stiffnessrequired at higher densities to achieve M max ≥ . (cid:12) .Assuming χ EFT-N LO is valid up to n m = 2 . n sat (1 . n sat ), to reach 2 . (cid:12) the “averaged” c s above2 . n sat (1 . n sat ) has to be greater than ∼ . ∼ . c s is gradually increasing) without violating causality belowthe central density of the maximum mass star. On theother hand, the low-density behavior of the EOS con-trols predictions on the radii and tidal deformabilities ofcanonical-mass neutron stars 1 . − . (cid:12) , for which var-ious astrophysical and gravitational-wave constraints areavailable. The relatively soft N LO EOS up to 2 . n sat guarantees that the typical NS radius (cid:46)
13 km, even withvery stiff matter at higher densities ( c s, match ≈ . . n sat ) that leads to M max ≈ . (cid:12) ; see redcurves in Fig. 4 (a).Calculations of tidal deformabilities [71–73] are con-fronted with the constraints inferred from GW170817 [5–0 FIG. 6. Panel (a): the solid lines show contours of n m in the M max – c s, match plane, and the dashed lines bracket N LO ± σ uncertainties. The upper horizontal line indicates M max = 2 . (cid:12) ; see also examples in Fig. 4. The grey-shaded region isexcluded by the binary tidal deformability constraint ˜Λ . ≤
720 from GW170817 at the 90% crediblitiy level [6] if N LO-cenis assumed; the dot-dashed lines refer to constraints with the N LO ± σ boundaries. The thin dotted line indicates a lowerupper bound with N LO-cen and ˜Λ . ≤ c s, match are displayed inthe M max – n m plane; the upper-right grey-shaded region is excluded by causality. For n m ∈ [2 . , . n sat , extrapolations from χ EFT using ZL models with L = 50 MeV and L = 60 MeV are applied (see Fig. 7).
7] in Fig. 5. Gravitational waveform fitting using thestandard PhenomPNRT model [6, 74] directly sets con-straints on the binary chirp mass M = 1 . ± .
001 M (cid:12) and the tidal deformability ˜Λ ≤
720 (90% credibility). Inwhat follows, we will denote the constraint on ˜Λ at M as ˜Λ . ≤ n m results in a large radius for a given c s, match .An EOS stiffening drastically from N LO below 1 . n sat ends up violating ˜Λ . ≤
720 if M max (cid:38) . (cid:12) (greenband in Fig. 5 (a)). The ˜Λ– M constraint can be trans-lated to a constraint on Λ at the mass M , Λ M , but it issubject to small additional uncertainties from the poorlydetermined mass ratio q of GW170817 and EOS system-atics. Using the quasi-universal EOS relation Λ = q Λ ,which is valid to 10% −
20% for M = 1 .
186 M (cid:12) and q > . M (cid:39) / ( M /M ) ˜Λ M , (24)valid to a few percent. Even considering the q andEOS uncertainties, one sees that n m (cid:46) . n sat violatesthe GW170817 constraint (Fig. 5 (b)). There exists aminimum n m ≈ . n sat for N LO-NSM to survive the˜Λ ≤
720 constraint (when M max ≥ . (cid:12) is assumed),and an even smaller upper bound e.g., ˜Λ (cid:39)
600 [76–78]which would increase the minimum required n m . Itis noteworthy that the posteriors of ˜Λ for GW170817suggest a peak value around ≈ . − . (cid:12) , increases the allowed range for n m and c s, match to be consistent with data; the generic trendis shown in Fig. 6. Specifically, Fig. 6 (a) demonstrateshow M max scales with c s, match using the N LO-NSMEOS for n m = 1 . , . , . n sat . The solid curves corre-spond to results for N LO-cen and the dashed ones with ± σ uncertainties. The dots indicate the intersectionsof the central curves with M max = 2 . (cid:12) for thesame EOSs as in Fig. 4 (b). The χ EFT uncertaintiesat the respective densities only slightly broaden thesecorrelations. Together with GW170817, the constraint M max ≥ . (cid:12) rules out very weakly-interacting matter( c s ≈ .
33) at high densities, whereas M max ≥ . (cid:12) rules out matter with c s (cid:46) . LO-cen is nearly parallelto the n m contours. The dot-dashed lines bracket ± σ uncertainties in N LO, and the softer N LO-1 σ (stifferN LO +1 σ ) corresponds to the upper (lower) boundary.This feature is more clearly displayed in panel (b), where M max is plotted against n m for fixed values of c s, match .For matching densities (cid:46) . − . n sat , all constructedEOSs result in ˜Λ . >
720 and can be therefore con-sidered ruled out by GW170817 (see also examples in1
FIG. 7. M – R relations for EOSs extrapolated to high densi-ties beyond the χ EFT calculations using the ZL parametriza-tion [79] for NSM at n B ≥ . n sat ; the thin black line indi-cates where the NS central densities are 3 . n sat . From leftto right, the colored, dotted curves represent L = 45 MeVto L = 75 MeV in increments of 5 MeV, and the black-solid(black-dashed) curves refer to χ EFT-N LO (N LO), includ-ing ± σ uncertainties, up to ∼ . n sat . The L = 50 MeV(red) and L = 60 MeV (green) ZL EOSs are used in Fig. 6 toinvestigate the trend of M max ( n m ) scaling relations for match-ing densities between 2 . n sat (marked with horizontal boxes)and 3 . n sat . Fig. 5 (a)). If an even lower upper bound on ˜Λ . were to be established, the excluded region would becomelarger, increasing the threshold of minimally allowed n m .Compatibility with GW170817 is readily satisfied if the χ EFT calculations (with uncertainties) are assumed validup to 2 . n sat consistent with previous studies [80]. Theevolution of M max with c s, match has been known [56, 81–83], but it was unclear how the uncertainty in the low-density EOS translates to an uncertainty in the de-rived upper bound. As shown in Fig. 6 (b), we findthat for n m = 2 . n sat , the uncertainty in M max rangesfrom ≈ . (cid:12) for c s, match = 0 .
33 (blue-dashed line)to ≈ .
05 M (cid:12) for c s, match = 1 (black-dashed line) withN LO ± σ inputs at low densities.Extrapolating the χ EFT EOS to densities somewhathigher than 2 . n sat will be useful for the subsequent dis-cussion. We find that the ZL parametrization [79] ofNSM matter, which has a single parameter correspond-ing to the symmetry energy coefficient L , is a convenientextrapolation tool. Figure 7 shows M – R curves for theZL EOSs together with a standard crust. For example, L = 45 MeV (65 MeV) successfully tracks N LO, while L = 45 MeV (75 MeV) tracks N LO, for − σ (+ σ ). Using the ZL parametrization, we show the conse-quences of extending the nucleonic EOS to 3 . n sat inFig. 6 (b). We find that the ZL EOSs corresponding to L = 60 MeV and 50 MeV, respectively, smoothly join the M max ( n m ) relations for the N LO +1 σ and N LO-cenEOSs ; the reason is that the masses of stars with cen-tral density n cent = 2 . n sat are similar in both cases. Thederived bounds on n m and c s, match illuminate the impor-tance of including nuclear-matter calculations in the den-sity range 1 − n sat . Standard extrapolations based onnucleonic models, similar to the ZL parametrization, areusually associated with a more gradual profile of c s ( n B )at low-to-intermediate densities, which cannot reconcilethe small radii and/or small tidal deformabilities inferredfor canonical-mass NSs with large maximum masses. Thenecessary rapid change in the sound speed guided by thesimple matching scheme serves to indicate the breakdownof such extrapolations at high densities. A very highNS mass, e.g., (cid:38) .
45 M (cid:12) (2 . (cid:12) ), would be in conflictwith causality and standard extrapolation up to 3 . n sat (2 . n sat ); therefore indicating something unusual in theEOS should should be taking place near this density.This is consistent with the findings of Refs. [45, 46].It is worth mentioning that so far we have intention-ally avoided finite discontinuities in the energy density ε ,which would otherwise introduce an additional parame-ter that characterizes the strength of a sharp first-orderphase transition. In that scenario, the M max bounds willbe shifted downwards due to the softening induced by thephase transition, while GW170817 boundaries may be-come more complicated depending on the possible forma-tion of disconnected branches at intermediate densitieson the M – R diagram [85, 86]. However, given the sys-tematic uncertainties involved in obtaining ˜Λ from grav-itational waveform data, the previously inferred boundsshould still apply [75]. In any case, as discussed in Sec. II,useful information on the minimal radii R min ( M ) can beobtained from matching to the causal EOS with a dis-continuity ∆ ε m specified by M max , and we will elaborateon these lower bounds on R with χ EFT inputs up to n m later in Sec. IV C. Note that M max ( n m ) for the extrapolated EOSs will eventuallybend upwards at sufficiently large n m , which is a generic fea-ture whenever a “standard” nucleonic-like EOS (i.e. graduallyincreasing c s without kinks or discontinuities that naturally ex-tends from low-density e.g. χ EFT calculations) is switched toa linear EOS at some critical density, with or without discon-tinuities in ε (see, e.g., Fig. 5 in Ref. [84]). However, we limitour studies to n m (cid:46) . n sat , as there is little guidance for thevalidity of nucleonic degrees of freedom at higher densities fromtheory. FIG. 8. Panel (a): scaling relations between M max and n max ; panel (b): scaling relations between M max and R M max . Bothrelations, shown as dot-dashed lines, follow from the maximally compact EOS (see Sec. II). The black dashed curves correspondto the presence of a low-density nuclear mantle (crust + N LO EOS) for n B ≤ n m , with fixed sound speeds c s, match = 0 . c s, match = 1 . n B > n m . The grey-shaded region is excluded by GW170817 (˜Λ . ≤
720 and N LO-cen). The solidcolored curves show contours of n m = 1 . , . , . n sat for N LO-cen; dashed colored curves show ± σ uncertainties. For EOSsthat accommodate M max ≥ . (cid:12) , the permitted ranges of n max and R M max are severely restricted. B. Refinements of the limits from the maximallycompact EOSs
The “maximally compact” EOSs [53, 56, 87, 88] forself-bound stars with P = 0, P = 0 for ε ≤ ε in Eq. (1)lead to absolute limits on the central density n max andthe radius R M max of the maximum-mass star. In this sec-tion, we refine these limits including χ EFT uncertaintiesat low densities and tidal deformability constraints.In Fig. 8 (a), the dot-dashed boundary represents theabsolute upper limit on M max as a function of n max ,the highest possible baryon density from the maximallycompact EOSs; see Eq. (7). The slightly lower blackdashed boundary matches the maximally compact EOSto a low-density nuclear EOS at some density n m vary-ing from n sat to about 3 . n sat (from left to right). Therelatively small difference between these two boundariessuggests that effects on the absolute upper bound on n max and M max from the low density EOS is small,and for M max ≥ . (cid:12) , n max should be smaller than5 . − . n sat . This is in good agreement with ≈ n sat ob-tained in Ref. [45]. For n m ≤ . n sat , we employ χ EFTcalculations with uncertainties, and the ZL parametriza-tions (see Fig. 7) are applied for n m between 2 . − . n sat .If the high-density matter is assumed to be much softerwith c s, match = 0 .
33, matching it to the nuclear EOSat different matching densities n m gives rise to the pre-dicted M max – n max relation shown by the lower dashedcurve. The grey-shaded region is ruled out by tidal deformability constraints inferred from GW170817, pro-hibiting small values of n m below 1 . − . n sat (see alsoFig. 5). As a result, c s, match (cid:46) .
33 is incompatible with M max (cid:38) . (cid:12) ; see also Fig. 6. Furthermore, imposing M max ≥ . (cid:12) leads to 5 . < n max /n sat < . n m /n sat = 1 . . . c s, match from 1 to below 0 .
33. In each case, thehighest M max as well as the smallest n max correspondto where they end at the c s, match = 1 upper bound-ary (black dashed line). The N LO ± σ uncertaintyat 2 . n sat translates to ≈ . n sat uncertainty in n max (5 . − . n sat ) if M max = 2 . (cid:12) , and ≈ . n sat uncer-tainty for M max = 2 . (cid:12) . Beyond n m (cid:38) . n sat , ex-trapolation of the χ EFT calculations is needed for whichthe curves would move to the lower-right while remainingunder the c s, match = 1 bound. Using the ZL parametriza-tion to extrapolate up to 3 . n sat (not shown), we obtain n max ≤ . − . n sat and M max ≤ . − .
48 M (cid:12) .As discussed in Sec. II, the maximally compact EOSwith c s = 1 determines the smallest possible radius at agiven mass. Figure 8 (b) displays the absolute bound onthe radius of the maximum mass star, R M max , as well as amore realistic bound taking into account the low-densityEOS below n m . Assuming χ EFT up to n m = 2 . n sat and M max = 2 . (cid:12) , the N LO ± σ uncertainties induce anuncertainty ≈ . R M max = 11 . − .
66 km. For M max = 2 . (cid:12) , an uncertainty ≈ . R M max = 12 . − n m (cid:38) . n sat , M max ≥ . (cid:12) leads to R M max ≥ .
49 km. The tidal deformability constraint inferredfrom GW170817 instead corresponds to limits on theradii of canonical-mass stars. With the simple matchingcondition used here, that constraint simultaneously rulesout too large R M max , e.g., R M max ≤ .
18 km if M max =2 . (cid:12) and R M max ≤ .
79 km if M max = 2 . (cid:12) .On the other hand, introducing a finite discontinuity in ε would decrease R M max and increase n max , but to reachthe same M max necessitates the transition density to besmaller than the matching density n m when there is nodiscontinuity [89]. The overall effect is that larger n max and smaller R M max are possible but must still lie withinthe bounds set by the maximally compact EOSs. C. Bounds on the neutron star radius and tidaldeformability
Here, we address how a specification of the nuclearEOS up to ∼ . n sat including EFT uncertainties, pro-vides robust upper and lower bounds on the radius of aNS of any mass. As we discuss below, these bounds areinsensitive to the details of dense matter physics beyondthe densities accessible by χ EFT.Earlier work has shown that the radii of canonical NSswith masses in the range 1 . − . (cid:12) are most sensitiveto the EOS in the density interval 1 . − . n sat [2], andthis is further quantified in Sec. IV D. As a result, calcu-lations up to (cid:46) . n sat are adequate to place stringentbounds on the NS radius [1, 80, 90]. These studies com-bined information about the EOS at n B (cid:46) . n sat andobservational evidence for NSs with mass M (cid:38) . (cid:12) to derive bounds on the radius of NSs in the observedmass range. Here, we follow a similar procedure but withN LO- χ EFT EOS and its associated EFT truncation-error estimates discussed in Sec. III. The upper boundon the radius is obtained by matching smoothly both theenergy density and pressure to the maximally stiff EOSwith c s = 1 at 2 . n sat . This gives an M – R relation with M max (cid:39) . (cid:12) , as predicted by Eq. (4) with ε (cid:39) ε sat .The minimum radius is determined by introducing afinite discontinuity in the energy density ∆ ε m at n m thatproduces a specified value of M max . Above the density ε m +∆ ε m , the EOS is assumed to be the causal EOS with c s = 1, and the EOS below the density ε m is assumed tobe the original χ EFT EOS. If the pressure at n m was van-ishingly small, this would effectively give the R min ( M )relation for the maximally compact EOS of self-boundstars as described in Sec. II but with ε = ε m + ∆ ε m .With finite pressure at n m based on χ EFT calculations, R min ( M ) will be larger and defines the minimum radiifor normal NSs. The magnitude of the discontinuity∆ ε m at the matching point n m would then be determinedsolely by the maximum mass according to Eq. (4). For M max = 2 . (cid:12) , we find that ∆ ε m ≈ ε nuc ( n m ) (cid:39) ε sat where ε nuc ( n m ) is the energy density in the nuclear- matter phase calculated using χ EFT. To accommodatea maximum mass of 2 . (cid:12) requires a much smaller dis-continuity, ∆ ε m ≈ . ε nuc ( n m ). Furthermore, all thetrajectories within any ± σ band for each value of M max have nearly identical values of ∆ ε m resulting from thefact that P m (cid:28) ε m . Thus the relation between ∆ ε m at2 . n sat and M max is indeed relatively insensitive to thelow-density EOS. R [km]0 . . . . . . . M a ss M [ M (cid:12) ] n m = 2 n sat χ EFT at N LO + causal EOSmax. R min. R ( M max = 2 . (cid:12) )min. R ( M max = 2 . (cid:12) ) FIG. 9. Radius bounds obtained by combining N LO- χ EFTpredictions up to n m = 2 . n sat and maximum-mass infor-mation is shown. The orange bands show the upper boundon the NS radius, while the black and purple bands depictthe lowers bounds corresponding to M max = 2 . (cid:12) and M max = 2 . (cid:12) , respectively. Figure 9 shows the minimum radius and maximum ra-dius bounds. The central values of the minimum radii R min ( M ) for M max = 2 . (cid:12) and M max = 2 . (cid:12) areshown as black and purple solid curves, respectively,while the darker and lighter bands reflect 1 σ and 2 σ un-certainties, respectively. To 2 σ confidence, the minimumradius of a 1 . (cid:12) star ranges from 9 . − . M max is varied from 2 .
00 M (cid:12) to 2 .
93 M (cid:12) ; roughly, the minimumvalue of R . ∝ M / . Similarly, the minimum values of R M max vary from 9 . − . R . > .
68 km and R M max > . R , M max , and the threshold bi-nary mass M thres for prompt collapse of a merger rem-nant. We can therefore provide a more restrictive boundfor R M max since M max is believed to be ≥ . (cid:12) .The maximum radius curves in Fig. 9 are identical tothose presented in Fig. 3 (a). Figure 9 demonstrates howfuture discoveries of NSs with large masses can constrainthe radii of all NSs. Several interesting insights can begleaned from this figure. A striking, albeit expected, fea-ture is the convergence of the upper and lower radiusbounds with increasing M max . This is in accordance withthe facts that the discontinuity ∆ ε m leading to the mini-4 R [km]0 . . . . . . . M a ss M [ M (cid:12) ] n m = 3 n sat χ EFT at N LO + polytrope + causal EOSmax. R min. R ( M max = 2 . (cid:12) ) FIG. 10. Similar to Fig. 9 but obtained using the polytropicextrapolation of the χ EFT EOS up to n m = 3 . n sat . mum radii has to decrease to achieve a higher M max [84]and that the limit ∆ ε m → . (cid:12) NS would be reduced fromabout 3 km when M max = 2 . (cid:12) to about 0 . M max = 2 . (cid:12) . Another feature worth noting is theevolution of the 2 σ lower bound on the NS radius. It in-creases by about 2 km, from 9 . M max = 2 . (cid:12) to11 . M max = 2 . (cid:12) . Comparing the black andpurple bands shows that the radii of heavier neutron starsare even more tightly constrained with increasing M max .Future observational constraints on NS radii in the massrange 1 . − . (cid:12) could be valuable in this regard sinceX-ray and GW observations are best suited to provideradius information at the level of 5% uncertainty in thismass range [92]. Results in Fig. 9 also demonstrate thatan upper bound of about 13 km for R . obtained fromGW170817 is consistent with NSs with M max (cid:39) . (cid:12) .The trends seen in Fig. 9 also have important impli-cations for the EOS of matter at the highest densitiesencountered in the NS inner core. Our results imply that M max > . (cid:12) and/or radii > . (cid:39) . (cid:12) can only be achieved if c s (cid:39) LO- χ EFT calculations. Im-proving the EOS, especially the EFT truncation errorsin the vicinity of n B (cid:39) . n sat , will be critical in extract-ing better constraints on the EOS at higher densities inthe core if future observations favor these large radii ormasses. Supporting c s (cid:39) − n sat requires a formof strongly interacting relativistic matter that poses sig-nificant challenges for dense-matter theory and QCD [93].Encouraged by the apparent convergence of χ EFT cal-culations over the density interval 1 − n sat , it is naturalto ask if a nuclear physics based description of dense mat-ter can be extended to higher density. Extrapolating the EOS from 2 . n sat to 3 . n sat will be model-dependent,even in the absence of phase transitions to non-nucleonicmatter, since we presently do not have reliable calcu-lations at higher densities. We climb this rung of thedensity ladder with some reservation to motivate and ex-plore the impact of future calculations of the EOS in thisdensity interval.We consider a polytropic model where P = κε γ inwhich the parameters κ and γ are determined by fittingto the behavior predicted by χ EFT calculations in thedensity interval 1 . − . n sat to extrapolate the EOS from2 . n sat to 3 . n sat . This choice is somewhat arbitrary andis chosen to approximately capture the key features of thedensity dependence of the EOS predicted by χ EFT. Theresulting radius bounds are shown in Fig. 10. We havechecked that alternate extrapolations using the ZL EOSs,with parameters chosen to suitably match the χ EFT re-sults at 2 . n sat , do not significantly alter our conclusions.A comparison between the results shown in Fig. 9 withthose in Fig. 10 reveals the following insights. First, theincrease in n m does not alter the bounds on R min ( M )(including R M max ), as a function of M max , except that inthe extrapolated case M and M max cannot exceed about2 . (cid:12) . These bounds are therefore particularly robustfor M < . (cid:12) . . . . . . . M [M (cid:12) ]10 T i d a l D e f o r m a b ili t y Λ n m = 2 n sat χ EFT at N LO + causal EOSmax. Λmin. Λ( M max = 2 . (cid:12) )min. Λ( M max = 2 . (cid:12) ) FIG. 11. Bounds on the tidal deformability Λ obtained using χ EFT N LO EOS up to n m = 2 . n sat . As in Fig. 9, theorange bands show the upper bound, while the lower boundscorresponding to M max = 2 . (cid:12) and M max = 2 . (cid:12) areshown by the black and purple bands, respectively. The ver-tical solid line depicts the constraint inferred from GW170817,70 ≤ Λ . ≤ Second, the increase in n m results in more stringentupper bounds on the NS radius for masses in the range1 . − . (cid:12) . For example, the polytropic extrapolationto 3 . n sat predicts R max (1 . (cid:12) ) = 11 . +0 . − . km, whichis to be contrasted with R max (1 . (cid:12) ) = 12 . +0 . − . km ob-tained using n m = 2 . n sat . This reduction has implica-tions for the interpretation of future radius measurements5which aim for an accuracy of better than 5% [92]. If theseobservations favor NSs in this mass range to have radii >
12 km, it would require new mechanisms to rapidlystiffen the EOS below 3 . n sat .Third, if the secondary component in GW190814 wasconfirmed to be a massive NS, new mechanisms wouldalso be implicated at a lower density, since the extrap-olated EOS up to 3 . n sat predicts M max in the range2 . − .
53 M (cid:12) at ± σ .Upper bounds on the tidal deformability Λ can also bederived by smoothly matching the causal EOS to a low-density EOS [94], whereas lower bounds on Λ are deter-mined by matching them to a causal EOS with a discon-tinuity in the energy density determined by M max [75, 86](the same procedure used to determine the minimum ra-dius in Figs. 9 and 10). The bounds for the N LO- χ EFTEOS with n m = 2 . n sat are shown in Fig. 11. The roleof M max is clear from comparison of the M max = 2 . (cid:12) and M max = 2 . (cid:12) cases.In Fig. 12 we show the upper and lower bounds ob-tained by using χ EFT up to 2 . n sat and the polytropicextrapolations from χ EFT up to 3 . n sat . The fact thatuncertainties in the GW170817 constraint of Λ extendalmost precisely between the lower ( M max = 2 . (cid:12) with a large discontinuity ∆ ε m at n m ) and upper bounds( c s, match = 1 without discontinuity) to within 2 σ for both n m = 2 . n sat and n m = 3 . n sat cases is not a coinci-dence. It is a consequence of the fact that for those val-ues of n m , ˜Λ . <
720 is always satisfied for all valuesof c s, match ≤ . . . . . . M [M (cid:12) ]10 T i d a l D e f o r m a b ili t y Λ n m = 3 n sat χ EFT at N LO + polytrope + causal EOSmax. Λmin. Λ( M max = 2 . (cid:12) ) FIG. 12. Upper and lower limits on the tidal deformabilityas a function of mass based on the extrapolated EOS used inFig. 10. The vertical solid line depicts the constraint inferredfrom GW170817, 70 ≤ Λ . ≤ A comparison between the results shown in Fig. 11 andFig. 12 provides quantitative insights into how access tothe EOS at higher density will impact predictions for thetidal deformability Λ, especially for more massive NSs. Itillustrates how constraints on Λ from future GW detec- tions from binaries with massive NSs can provide insightson the evolution of c s in the density interval 2 − n sat .For example, if Λ . (cid:38) χ EFT predictions even in the density interval1 − n sat , and Λ . (cid:38)
50 would be difficult to accom-modate without new mechanisms to significantly stiffenthe EOS in the density interval 2 − n sat . On the otherhand, if Λ . (cid:46) − n sat , a near-causal EOS at higher densities, and M max not significantly larger than 2 M (cid:12) . D. Sensitivity to EOS density ranges
It is apparent that the limits to NS radii and tidal de-formabilities are sensitive to the EOS in the density range1 − n sat , precisely where the restrictions from χ EFT areimportant. This is not surprising given the tight corre-lation between R . and the NSM pressure for 1 − n sat discovered by Ref. [2]. However, up to this point, we haveassumed fixed sound speeds above n m . In this section,we demonstrate that this correlation is insensitive to thedetails of the assumed EOS at all relevant densities; fur-thermore, we quantify this correlation and extend it toinclude the quantities R . and M max .We evaluate these correlations by considering severalparametrization schemes to construct families of high-density NSM EOSs at densities larger than about 0 . n sat ,the assumed core-crust boundary. All configurations areassumed to have a crust modeled with the SLy4 EOS [95].Each EOS is given as a function of n B only and is im-plicitly considered to represent beta-equilibrium matter.The parameters for each parametrization scheme are con-strained to ensure causality, c s ≥
0, a minimum value M max = 2 . (cid:12) , a lower limit to the neutron-matter en-ergy and pressure suggested by the unitary-gas conjec-ture [96] at all supra-nuclear densities, and upper limitsto the NSM energy and pressure at n sat implied by ex-perimental limits of S v = 36 MeV and L = 80 MeV [97].Note that the latter two constraints are broader thanthe NSM- χ EFT ± σ constraints, so that the correla-tions we find are conservatively expressed. Also, for eachparametrization, we have ensured a minimum of 15,000realizations that satisfy our constraints. We quantify acorrelation in terms of the covariance between two quan-tities A and B ,cov( A, B ) = (cid:88) i,j ( A i − ¯ A )( B j − ¯ B ) σ A σ B . (25)The σ ’s represent standard deviations. We take A = P ( n B ) and B = R . , R . , or M max . Here, j ranges overall realizations of a given parameterized EOS and i overall values of n B smaller than the central density of therelevant configuration for B .Figure 13 shows the correlations between the pres-sure P ( n B ) and R . , R . and M max as functions ofthe baryon number density n B for a variety of NSM6 n B (fm ) c o v ( P ( n B ) , R . ) n-EXPk-EXPSpectral4QuarkyonicPP3+1RMFAverage n B (fm ) c o v ( P ( n B ) , R . ) n-EXPk-EXPSpectral4QuarkyonicPP3+1RMFAverage n B (fm ) c o v ( P ( n B ) , M m a x ) n-EXPk-EXPSpectral4QuarkyonicPP3+1RMFAverage FIG. 13. Correlations among P ( n B ), R . , R . and M max for 6 EOS parametrizations (see text for details). “Average”refers to the mean of all models. Blue histograms show thesummed distributions of the central densities of the relevantstars. parametrizations in common use. The parametrizations“n-EXP” and “k-EXP” are three-parameter Taylor ex-pansions of the NSM energy in terms of n B and n / [96],respectively. “n-EXP” is commonly used to model thenuclear energy around saturation; we take a Taylor ex-pansion up to the fourth-order term [( n B − n sat ) /n sat ] . Two of the coefficients are set to match the crust EOS,leaving three free parameters. “k-EXP” contains a ki-netic term ∝ ( n B /n sat ) / and a higher-order term up to( n B /n sat ) / . It also has three free parameters after usingtwo coefficients to match the crust EOS. ‘ ‘Spectral4” isthe four-parameter spectral decomposition method [98–100]. “Quarkyonic” has two parameters, Λ and κ , speci-fying the quarkyonic momentum shell thickness and thetransition density, and one parameter (effectively con-trolling L ) for the nucleon potential [93]. “PP3+1” is afour-parameter piecewise-polytrope with three segmentsappended to the crust [101]. The density n separatingthe first two segments is a parameter, while n and n arechosen to scale as n = 2 n and n = 2 n . The corre-sponding bounding pressures P , P , and P are the otherthree free parameters . “RMF” is a relativistic mean fieldmodel based on the FSU2 EOS [102] and contains σ , ω ,and ρ meson exchanges. It has seven coupling constants,of which three are fixed by saturation properties of SNM;the remaining four free parameters can be mapped to S v , L , the effective nucleon mass at the saturation density, M ∗ , and the ω self-interaction coupling ζ .The covariance parameter cov( P ( n B ) , R . ) peaksaround n B = 1 . +1 . − . n sat , whereas cov( P ( n B ) , R . )and cov( P ( n B ) , M max ) peak around n B = 2 . +2 . − . n sat and n B = 3 . +2 . − . n sat , respectively. The uncertaintiescorrespond to 50% of the peak covariance. Figure 13 alsoquantifies the extent to which the central baryon densi-ties, and the width of their distributions, increase withthe NS mass. Notably, the central baryon number densi-ties peak at about 30% higher density than do the peakcovariance in all three cases, but the widths of the centraldensity distributions rapidly increase with NS mass.The correlation between the pressure P ( n B ) and R . isstrongest between n sat and 3 . n sat , as expected, and thatbetween the pressure and R . is strongest at about 40%higher densities. Significantly, these results appear to berelatively insensitive to the details of the parametriza-tions. The standard deviations of both cov( P ( n B ) , R . )and cov( P ( n B ) , R . ) for the six parametrizations aresmall, being σ cov ,R < . σ cov ,R < .
05 near the covariance peaks. The bottom line is theseresults demonstrate, at present, that χ EFT greatly con-strains R . and, to a slightly lesser degree, R . . Thesituation is somewhat different for M max , where pres-sures at densities between 2 . n sat and 6 . n sat dominate.In addition, the standard deviation of cov( P ( n B ) , M max )among the six parametrizations are somewhat larger, be-ing σ cov ,M max < .
25 at all densities and σ cov ,M max < . M max resultsare more model-dependent, and the significant densitieslikely lie above the validity range for χ EFT. However,further refinement of EFT techniques at high densities The additional parameter n greatly increases the flexibility ofPP3+1 compared to the three-parameter ( P , P , P ) set PP3often employed [101]. n sat . V. DISCUSSION AND OUTLOOKA. Current and future constraints
To shed light on the properties of dense matter, the ob-servational constraints used in this work are taken from(i) a handful of well measured NS masses from radio ob-servations [8–12], (ii) the chirp and combined masses aswell as bounds on tidal deformabilities of NSs deducedfrom GW detections in the binary NS-NS merger eventGW170817 [5–7], and (iii) radius estimates from NICERfor a NS of mass (cid:39) . (cid:12) [103, 104]. An upper bound of M max (cid:46) . (cid:12) on the maximum gravitational mass of acold, spherical NS was inferred from several studies usingEM and GW data from GW170817 [3, 4, 105–107], butan upper bound on M max itself does not provide furtherlimits on the sound speed or bounds to NS radii since theEOS could suddenly soften above n m .The NICER M – R constraints on J0030+0451, namely, R = 13 . +1 . − . km, M = 1 . +0 . − . M (cid:12) [104] and R =12 . +1 . − . km, M = 1 . +0 . − . M (cid:12) [103], and some EMobservations of GW170817 [108–110] favor larger radiithan indicated by GW observations from GW170817,10 −
13 km [6, 7], but the degree of tension is slight. Jointanalyses of these data yield tighter but still consistentconstraints on the typical NS radius ∼ . . +1 . − . km to 68.3% confidence.It is fortunate that NICER targets also include severalpulsars for which the masses are independently measuredto high precision, e.g., PSR J1614-2230 (cid:39) .
91 M (cid:12) andPSR J0740+6620 (cid:39) .
14 M (cid:12) , and PSR J0437-4715 [116]with mass ≈ .
44 M (cid:12) . The possibility to measure radiiof both intermediate as well as very massive NSs opensup the possibility to contrast the radii of ∼ . (cid:12) stars, R . , and more typical ∼ . (cid:12) stars, R . , to furtherconstrain the EOSs [89].We show in Fig. 14 the difference ∆ R = R . − R . for stars with the N LO EOS up to 2 . n sat , ZL EOS ex-trapolations up to a range of matching densities n m =2 . − . n sat , and various linearly matched EOSs withdifferent c s, match at higher densities. The ZL extrapo-lation with L = 50 MeV indicates that roughly above n m (cid:38) . n sat , all values of c s, match lead to R . ≤ R . .The boundary between positive and negative ∆ R shifts abit when using the slightly stiffer ZL extrapolation with L = 60 MeV: in this case n m (cid:38) . n sat will guarantee R . ≤ R . ; note that 2 . n sat is already the central den-sity of a 1 . (cid:12) star. We also checked radii differencesbetween 2 . (cid:12) and 1 . (cid:12) stars, ∆ R (cid:48) = R . − R . ,and found that ∆ R (cid:48) is generally less than ∆ R , with thelargest decreases of a few tenths of a km occurring for the FIG. 14. Radius differences ∆ R = R . − R . using the ZLextrapolations with L = 50 MeV and L = 60 MeV joined con-tinuously to linear EOSs at n m between 2 . n sat and 3 . n sat . smaller values of c s, match . For c s, match (cid:38) .
7, there arenegligible differences. ∆ R or ∆ R (cid:48) being negative is typ-ical when extrapolations to even higher densities are ap-plied, or if there is additional softening in the EOS beforereaching the central density of the maximum-mass star.Should observations suggest R . > R . or R . > R . ,standard extrapolations such as ZL-models predict someunusual stiffening should occur below (cid:46) . − . n sat .Furthermore, if ∆ R turns out to be greater than 0 . n m (cid:46) . n sat , which suggests a very high M max and lesscompatibility with radius constraints from GW170817;see Fig. 6 (b). However, NICER observations may notachieve the needed O (0 . ∼ . n sat correspond to0 . − . (cid:12) within 1 σ uncertainties of χ EFT calcula-tions (Fig. 9), it will be greatly helpful if radii of verylow-mass NSs ∼ . (cid:12) can be obtained through X-rayobservations, or tidal deformability measurements of bi-nary systems with very low chirp masses.From a different perspective, more accurate experi-mental determinations of S v and L at n sat from e.g.,PREX, CREX, and FRIB/MSU, will be important totest χ EFT predictions of properties of neutron-rich mat-ter. At the present time, S v and L are believed to beunderstood to the 10% and 40% levels, respectively [97].For n B > n sat , constraints from the analyses of the col-lective flow of matter in HICs could be informative.The best available information for the present comesfrom the analysis of HICs of Au nuclei using Boltzmann-type kinetic equations. The elliptic and sideways flow ob-servables from these collisions are sensitive to the mean-8field potential and to in-medium NN collisions at cen-tral densities of 2 − n sat , and suggest SNM pressures of7 . − to 14 MeV fm − at 2 . n sat [117]. In com-parison, N LO calculations for SNM predict somewhatlarger pressures of 10 . − to 18 . − at2 . n sat [42], which are, nevertheless, consistent withintheir stated 1 σ uncertainties. However, the predictionsfrom HICs involve model-dependent assumptions con-cerning the density- and momentum-dependencies of theassumed nuclear interactions, which have not been sys-tematically explored; see Ref. [118] and references thereinfor the relevance of single-particle potentials in HICs. Inaddition to these uncertainties, HICs probe nearly sym-metric matter, and to apply their observables to NSMrequires an additional extrapolation involving the sym-metry energy at supra-nuclear densities.To improve the current status, heavy-ion facilitiesacross the world, such as RHIC, FAIR, NICA, J-PARC,and HIAF, have launched programs to map out the QCDphase diagram of strongly interacting matter. The studyof more neutron-rich matter in HICs, together with im-proved, systematic, modeling would be very valuable fordense-matter physics, not only for cold neutron stars,but also for understanding mergers involving NSs. Asthe analyses of HIC data have largely been done withnucleonic degrees of freedom, it would be also interestingand desirable to extend such analyses to include quarkdegrees of freedom and their subsequent hadronizationas in RHIC and CERN experiments at higher energies. B. . (cid:12) neutron stars and the nature of thecomponents of GW190425 and GW190814 The ranges of 1 . − . (cid:12) in GW190425 [119] and2 . +0 . − . M (cid:12) in GW190814 [120] for one of the com-ponents in these merger events have raised the possi-bility that those compact objects could be NSs as op-posed to being low-mass BHs. The data from GW190425was inconclusive concerning the nature of the inspirallingbinary [119], but some works favored the scenario inwhich the more massive component is a BH instead ofa very heavy NS [121]. If it is a priori assumed that M max (cid:46) . (cid:12) , a possibility motivated by EM and GWdata from GW170817, the interpretation that it was aBNS merger instead statistically favors masses of approx-imately 1 . ± . (cid:12) and 1 . ± . (cid:12) , while a neutron-star-black-hole (NSBH) merger interpretation favors a1 . ± . (cid:12) NS and a 2 . ± . (cid:12) BH [121]. Whileboth scenarios are statistically equally likely, the fact thatthe BNS masses are incompatible with those of observedgalactic BNS systems, while the NS mass in the NSBHscenario is compatible, seems to favor the NSBH interpre-tation. However, in either scenario according to this anal-ysis, GW190425 would likely not contain a NS > . (cid:12) .In the case of GW190814, there is no additional informa-tion, aside from one’s assumption about M max , to decideif the primary is a high-mass NS or a low-mass BH. How- ever, statistical analyses suggest that the probability ofit’s secondary being a NS is very low [47, 48, 120]. If ei-ther GW190425 or GW190814 contains a ∼ . − . (cid:12) NS, questions to address are:
What is the physical stateof dense matter that could support such a heavy NS, andwhat radius constraints would follow?
The scenario that GW190814’s secondary componentwas an approximately 2 . (cid:12) NS does not itself vio-late theoretical limits from causality and the GW170817constraint that ˜Λ <
720 for M = 1 .
186 M (cid:12) , but chal-lenges remain finding physical mechanisms that can con-nect very stiff high-density matter with the relativelysoft nuclear matter at (cid:46) . n sat predicted from mod-ern χ EFT calculations. As Fig. 6 shows, the conformallimit c s ≤ / M max (cid:38) . (cid:12) and fromGW170817’s tidal deformability constraint. Standard ex-trapolations that assume gradually increasing c s profilesare unlikely to be compatible with M max ≥ . (cid:12) [120].In particular, the requirement that c s remains above ∼ . (cid:38) . n sat is hard to ex-plain. Extrapolations of non-relativistic potential mod-els generally result in steadily increasing sound speedswith density, and it becomes problematic to prevent themfrom becoming acausal within NSs. At densities relevantto the center of very massive NSs, it is reasonable toexpect the emergence of exotic degrees of freedom. Asharp first-order transition to stiff quark matter at someintermediate density is capable of reconciling small radiiand high masses (cid:38) . (cid:12) (see examples of R min ( M )in Sec. IV C). With an increasing lower bound on M max and/or smaller assumed values of c s at high densities,the transition threshold has to be pushed downward ap-proaching 1 . − . n sat (similar to the results shown inFig. 6 but involving a discontinuity ∆ ε that further de-creases M max and favors lower values of n m [89]).For most microscopic quark-matter models, for ex-ample the original MIT bag model [122], the originalNambu-Jona–Lasinio (NJL) model [123], and their varia-tions, perturbative QCD matter [124], and quartic poly-nomial parametrizations [125], the speed of sound turnsout to be weakly density-dependent. To be consistentwith massive pulsars ∼ (cid:12) , strong repulsive inter-actions that stiffen the quark EOS, possibly reaching c s ≥ .
4, have been implemented [126–128]. The max-imally achievable c s is model-dependent, and requiring c s (cid:38) . c s by restrict-ing the nucleonic momentum phase space when quarksappear, and in some cases are capable of simultaneouslyreaching > . (cid:12) and satisfying the GW170817 con-straint ˜Λ . < c s peak behavior, cannot jointly sat-isfy these conditions, reaching at most M max (cid:39) . (cid:12) .However, we find that other versions [93, 129], in whichthe quark abundances grow more slowly and that can re-tain large abundances of nucleons at high density, cansimultaneously achieve these conditions.Using extrapolation functions in terms of c s and µ ,Annala et al. [131] found that the risk of hadronic EOSsviolating causality at high-enough densities ( (cid:38) . n sat )to achieve high masses is remedied if a transition to per-turbative QCD-like (soft c s ≈ /
3) quark matter occursat high densities. However, considering that current cal-culations in perturbative QCD itself are only valid atdensities n B (cid:38) n sat , interpolations down to NS densi-ties are problematic. The main feature of such a transi-tion can be reproduced by simply requiring c s → / n B (cid:38) n sat , but at intermediate densities the conformallimit c s ≤ / M – R relation or tidal deformabilitiesdue to the masquerade problem [125]. C. Theoretical aspects of the dense matter EOS
In this study, we have used microscopic calculations ofthe EOS of PNM and SNM with χ EFT NN and 3N inter-actions up to N LO to construct the EOS of NSM. N LOis currently the highest order in the χ EFT expansion atwhich all two- and many-body interactions have been de-rived. This study was possible because recent advancesin MBPT [35] have enabled significantly improved PNMpredictions at this order, and, for the first time, providedorder-by-order calculations of the SNM EOS with NNand 3N interactions up to N LO [35, 42, 44]. The un-derlying assumptions are that nucleons are the relevantdegrees of freedom in this density range, and that χ EFTprovides a systematic expansion for the EOS to calcu-late EFT truncation errors. While these are reasonableassumptions, several questions remain and more work isneeded to address them.Experimental validation of χ EFT predictions for theEOS of bulk matter relies on comparisons to the empiri-cal saturation point, and constraints on the nuclear sym-metry energy and its derivative with respect to density at n sat . While the MBPT calculations used in this work arewell within the joint experimental constraint in the S v – L plane [42], the Λ = 500 MeV Hamiltonians—as discussedin Ref. [35]—actually do not saturate inside the empiricalrange for the saturation point, n sat = 0 . ± .
007 fm − with ( E/A ) sat = − . ± .
57 MeV. Note, however,that this empirical range was obtained in Refs. [35, 67]from a set of energy density functionals, and thus onlyhas limited statistical meaning. The predicted 2 σ confi- dence ellipses for the nuclear saturation point at N LOand N LO are shown in Fig. 9 of Ref. [43].Nuclear saturation in SNM emerges from a delicatecancellation of different contributions in the nuclearHamiltonian, in contrast to the properties of neutron-richNSM EOS. This cancellation is sensitive to the short-and intermediate-range 3N interactions at N LO that donot contribute to the PNM EOS; e.g., the 3N contactinteraction ( ∝ c E ) is Pauli-blocked in PNM [63]. To-gether with the fact that the proton fraction is small, thismeans that the nuclear saturation properties are of rela-tively minor importance for constructing the NSM EOS.Nonetheless, a better understanding of nuclear saturationproperties may help identify and quantify systematic un-certainties in the nuclear interactions. This might alsolead to a better understanding of the link between (satu-ration) properties of infinite matter and medium-mass toheavy nuclei [132, 133] to explain why χ EFT potentialsgenerally tend to underestimate charge radii [134–136].In this context, it is worth noting that systematic EFTcalculations of the EOS of NSM, which is characterizedby a small proton fraction, would obviate the need to relyon the quadratic expansion Eq. (20) (see, e.g., Ref. [31] inwhich the energy of adding a proton to PNM was calcu-lated). When such calculations become available one cangauge the extent to which the EOS of NSM is correlatedwith the empirical properties of SNM.Developing and applying improved order-by-order χ EFT NN and 3N potentials up to N LO within differentregularization schemes to finite nuclei and infinite mat-ter is an important task for future work [132, 133, 136–138]. While this work has already accounted for cor-related EFT truncation errors (as recently derived bythe BUQEYE collaboration), uncertainties in the low-energy couplings at a given χ EFT order, e.g., from fit-ting the nuclear interactions, have not been considered.Bayesian parameter estimation can be used to obtainlow-energy couplings with statistically robust uncertain-ties [139–141]. A full Bayesian analysis of the EOS wouldthen enable a consistent quantification and propagationof theoretical uncertainties, from those in the fits of thenuclear interactions to EFT truncation errors, via MonteCarlo sampling and the Jupyter notebooks provided bythe BUQEYE collaboration [41].Finally, we note that theoretical studies to access theEOS of matter at n B (cid:38) n sat can have a significant im-pact on NS properties, especially on the correlation be-tween M max and the NS radii. Detailed studies of EFTtruncation errors at these higher densities, especially in χ EFT formulations with explicit ∆ intermediate stateswould be valuable. Models that include additional de-grees of freedom such as pions, hyperons, and quarkswhile still being able to accommodate massive NSs canprovide new insights. Further work is needed to im-prove these models. Together with advances in nuclear-matter calculations from χ EFT at low densities (see, e.g.,Refs. [35, 36]), these will enable astrophysical applica-0tions based on microscopic calculations over a wide rangein density and proton fraction.
D. Comparison with other works
As noted earlier, the uncertain nature of the less com-pact object in GW190814 with mass (cid:39) . (cid:12) has piquedthe interest of the dense-matter and nuclear-physics com-munities. Below we briefly discuss how our study differsfrom or complements the findings of several other recentarticles [45–52] that have addressed the implications ofthe possible existence of NSs with such high masses.Several of these articles, including Refs. [45, 47, 49],have relied on nuclear physics based EOSs to describematter in the crust and outer core to show that the exis-tence of a 2 . (cid:12) NS would require c s ≥ . c s in the inner core. Mostnotably, Ref. [51] derives strict upper bounds on the max-imum mass of NSs that depend only on bulk propertiesof NSs, such as the radii and the tidal deformabilities tofind that a NS in GW190814 would not be inconsistentwith present astronomical constraints if c s is large in theinner core. Our finding suggests that a 2 . (cid:12) NS wouldrequire c s ≥ . − . LO- χ EFT EOS that allows us to properly incorporate EFTtruncation errors at n B ≤ . n sat .Lim et al. [46] combine nuclear models valid in thevicinity of normal nuclear densities and a maximally stiffEOS at higher density to show that 2 . − . (cid:12) NScan exist without strongly affecting the properties suchas radius, tidal deformability, and moment of inertia ofcanonical NSs with mass ∼ . (cid:12) . They argue thatproperties of NSs with masses ∼ (cid:12) such as R ∼ . would be significantly different depending on whether thesecondary component of GW190814 was a black hole or aNS. Our results support these findings, but go beyond bydelineating how the lower and upper bounds on the radiiof NSs in the mass range 1 . − (cid:12) would be constrainedif future observations were to confirm the existence of NSswith masses (cid:39) . − . (cid:12) .Using FSU-type relativistic mean field-theoretical(RMFT) models, Fattoyev et al. [50] found that the rapidincrease in pressure with density required to support a2 . (cid:12) NS, while barely accommodating the deformabil-ity constraint from the first analysis of GW170817 datathat indicates Λ . ≤
800 [5] but not the updated bounds70 ≤ Λ . ≤
580 [7] (see Ref. [142] for a similar study), isinconsistent with energy density functionals tuned to re-produce properties of nuclei and flow data from HICs.Note that Fattoyev et al. [50] only applied Λ . con-straint without a comparison of the binary tidal deforma-bility ˜Λ. We have confirmed that FSU-like RMFT in- teractions cannot accommodate both ˜Λ . ≤
720 and M max ≥ .
54 M (cid:12) [115].Other recent works studied hyperonic matter in theEOS and/or rapid rotations that stabilize more massivestars than non-rotating configurations, which may or maynot be consistent with GW190814 [143–146]; we do notconsider these effects in the present paper.
VI. CONCLUSIONS
In this paper, we have highlighted the important role M max plays in determining bounds on the radii of neutronstars. In the most extreme case, in which only causalityis assumed with the EOS ε = ε + P , absolute upperbounds on M max (cid:39) .
09 M (cid:12) and R M max (cid:39) . ε > ε sat (Eq. (4) and Eq. (5)). Firm lowerbounds on R min ( M ) and R M max that scale with M max canalso be established. For the case that M max = 2 . (cid:12) , R min (1 . (cid:12) ) = 8 . R M max = 8 . M max ≥ . (cid:12) , we find R min (1 . (cid:12) ) > .
75 kmand R min (2 . (cid:12) ) > . FIG. 15. The same plot as Fig. 6 but showing M max contourson the ( c s, match , n m ) plane. Note that the GW170817 bound-aries excluding n m (cid:46) . − . n sat will be shifted downwardsif there is first-order transition at such low densities. If instead an upper limit c s < ε = ε + P/c s , then R min ( M ) and M max depend sen-sitively on c s and decrease with it. For the case c s =1 / ε = ε sat , for example, M max = 2 .
48 M (cid:12) , R min (1 . (cid:12) ) = 12 . R M max = 13 . LO in χ EFT. For a given n B , the NSM EOS always has a lower1 m /n sat R . ( k m ) (a) R max R M max = 2.61.4, min R M max = 2.31.4, min R M max = 2.01.4, min m /n sat R . ( k m ) (b) R max R M max = 2.62.0, min R M max = 2.32.0, min R M max = 2.02.0, min R max R M max = 2.61.4, min R M max = 2.31.4, min R M max = 2.01.4, min R max R M max = 2.62.0, min R M max = 2.32.0, min R M max = 2.02.0, min FIG. 16. Similar to Figs. 9 and 10, but displaying the minimum and maximum radii of 1 . (cid:12) (panel (a)) and 2 . (cid:12) (panel(b)) stars as a function of the matching density n m = 1 . − . n sat . Additionally, R min contours and uncertainty bands for thecase M max = 2 . (cid:12) are shown. ε than the PNM EOS. The pressure of NSM is less thanPNM at the same n B , typically by < − , exceptfor n B (cid:38) . n sat when it becomes greater (Fig. 2 (b)).The proton fraction below 2 . n sat never exceeds the crit-ical minimum value required for the direct URCA processof enhanced neutrino emission [147, 148].The existence of the NS crust together with a nucleonicEOS below a matching density n m establishes R max ( M ).Extremes are again found by assuming c s, match = 1 fordensities above n m , for which the EOS is now ε = ε m + P − P m . Assuming ε m = ε sat , and that P m is given by χ EFT-N LO, the upper bounds are R . , max ≈ . R . , max ≈ . n m = n sat ),which are nearly identical to the case shown in Fig. 1 witha slightly different value of P m at ε m = ε sat . These valuesare not in tension with observations, and with increasing n m , the corresponding upper bounds on R . and R . decrease. For the same ε m or n m , M max is not sensitiveto the value of P m or the nucleonic EOS between thecrust and n m , and is close to that of the case P = 0(self-bound stars) for the causal EOS; see also Fig. 3.Using the N LO EOS between the crust and n m ,we determined how M max depends on n m and c s, match (Fig. 15); we find that M max ≥ . (cid:12) requires c s, match > .
35 (i.e., the EOS violates the conformal limit c s ≤ / n m = n sat , and c s, match > . n m = 2 . n sat . Theconformal limit is also violated for n m > . n sat , evenif M max is as low as 2 . (cid:12) . If M max > .
45 M (cid:12) , n m must not exceed 3 . n sat no matter what the value of c s, match is. The calibrated uncertainties in χ EFT-N LOlead to relatively small uncertainties, less than 0 . (cid:12) ,in M max ( n m , c s, match ). Additionally, satisfying the GW170817 tidal deforma-bility constraint ˜Λ . <
720 and imposing M max > . (cid:12) requires n m > . n sat and c s > .
35. This limitis not very sensitive to M max . Even if M max > . (cid:12) , itis required that 1 . < n m /n sat < . c s, match > . . (cid:12) star evidently requires a sig-nificant change from normal hadronic EOSs to a muchstiffer EOS between 1 . n sat and 2 . n sat . In the presenceof a discontinuity in ε , the lower bound n m (cid:38) . n sat can decrease, whereas the upper bound n m (cid:46) . n sat remains unaffected as it is imposed by causality.For stars with a normal crust, refined upper limits to R M max can be found using the GW10817 constraint andan assumed value for M max , while lower limits follow fromthe causal EOS: 9 km < R M max < . M max ≤ . (cid:12) and 11.3 km < R M max < . M max ≤ . (cid:12) .The minimum radii R min ( M ) established for the self-bound configurations are relaxed in the case of stars withnormal nuclear-matter crusts, by using N LO up to n m and applying a phase transition into the causal EOS( c s, match = 1) involving a finite discontinuity ∆ ε m at n m .These minimum radius bounds will be highly sensitiveto the assumed M max and steadily increase with it; see,e.g., Fig. 9. For the same value of n m , the correspond-ing maximum radii R max ( M ) follow by attaching to thecausal EOS but without a density discontinuity, whichgives rise to a unique value of M max at a given n m andtherefore the maximum radius bounds do not rely on theassumptions about M max .The results shown in Figs. 9 and 10 are summarizedin Fig. 16 for the specific cases of R . and R . , with abroader range of n m explored between 1 . − . n sat . If2 m /n sat R . ( k m ) (a) R c s = 1/31.4, max R M max = 2.0, c s = 1/31.4, min m /n sat R . ( k m ) (b) R c s = 1/32.0, max R M max = 2.0, c s = 1/32.0, min R c s = 1/31.4, max R M max = 2.0, c s = 1/31.4, min R c s = 1/32.0, max R M max = 2.0, c s = 1/32.0, min FIG. 17. The maximum (orange) and minimum (black) bounds on R . and R . assuming c s ≤ / n B = n m ; χ EFT-N LO uncertainties are indicated (darker bands: ± σ ; lighter bands: ± σ ). The bands merge and terminate at critical matchingdensities above which M max < . (cid:12) . χ EFT-N LO is assumed valid up to 2 . n sat , the upperand lower bounds on NS radii are substantially tightenedin comparison with using χ EFT-N LO only up to n sat :for example if M max ≥ . (cid:12) then R . must lie be-tween 12 . +0 . − . km and 9 . +0 . − . km at the 1 σ level, whichis consistent with earlier studies in Ref. [80].Figure 16 also conveniently illustrates the dramatic ef-fect of increasing the lower bound on M max for the al-lowed ranges between R min ( M ) and R max ( M ), whichimproves (shrinks) the R . ( R . ) bounds by an aver-age 3 km / M (cid:12) (5 km / M (cid:12) ); these limits could be furtherrestricted by forthcoming observations. We note that R . or R . < . M max > . (cid:12) (assuming n m = 2 . n sat ). In addition, iffuture measurements from different sources and messen-gers, e.g., X-ray data from QLMXBs or PREs (or GWdetections of mergers by LIGO) vs NICER targets, wereto exhibit discrepancies in the radius inference close toor larger than the gaps between the minimum and maxi-mum bands shown on this figure, then these are hints ofa large energy-density discontinuity ∆ ε in the EOS (ac-companied with high-density stiff matter) occurring at n B (cid:46) n m [89].There has been speculation that the speed of soundin QCD at finite baryon density may be bounded by theconformal limit which requires c s < / N c ) gauge theories for which a holographic or gravitydual exist. In these theories the speed of sound can becalculated at finite baryon density in the large- N c limitusing classical supergravity methods in a curved space-time [150], and for a large class of such theories (for ex- ceptions, see Refs. [151, 152]) c s < / c s < / c s (cid:39) .
2, then decreasesacross hadron-quark cross-over region, corresponding totemperatures in the range 100 −
200 MeV, and eventuallyincreases again to reach its asymptotic value of c s (cid:39) / T (cid:39)
500 MeV [154].Motivated by the discussion above, we briefly commenton the astrophysical implications of the conjecture that c s < / c s < / M max > . (cid:12) while still allowing for a soft EOS at intermediate den-sity needed to ensure that R . <
13 km. This is alsoevident from Fig. 15 which shows that when c s < / σ level, to simultaneously satisfythe tidal deformability constraint from GW170817 and M max > . (cid:12) if χ EFT-N LO is valid beyond 1 . n sat .Figure 17 shows how the bounds on the radius areinfluenced when c s < / R . with n m is strikingand implies that if c s < / χ EFT-N LO is validup to 1 . n sat , then R . must lie between 12 . +0 . − . kmand 13 . +0 . − . km at the 1 σ level. Further, requiring that M max > . (cid:12) excludes a significant fraction of the χ EFT-N LO predicted range for the pressure for densi-ties between 1 . − . n sat . A tiny sliver of high pressureclose to the edge of the 2 σ boundary remains, and im-3 m /n sat R . ( k m ) (a) R c s = 1/21.4, max R M max = 2.0, c s = 1/21.4, min m /n sat R . ( k m ) (b) R c s = 1/22.0, max R M max = 2.0, c s = 1/22.0, min R c s = 1/21.4, max R M max = 2.0, c s = 1/21.4, min R c s = 1/22.0, max R M max = 2.0, c s = 1/22.0, min FIG. 18. Similar to Fig. 17, but assuming c s ≤ / n B = n m . plies that R . = 13 . ± . R . areshown in the right panel.In Fig. 18 we show the maximum and minimum boundson R . and R . obtained by imposing an intermedi-ate limit of c s ≤ /
2. In this case for n m = 2 . n sat ,we find that 11 . +0 . − . km < R . < . +0 . − . km and M max < . ± .
04 M (cid:12) (Fig. 6 (b)), to 1 σ confidence.The corollary to this implies that measurements of R ∼ . that are smaller than 11 . c s ≥ / . n sat , or that n m < . n sat .This is particularly interesting because a recent analysisof the tidal deformability constraints from GW170817 inRef. [78] suggests 11 . +0 . − . km (90% credible interval).We also quantified the sensitivity of the key observables R . , R . , and M max to the pressure as a function ofdensity P ( n B ), and found that the highest correlations,i.e., the most sensitive regions, involve the density ranges1 . − . n sat , 1 . − . n sat , and 2 . − . n sat , respectively.We showed that positive values of ∆ R = R . − R . ,potentially possible with NICER, would indicate lowmatching densities (cid:46) . − . n sat and relatively largevalues of c s, match (cid:38) . − .
6, which would also implylarge values of M max . In the absence of a dramatic stiff-ening of the EOS near 2 . n sat , the expectation is that∆ R <
0. This is usually the case if extrapolations basedon nucleonic-like models are used up to even higher den-sities and/or there is extra softening below M max .The merger events GW190425 and GW190814 are eachconsistent with at least one component (cid:38) . (cid:12) whichcould be either a massive NS or a low-mass BH, althoughGW190425 could instead involve two ∼ . (cid:12) NSs.Should either system contain a NS with M (cid:38) . (cid:12) , theimplications would be that the conformal limit c s ≤ / n m is likely larger than n sat ); if n m > . n sat (2 . n sat ), the average c s above n m should be > . . n m (cid:38) . n sat (could be low-ered if there is sudden softening in the EOS induced by astrong first-order transition) and c s, match (cid:38) . M max ,the radii of NSs, and the role of the nucleonic EOS fordensities beyond n sat . We have illustrated that χ EFTcalculations provide an EOS for NSM up to ∼ . n sat with quantifiable EFT uncertainties [42, 43], and more-over, these uncertainties are small enough to have rela-tively minor influence on our major conclusions. Nev-ertheless, our results reveal that a reliable extension ofthe EOS to densities moderately in excess of 2 . n sat isgreatly desired, and could provide important additionalconstraints on the relations among M max , R min ( M ), and R max ( M ). These theoretical relationships would be soonconfronted with X-ray, radio, and GW observations. ACKNOWLEDGMENTS
We thank R. J. Furnstahl, B.-A. Li, J. A. Melendez,and D. R. Phillips for useful discussions, and the Net-work for Neutrinos, Nuclear Astrophysics, and Symme-tries (N3AS) for encouragement and support. C.D. ac-4knowledges support by the Alexander von HumboldtFoundation through a Feodor-Lynen Fellowship and theU.S. Department of Energy, the Office of Science, theOffice of Nuclear Physics, and SciDAC under awardsDE-SC00046548 and DE-AC02-05CH11231. S.H. is sup-ported by the National Science Foundation, Grant PHY-1630782, and the Heising-Simons Foundation, Grant 2017-228. J.M.L. and T.Z. acknowledge support by theU.S. DOE under Grant No. DE-FG02-87ER40317 andby NASA’s NICER mission with Grant 80NSSC17K0554.M.P.’s research was supported by the Department of En-ergy, Grant No. DE-FG02-93ER40756. The work of S.R.was supported by the U.S. DOE under Grant No. DE-FG02-00ER41132. [1] S. Gandolfi, J. Carlson, and S. Reddy, Phys. Rev. C , 032801 (2012), arXiv:1101.1921 [nucl-th].[2] J. Lattimer and M. Prakash, Astrophys. J. , 426(2001), arXiv:astro-ph/0002232.[3] B. Margalit and B. D. Metzger, Astrophys. J. Lett. ,L19 (2017), arXiv:1710.05938 [astro-ph.HE].[4] M. Shibata, E. Zhou, K. Kiuchi, and S. Fujibayashi,Phys. Rev. D , 023015 (2019), arXiv:1905.03656[astro-ph.HE].[5] B. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev.Lett. , 161101 (2017), arXiv:1710.05832 [gr-qc].[6] B. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X , 011001 (2019), arXiv:1805.11579 [gr-qc].[7] B. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev.Lett. , 161101 (2018), arXiv:1805.11581 [gr-qc].[8] P. Demorest, T. Pennucci, S. Ransom, M. Roberts, andJ. Hessels, Nature , 1081 (2010), arXiv:1010.5788[astro-ph.HE].[9] J. Antoniadis et al. , Science , 6131 (2013),arXiv:1304.6875 [astro-ph.HE].[10] E. Fonseca et al. , Astrophys. J. , 167 (2016),arXiv:1603.00545 [astro-ph.HE].[11] Z. Arzoumanian et al. (NANOGrav), Astrophys. J.Suppl. , 37 (2018), arXiv:1801.01837 [astro-ph.HE].[12] H. T. Cromartie et al. , Nature Astron. , 72 (2019),arXiv:1904.06759 [astro-ph.HE].[13] F. ¨Ozel and P. Freire, Ann. Rev. Astron. Astrophys. ,401 (2016), arXiv:1603.02698 [astro-ph.HE].[14] G. B. Rybicki, C. O. Heinke, R. Narayan, and J. E.Grindlay, Astrophys. J. , 1090 (2006), arXiv:astro-ph/0506563.[15] F. ¨Ozel, T. Guver, and D. Psaltis, Astrophys. J. ,1775 (2009), arXiv:0810.1521 [astro-ph].[16] S. Bogdanov, G. B. Rybicki, and J. E. Grindlay, Astro-phys. J. , 668 (2007), arXiv:astro-ph/0612791.[17] M. C. Miller, Astrophys. J. , 27 (2016),arXiv:1602.00312 [astro-ph.HE].[18] A. Watts et al. , PoS AASKA14 , 043 (2015),arXiv:1501.00042 [astro-ph.SR].[19] E. Epelbaum, H.-W. Hammer, and U.-G. Meissner,Rev. Mod. Phys. , 1773 (2009), arXiv:0811.1338[nucl-th].[20] R. Machleidt and D. Entem, Phys. Rept. , 1 (2011),arXiv:1105.2919 [nucl-th].[21] H.-W. Hammer, S. Knig, and U. van Kolck, Rev. Mod.Phys. , 025004 (2020), arXiv:1906.12122 [nucl-th].[22] I. Tews, Z. Davoudi, A. Ekstrm, J. D. Holt, and J. E.Lynn, J. Phys. G , 103001 (2020), arXiv:2001.03334[nucl-th].[23] I. Tews, T. Krger, K. Hebeler, and A. Schwenk, Phys.Rev. Lett. , 032504 (2013), arXiv:1206.0025 [nucl- th].[24] K. Hebeler, J. Lattimer, C. Pethick, and A. Schwenk,Astrophys. J. , 11 (2013), arXiv:1303.4662 [astro-ph.SR].[25] G. Baardsen, A. Ekstrm, G. Hagen, andM. Hjorth-Jensen, Phys. Rev. C , 054312 (2013),arXiv:1306.5681 [nucl-th].[26] C. Drischler, V. Soma, and A. Schwenk, Phys. Rev. C , 025806 (2014), arXiv:1310.5627 [nucl-th].[27] G. Hagen, T. Papenbrock, A. Ekstrm, K. Wendt,G. Baardsen, S. Gandolfi, M. Hjorth-Jensen, andC. Horowitz, Phys. Rev. C , 014319 (2014),arXiv:1311.2925 [nucl-th].[28] A. Carbone, A. Rios, and A. Polls, Phys. Rev. C ,044302 (2013), arXiv:1307.1889 [nucl-th].[29] L. Coraggio, J. Holt, N. Itaco, R. Machleidt, L. Mar-cucci, and F. Sammarruca, Phys. Rev. C , 044321(2014), arXiv:1402.0965 [nucl-th].[30] C. Wellenhofer, J. W. Holt, N. Kaiser, and W. Weise,Phys. Rev. C , 064009 (2014), arXiv:1404.2136 [nucl-th].[31] A. Roggero, A. Mukherjee, and F. Pederiva, Phys. Rev.Lett. , 221103 (2014), arXiv:1402.1576 [nucl-th].[32] J. Holt and N. Kaiser, Phys. Rev. C , 034326 (2017),arXiv:1612.04309 [nucl-th].[33] C. Drischler, A. Carbone, K. Hebeler, and A. Schwenk,Phys. Rev. C , 054307 (2016), arXiv:1608.05615 [nucl-th].[34] A. Ekstrm, G. Hagen, T. Morris, T. Papenbrock,and P. Schwartz, Phys. Rev. C , 024332 (2018),arXiv:1707.09028 [nucl-th].[35] C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev.Lett. , 042501 (2019), arXiv:1710.08220 [nucl-th].[36] D. Lonardoni, I. Tews, S. Gandolfi, and J. Carlson,Phys. Rev. Res. , 022033 (2020), arXiv:1912.09411[nucl-th].[37] M. Piarulli, I. Bombaci, D. Logoteta, A. Lovato,and R. Wiringa, Phys. Rev. C , 045801 (2020),arXiv:1908.04426 [nucl-th].[38] K. Hebeler, J. Holt, J. Menendez, and A. Schwenk,Ann. Rev. Nucl. Part. Sci. , 457 (2015),arXiv:1508.06893 [nucl-th].[39] C. Drischler, W. Haxton, K. McElvain, E. Mereghetti,A. Nicholson, P. Vranas, and A. Walker-Loud (2019)arXiv:1910.07961 [nucl-th].[40] F. Sammarruca and R. Millerson, Front. in Phys. , 213(2019).[41] BUQEYE collaboration (2020) https://buqeye.github.io/software/ .[42] C. Drischler, R. Furnstahl, J. Melendez, and D. Phillips(2020) arXiv:2004.07232 [nucl-th]. [43] C. Drischler, J. Melendez, R. Furnstahl, and D. Phillips(2020) arXiv:2004.07805 [nucl-th].[44] M. Leonhardt, M. Pospiech, B. Schallmo, J. Braun,C. Drischler, K. Hebeler, and A. Schwenk (2019)arXiv:1907.05814 [nucl-th].[45] H. Tan, J. Noronha-Hostler, and N. Yunes (2020)arXiv:2006.16296 [astro-ph.HE].[46] Y. Lim, A. Bhattacharya, J. W. Holt, and D. Pati(2020) arXiv:2007.06526 [nucl-th].[47] I. Tews, P. T. Pang, T. Dietrich, M. W. Coughlin,S. Antier, M. Bulla, J. Heinzel, and L. Issa (2020)arXiv:2007.06057 [astro-ph.HE].[48] R. Essick and P. Landry (2020) arXiv:2007.01372 [astro-ph.HE].[49] A. Tsokaros, M. Ruiz, and S. L. Shapiro (2020)arXiv:2007.05526 [astro-ph.HE].[50] F. Fattoyev, C. Horowitz, J. Piekarewicz, and B. Reed(2020) arXiv:2007.03799 [nucl-th].[51] D. A. Godzieba, D. Radice, and S. Bernuzzi (2020)arXiv:2007.10999 [astro-ph.HE].[52] A. Kanakis-Pegios, P. Koliogiannis, and C. Moustakidis(2020) arXiv:2007.13399 [nucl-th].[53] S. Koranda, N. Stergioulas, and J. L. Friedman, Astro-phys. J. , 799 (1997), arXiv:astro-ph/9608179.[54] R. C. Tolman, Phys. Rev. , 364 (1939).[55] J. Oppenheimer and G. Volkoff, Phys. Rev. , 374(1939).[56] J. M. Lattimer, M. Prakash, D. Masak, and A. Yahil,Astrophys. J. , 241 (1990).[57] J. W. Negele and D. Vautherin, Nucl. Phys. A , 298(1973).[58] G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. , 299 (1971).[59] J. Melendez, S. Wesolowski, and R. Furnstahl, Phys.Rev. C , 024003 (2017), arXiv:1704.03308 [nucl-th].[60] J. Melendez, R. Furnstahl, D. Phillips, M. Pratola,and S. Wesolowski, Phys. Rev. C , 044001 (2019),arXiv:1904.10581 [nucl-th].[61] T. Krger, I. Tews, K. Hebeler, and A. Schwenk, Phys.Rev. C , 025802 (2013), arXiv:1304.2212 [nucl-th].[62] D. Entem, R. Machleidt, and Y. Nosyk, Phys. Rev. C , 024004 (2017), arXiv:1703.05454 [nucl-th].[63] K. Hebeler and A. Schwenk, Phys. Rev. C , 014314(2010), arXiv:0911.0483 [nucl-th].[64] P. Bevington and D. Robinson, Data Reduction and Er-ror Analysis for the Physical Sciences , 3rd ed. (McGraw-Hill, 2003).[65] J. D. Evans,
Straightforward Statistics for the Behav-ioral Sciences (Brooks/Cole Publishing, Pacific Grove,Calif., 1996).[66] A. G. Asuero, A. Sayago, and A. G. Gonz´alez, Crit.Rev. Anal. Chem. , 41 (2006).[67] C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev.C , 054314 (2016), arXiv:1510.06728 [nucl-th].[68] N. Kaiser, Phys. Rev. C , 065201 (2015),arXiv:1504.00604 [nucl-th].[69] C. Wellenhofer, J. W. Holt, and N. Kaiser, Phys. Rev.C , 055802 (2016), arXiv:1603.02935 [nucl-th].[70] R. Somasundaram, C. Drischler, I. Tews, and J. Mar-gueron (2020) arXiv:2009.04737.[71] T. Damour and A. Nagar, Phys. Rev. D , 084035(2009), arXiv:0906.0096 [gr-qc].[72] T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read,Phys. Rev. D , 123016 (2010), arXiv:0911.3535 [astro- ph.HE].[73] S. Postnikov, M. Prakash, and J. M. Lattimer, Phys.Rev. D , 024016 (2010), arXiv:1004.5098 [astro-ph.SR].[74] T. Dietrich et al. , Phys. Rev. D , 024029 (2019),arXiv:1804.02235 [gr-qc].[75] T. Zhao and J. M. Lattimer, Phys. Rev. D , 063020(2018), arXiv:1808.02858 [astro-ph.HE].[76] S. De, D. Finstad, J. M. Lattimer, D. A. Brown,E. Berger, and C. M. Biwer, Phys. Rev. Lett. ,091102 (2018), [Erratum: Phys.Rev.Lett. 121, 259902(2018)], arXiv:1804.08583 [astro-ph.HE].[77] P. Landry and R. Essick, Phys. Rev. D , 084049(2019), arXiv:1811.12529 [gr-qc].[78] 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), arXiv:1908.10352 [astro-ph.HE].[79] T. Zhao and J. M. Lattimer, Phys. Rev. D , 023021(2020), arXiv:2004.08293 [astro-ph.HE].[80] I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C ,045804 (2018), arXiv:1804.02783 [nucl-th].[81] J. M. Lattimer and M. Prakash, Phys. Rept. , 127(2016), arXiv:1512.07820 [astro-ph.SR].[82] C. C. Moustakidis, T. Gaitanos, C. Margaritis, andG. Lalazissis, Phys. Rev. C , 045801 (2017), [Erratum:Phys.Rev.C 95, 059904 (2017)], arXiv:1608.00344 [nucl-th].[83] C. Margaritis, P. Koliogiannis, and C. Moustakidis,Phys. Rev. D , 043023 (2020), arXiv:1910.05767[nucl-th].[84] M. G. Alford and S. Han, Eur. Phys. J. A , 62 (2016),arXiv:1508.01261 [nucl-th].[85] K. Chatziioannou and S. Han, Phys. Rev. D , 044019(2020), arXiv:1911.07091 [gr-qc].[86] S. Han and A. W. Steiner, Phys. Rev. D , 083014(2019), arXiv:1810.10967 [nucl-th].[87] P. Haensel and J. Zdunik, Nature , 617 (1989).[88] J. M. Lattimer and M. Prakash, “What a Two So-lar Mass Neutron Star Really Means,” in From Nu-clei to Stars: Festschrift in Honor of Gerald E Brown ,edited by S. Lee (World Scientific, 2011) pp. 275–304,arXiv:1012.3208 [astro-ph.SR].[89] S. Han and M. Prakash, Astrophys. J. , 2 (2020),arXiv:2006.02207 [astro-ph.HE].[90] K. Hebeler, J. Lattimer, C. Pethick, and A. Schwenk,Phys. Rev. Lett. , 161102 (2010), arXiv:1007.1746[nucl-th].[91] A. Bauswein, O. Just, H.-T. Janka, and N. Stergioulas,Astrophys. J. Lett. , L34 (2017), arXiv:1710.06843[astro-ph.HE].[92] A. L. Watts et al. , Rev. Mod. Phys. , 021001 (2016),arXiv:1602.01081 [astro-ph.HE].[93] L. McLerran and S. Reddy, Phys. Rev. Lett. ,122701 (2019), arXiv:1811.12503 [nucl-th].[94] E. D. Van Oeveren and J. L. Friedman, Phys. Rev. D , 083014 (2017), arXiv:1701.03797 [gr-qc].[95] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, andR. Schaeffer, Nuclear Physics A , 231 (1998).[96] I. Tews, J. M. Lattimer, A. Ohnishi, andE. E. Kolomeitsev, Astrophys. J. , 105 (2017),arXiv:1611.07133 [nucl-th].[97] J. M. Lattimer and Y. Lim, Astrophys. J. , 51(2013), arXiv:1203.4286 [nucl-th]. [98] L. Lindblom, Phys. Rev. D , 103011 (2010),arXiv:1009.0738 [astro-ph.HE].[99] L. Lindblom and N. M. Indik, Phys. Rev. D , 084003(2012), arXiv:1207.3744 [astro-ph.HE].[100] L. Lindblom and N. M. Indik, Phys. Rev. D , 064003(2014), [Erratum: Phys.Rev.D 93, 129903 (2016)],arXiv:1310.0803 [astro-ph.HE].[101] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Fried-man, Phys. Rev. D , 124032 (2009), arXiv:0812.2163[astro-ph].[102] C. J. Horowitz and J. Piekarewicz, Phys. Rev. C ,062802 (2001), arXiv:nucl-th/0108036.[103] T. E. Riley et al. , Astrophys. J. Lett. , L21 (2019),arXiv:1912.05702 [astro-ph.HE].[104] M. Miller et al. , Astrophys. J. Lett. , L24 (2019),arXiv:1912.05705 [astro-ph.HE].[105] M. Ruiz, S. L. Shapiro, and A. Tsokaros, Phys. Rev. D , 021501 (2018), arXiv:1711.00473 [astro-ph.HE].[106] L. Rezzolla, E. R. Most, and L. R. Weih, Astrophys. J.Lett. , L25 (2018), arXiv:1711.00314 [astro-ph.HE].[107] B. P. Abbott et al. (LIGO Scientific, Virgo), Class.Quant. Grav. , 045006 (2020), arXiv:1908.01012 [gr-qc].[108] D. Radice, A. Perego, F. Zappa, and S. Bernuzzi, Astro-phys. J. Lett. , L29 (2018), arXiv:1711.03647 [astro-ph.HE].[109] M. W. Coughlin et al. , Mon. Not. Roy. Astron. Soc. , 3871 (2018), arXiv:1805.09371 [astro-ph.HE].[110] K. Kiuchi, K. Kyutoku, M. Shibata, and K. Taniguchi,Astrophys. J. Lett. , L31 (2019), arXiv:1903.01466[astro-ph.HE].[111] P. Landry, R. Essick, and K. Chatziioannou, Phys. Rev.D , 123007 (2020), arXiv:2003.04880 [astro-ph.HE].[112] R. Essick, I. Tews, P. Landry, S. Reddy, and D. E. Holz(2020) arXiv:2004.07744 [astro-ph.HE].[113] J.-L. Jiang, S.-P. Tang, Y.-Z. Wang, Y.-Z. Fan, and D.-M. Wei, Astrophys. J. , 1 (2020), arXiv:1912.07467[astro-ph.HE].[114] M. Al-Mamun, A. W. Steiner, J. Nttil, J. Lange,R. O’Shaugnessy, I. Tews, S. Gandolfi, C. Heinke, andS. Han (2020) arXiv:2008.12817 [astro-ph.HE].[115] T. Zhao and J. M. Lattimer, in preparation (2020).[116] D. Reardon et al. , Mon. Not. Roy. Astron. Soc. ,1751 (2016), arXiv:1510.04434 [astro-ph.HE].[117] P. Danielewicz, R. Lacey, and W. G. Lynch, Science , 1592 (2002), arXiv:nucl-th/0208016.[118] C. Constantinou, B. Muccioli, M. Prakash, andJ. M. Lattimer, Phys. Rev. C , 025801 (2015),arXiv:1504.03982 [astro-ph.SR].[119] B. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J.Lett. , L3 (2020), arXiv:2001.01761 [astro-ph.HE].[120] R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. , L44 (2020), arXiv:2006.12611 [astro-ph.HE].[121] R. J. Foley, D. A. Coulter, C. D. Kilpatrick, A. L.Piro, E. Ramirez-Ruiz, and J. Schwab, Mon. Not. Roy.Astron. Soc. , 190 (2020), arXiv:2002.00956 [astro-ph.HE].[122] G. Baym and S. Chin, Phys. Lett. B , 241 (1976).[123] Y. Nambu and G. Jona-Lasinio, Phys. Rev. , 345(1961).[124] A. Kurkela, P. Romatschke, and A. Vuorinen, Phys.Rev. D , 105021 (2010), arXiv:0912.1856 [hep-ph].[125] M. Alford, M. Braby, M. Paris, and S. Reddy, Astro-phys. J. , 969 (2005), arXiv:nucl-th/0411016. [126] T. Kojo, P. D. Powell, Y. Song, and G. Baym, Phys.Rev. D , 045003 (2015), arXiv:1412.1108 [hep-ph].[127] T. Kl¨ahn and T. Fischer, Astrophys. J. , 134 (2015),arXiv:1503.07442 [nucl-th].[128] R. Gomes, P. Char, and S. Schramm, Astrophys. J. , 139 (2019), arXiv:1806.04763 [nucl-th].[129] S. Han, M. Mamun, S. Lalit, C. Constantinou,and M. Prakash, Phys. Rev. D , 103022 (2019),arXiv:1906.04095 [astro-ph.HE].[130] G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song,and T. Takatsuka, Rept. Prog. Phys. , 056902 (2018),arXiv:1707.04966 [astro-ph.HE].[131] E. Annala, T. Gorda, A. Kurkela, J. Nttil, andA. Vuorinen, Nature Phys. (2020), 10.1038/s41567-020-0914-9, arXiv:1903.09121 [astro-ph.HE].[132] J. Hoppe, C. Drischler, K. Hebeler, A. Schwenk,and J. Simonis, Phys. Rev. C , 024318 (2019),arXiv:1904.12611 [nucl-th].[133] T. Hther, K. Vobig, K. Hebeler, R. Machleidt,and R. Roth, Phys. Lett. B , 135651 (2020),arXiv:1911.04955 [nucl-th].[134] S. Binder, J. Langhammer, A. Calci, and R. Roth,Phys. Lett. B , 119 (2014), arXiv:1312.5685 [nucl-th].[135] V. Lapoux, V. Som, C. Barbieri, H. Hergert, J. Holt,and S. Stroberg, Phys. Rev. Lett. , 052501 (2016),arXiv:1605.07885 [nucl-ex].[136] E. Epelbaum, H. Krebs, and P. Reinert, Front. in Phys. , 98 (2020), arXiv:1911.11875 [nucl-th].[137] A. Dyhdalo, R. Furnstahl, K. Hebeler, and I. Tews,Phys. Rev. C , 034001 (2016), arXiv:1602.08038 [nucl-th].[138] J. Hoppe, C. Drischler, R. Furnstahl, K. Hebeler,and A. Schwenk, Phys. Rev. C , 054002 (2017),arXiv:1707.06438 [nucl-th].[139] B. Carlsson, A. Ekstrm, C. Forssn, D. Strmberg,G. Jansen, O. Lilja, M. Lindby, B. Mattsson,and K. Wendt, Phys. Rev. X , 011019 (2016),arXiv:1506.02466 [nucl-th].[140] S. Wesolowski, N. Klco, R. Furnstahl, D. Phillips,and A. Thapaliya, J. Phys. G , 074001 (2016),arXiv:1511.03618 [nucl-th].[141] S. Wesolowski, R. Furnstahl, J. Melendez, andD. Phillips, J. Phys. G , 045102 (2019),arXiv:1808.08211 [nucl-th].[142] K. Huang, J. Hu, Y. Zhang, and H. Shen (2020)arXiv:2008.04491 [nucl-th].[143] E. R. Most, L. J. Papenfort, L. R. Weih, and L. Rezzolla(2020) arXiv:2006.14601 [astro-ph.HE].[144] N.-B. Zhang and B.-A. Li (2020) arXiv:2007.02513[astro-ph.HE].[145] V. Dexheimer, R. Gomes, T. Kl¨ahn, S. Han, andM. Salinas (2020) arXiv:2007.08493 [astro-ph.HE].[146] A. Sedrakian, F. Weber, and J.-J. Li, Phys. Rev. D , 041301 (2020), arXiv:2007.09683 [astro-ph.HE].[147] J. Boguta, Phys. Lett. B , 255 (1981).[148] J. Lattimer, M. Prakash, C. Pethick, and P. Haensel,Phys. Rev. Lett. , 2701 (1991).[149] A. Cherman, T. D. Cohen, and A. Nellore, Phys. Rev.D , 066003 (2009), arXiv:0905.0903 [hep-th].[150] J. M. Maldacena, Int. J. Theor. Phys. , 1113 (1999),arXiv:hep-th/9711200.[151] C. Ecker, C. Hoyos, N. Jokela, D. Rodr´ıguez Fern´andez,and A. Vuorinen, JHEP , 031 (2017), arXiv:1707.00521 [hep-th].[152] T. Ishii, M. J¨arvinen, and G. Nijs, JHEP , 003(2019), arXiv:1903.06169 [hep-ph].[153] P. M. Hohler and M. A. Stephanov, Phys. Rev. D ,066002 (2009), arXiv:0905.0900 [hep-th]. [154] P. Romatschke and U. Romatschke, Relativistic FluidDynamics In and Out of Equilibrium , Cambridge Mono-graphs on Mathematical Physics (Cambridge UniversityPress, 2019) arXiv:1712.05815 [nucl-th].[155] P. Bedaque and A. W. Steiner, Phys. Rev. Lett. ,031103 (2015), arXiv:1408.5116 [nucl-th].[156] I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astro-phys. J.860