New equations of state constrained by nuclear physics, observations, and QCD calculations of high-density nuclear matter
NNew equations of state constrained by nuclear physics,observations, and high-density QCD calculations
S. Huth,
1, 2, ∗ C. Wellenhofer,
1, 2, † and A. Schwenk
1, 2, 3, ‡ Technische Universit¨at Darmstadt, Department of Physics, 64289 Darmstadt, Germany ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum f¨ur Schwerionenforschung GmbH, 64291 Darmstadt, Germany Max-Planck-Institut f¨ur Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany
We present new equations of state for applications in core-collapse supernova and neutron starmerger simulations. We start by introducing an effective mass parametrization that is fit to recentmicroscopic calculations up to twice saturation density. This is important to capture the predictedthermal effects, which have been shown to determine the proto-neutron star contraction in super-nova simulations. The parameter range of the energy-density functional underlying the equation ofstate is constrained by chiral effective field theory results at nuclear densities as well as by func-tional renormalization group computations at high densities based on QCD. We further implementobservational constraints from measurements of heavy neutron stars, the gravitational wave signalof GW170817, and from the recent NICER results. Finally, we study the resulting allowed rangesfor the equation of state and for properties of neutron stars, including the predicted ranges for theneutron star radius and maximum mass.
I. INTRODUCTION
The equation of state (EOS) of dense matter is of cen-tral interest in nuclear physics and astrophysics. In themulti-messenger era, gravitational wave observations ofneutron star mergers such as GW170817 [1, 2] as well asmeasurements of neutron star radii from NASA’s NICERmission [3–5] are providing novel constraints on the EOSof neutron star matter.In recent years, many efforts have been undertaken toconstruct EOS parametrizations (see, e.g., Refs. [6–10])based on microscopic calculations of pure neutron mat-ter (PNM) [11–14]. However, such calculations are onlyavailable up to 1 − n , such that extrapolations to high densities using, e.g.,piecewise polytropes [6] or a speed of sound parametriza-tion [9, 10] have to be applied. High-density constraintsfrom neutron star observations, in particular precise massmeasurements of heavy neutron stars [15, 16], reduce theparameter space of possible extensions significantly.Compared to cold isolated neutron stars, core-collapsesupernovae (CCSNe) and neutron star mergers (NSM)probe a much wider range of the ( n, x, T ) space of theEOS, where n is the baryon density, x the proton frac-tion, and T the temperature. The EOSs commonly usedin numerical CCSN and NSM simulations are based ei-ther on phenomenological Skyrme energy-density func-tionals such as the Lattimer-Swesty EOS [17] or rela-tivistic mean-field models like the Shen EOS [18]. Morerecently, new EOSs for astrophysical simulations basedon these two approaches have been constructed [19–22].However, these phenomenological EOSs are often notconsistent with microscopic calculations at zero temper-ature [6, 11–14, 23] and nuclear thermodynamics [24–26]. ∗ Email: [email protected] † Email: [email protected] ‡ Email: [email protected]
Thermal effects are well characterized by the nucleon ef-fective mass m ∗ t ( n, x ) [26], where t = n , p for neutronsand protons. The effective mass is a crucial quantity inCCSN simulations that governs the proto-neutron starcontraction [27, 28]. Furthermore, recently a first step to-wards systematic computations of the EOS at intermedi-ate to high densities was achieved by Leonhardt et al. [29]who calculated symmetric nuclear matter (SNM) at den-sities 3 (cid:46) n/n (cid:46)
10 by using functional renormaliza-tion group (fRG) methods based on QCD. Constructingnew EOS parametrizations that are consistent with state-of-the-art constraints from nuclear physics, astrophysicalobservations, and high-density QCD calculations will en-able important progress in nuclear astrophysics.In this work, we develop a novel EOS functional thatincorporates recent micropscopic results for the nucleoneffective mass [26]. The parameters of the EOS func-tional are then fit to theoretical calculations at low andhigh densities as well as observational constraints frommass measurements, GW170817, and NICER. The sys-tematics of different parameter choices is analyzed in de-tail, and based on our EOS functional we derive compre-hensive uncertainty bands for the EOS and for neutronstar properties. Our work provides the groundwork fornew EOSs of hot and dense matter for CCSN and NSMsimulations.This paper is organized as follows. In Sec. II we re-view current constraints on the EOS by nuclear physics,observations, and high-density QCD calculations. Theeffective mass parametrization together with the EOSfunctional and its characteristics are the topic of Sec. III.In Sec. IV we discuss the different fit procedures usedto constrain the parameters of the EOS functional andstudy in detail the impact of variations of the parame-ters. Based on this, in Sec. V we examine our predictionsfor various neutron star properties and thermal effects.Finally, in Sec. VI we conclude with a short summaryand outlook. a r X i v : . [ nu c l - t h ] S e p II. OVERVIEW OF EQUATION OF STATECONSTRAINTS
Here, we briefly summarize presently available con-straints on the EOS of dense nuclear matter. First, inSec. II A we examine constraints from nuclear physics onthe properties of neutron-rich matter at densities up to1 − A. Constraints from nuclear physics
In this section we summarize constraints on the EOSfrom nuclear theory and experiment. In particular, wediscuss constraints on various characteristic parametersof the EOS around saturation density n : the bindingenergy B , the incompressibility K , the symmetry energycoefficient E sym , and the slope parameter L . (The effec-tive mass is discussed in Sec. III A.) Further, we discussconstraints from theoretical calculations on the EOS ofPNM.
1. Approximate EOS parametrizations
The bulk of dense matter in CCSN and NSM consistsof nucleons and electrons (as well as a small fraction ofmuons). For isolated neutron stars the temperature isnegligible with respect to nuclear energy scales. Thus,the neutron star EOS is predominantly determined bythe ground-state energy of strongly-interacting nucleonicmatter E ( n, β ). In neutron stars, the isospin asymmetry β = ( n n − n p ) /n is for a given nucleon number den-sity n = n n + n p fixed by the conditions of beta equi-librium and charge neutrality. A very useful approxima-tive parametrization of the dependence of E ( n, β ) on β is given by EA ( n, β ) ≈ EA ( n,
0) + S ( n ) β . (1)Here, E ( n,
0) is the energy of SNM, E ( n,
1) correspondsto PNM, and S ( n ) is called the symmetry energy. It hasbeen validated in various microscopic nuclear many-bodycalculations [30–35] that Eq. (1) provides a very goodapproximation of the exact β dependence over the entirerange β ∈ [0 ,
1] for densities up to n (cid:46) (1 − n . Thesymmetry energy S ( n ) may be obtained as the quadratic The accuracy of Eq. (1) is however expected to decrease withincreasing density [34, 35]. coefficient in the expansion of E ( n, β ) in β . From Eq. (1)it follows that the quadratic coefficient is approximatelyequal to the difference in energy of PNM to SNM:12 ∂ E/A∂β (cid:12)(cid:12)(cid:12) β =0 ≈ EA ( n, − EA ( n, . (2)Throughout this work, we assume a quadratic depen-dence of the interaction part of E ( n, β ) on β , so forthe interaction part the two sides of Eq. (2) are iden-tified. The kinetic part is modeled as a noninteractinggas of neutron and protons with density- and isospin-dependent effective masses. Within this approximation,the isospin-asymmetry dependence of the kinetic part istreated exactly. This approach is backed up by the micro-scopic calculations of isospin-asymmetric nuclear matterof Refs. [25, 32, 34]. See Sec. III for details on the form ofour EOS functionals. For the symmetry energy, we thenchoose the definition S ( n ) ≡ EA ( n, − EA ( n, . (3)(See, e.g., Ref. [34] for an analysis of the difference be-tween Eq. (3) and the second-order Taylor coefficient inthe β expansion.)Another useful parametrization is obtained by expand-ing the energy per particle in terms of the relative densitydifference to saturation density η = ( n − n ) / n : EA ( n, β ) ≈ − B + 12 Kη + ( E sym + Lη ) β . (4)Here, B = − E/A | n = n ,β =0 is the binding energy of SNM.Since the pressure P = n ∂ n E/A vanishes in SNM at n , the leading density dependence in the expansion of E ( n,
0) in terms of η corresponds to the incompressibility K . The expansion of the symmetry energy S ( n ) intro-duces the symmetry energy coefficient E sym = S ( n ) andthe slope parameter L . The latter is proportional to thepressure of PNM at saturation density, i.e., L = 3 n ∂S∂n (cid:12)(cid:12)(cid:12) n = 3 n P ( n , β = 1) . (5)The incompressibility at saturation density is propor-tional to the second derivative of the energy per particlewith respect to density: K = 9 ∂P∂n (cid:12)(cid:12)(cid:12) n ,β =0 = 9 n ∂ E/A∂n (cid:12)(cid:12)(cid:12) n ,β =0 . (6)The four parameters B , K , E sym and L provide a thor-ough characterization of the neutron star EOS at densi-ties around saturation density. No odd powers of β occur in the expansion if charge-symmetrybreaking effects are neglected. Notably, beyond second orderlogarithmic terms ∼ β n (cid:62) ln | β | may appear at zero tempera-ture [34, 36]. .
05 0 . .
15 0 . n [fm − ]0510152025 E / N [ M e V ] Hebeler et al ., ApJ (2013)Tews et al ., PRL (2013)Lynn et al ., PRL (2016)Drischler et al ., PRL (2019)Drischler et al ., GP-BGezerlis , Carlson , PRC (2010)Unitary gas ( ξ = 0 . .
05 0 . .
15 0 . n [fm − ]012345 P [ M e V f m − ] Hebeler et al ., ApJ (2013)Tews et al ., PRL (2013)Lynn et al ., PRL (2016)Drischler et al ., PRL (2019)Drischler et al ., GP-B
FIG. 1. Energy per particle
E/N (left panel) and pressure P (right panel) of PNM as a function of density n from variousmany-body calculations with chiral EFT interactions [6, 12–14, 37], see text for details. In the left panel, we also show thelow-density quantum Monte Carlo results by Gezerlis and Carlson [38] as well as the conjectured lower bound given by theenergy per particle of a unitary Fermi gas of neutrons [39]. For the saturation density n and energy B , we use n = 0 . B = 15 . K = 215(40) [14, 33, 40]. Many effortshave been devoted to determining the symmetry param-eter E sym and L , which play a crucial role in neutronstar structure and dynamics. In particular, microscopiccalculations of the EOS of PNM put tight constraintson E sym and L . We discuss these PNM and symmetryenergy constraints in the two paragraphs below.The temperature dependence of the EOS is key forCCSN and NSM applications. A very useful character-istic of the temperature dependence is given by the so-called thermal index Γ th , which is defined as [41]Γ th ( T, n, β ) = 1 + P ( T, n, β ) − P (0 , n, β ) ε ( T, n, β ) − ε (0 , n, β ) , (7)with the internal energy density ε = E/V . The ther-mal index of a noninteracting nucleon gas with density-dependent effective mass m ∗ ( n ) is given by [42]Γ th ( n ) = 53 − nm ∗ ( n ) ∂m ∗ ( n ) ∂n . (8)It has been verified in microscopic nuclear-matter calcu-lations that the form given by Eq. (8) provides a veryprecise approximation of Γ th ( T, n, β ) for β ∈ { , } , n (cid:46) n and T (cid:46)
30 MeV [26, 43]. A crucial noveltyin our EOS functionals is therefore the accurate imple-mentation of the effective masses of neutrons and protons m ∗ n , p ( n, β ). We discuss this and the available constraintson m ∗ n , p ( n, β ) in Sec. III A. Our results for the thermalindex are then examined in Sec. V B.
2. Neutron matter constraints
The modern approach to the description of the stronginteraction at nuclear energy scales is based on chiraleffective field theory (EFT) and renormalization group(RG) methods [44–46]. From general EFT convergencerestrictions as well as regulator and many-body conver-gence considerations, the viability of this approach is re-stricted to densities n (cid:46) n . The theoretical uncertain-ties in current implementations of chiral interactions ina given many-body framework arise from the interplayof finite-regulator artifacts, many-body and EFT trun-cation errors, and parameter-fitting ambiguities.Because nuclear forces are weaker in PNM, the theo-retical uncertainties are under better control there com-pared to SNM. In Fig. 1 we compare the results for theenergy per particle and pressure of PNM obtained fromseveral recent nuclear many-body calculations with chi-ral EFT interactions. The results by Hebeler et al. [6],Tews et al. [12], and Drischler et al. [14, 37] are basedon many-body perturbation theory, while the results byLynn et al. [13] were obtained from auxiliary-field dif-fusion Monte Carlo computations using local chiral in-teractions. In each case, the results include uncertaintyestimates, shown as bands in Fig. 1. These are basedon EFT truncation errors and different regulators inRefs. [13, 14, 37], while they are mainly due to uncertain-ties in the low-energy couplings that enter three-nucleonforces in Refs. [6, 12]. The uncertainty band of Drischler et al. PRL (2019) [14] is based on simple EFT truncationerrors. The results of Drischler et al.
GP-B [37] are con-structed from the same calculations [from Ref. [14]] butbased on a Bayesian uncertainty analysis using Gaussianprocesses, which leads to a very similar band for the com-bined GP-B (450) and (500) results. One sees that whileoverall the results from these calculations are in goodagreement, the uncertainties become sizeable for densi-ties n (cid:38) n .At densities near and above saturation density the un-certainties associated with the effective description of thenuclear interactions dominate over many-body trunca-tion effects. At low densities n (cid:28) n , the nuclear inter-actions are less intricate, but here the many-body accu-racy may be inflicted by the sensitivity to large-scatteringlength physics. Still, as shown in Fig. 1, the various chiralEFT-based many-body calculations discussed above arein reasonable agreement with the low-density results fromprecise quantum Monte Carlo computations by Gezerlisand Carlson [38].Finally, in Fig. 1 we also show the energy per particleof a unitary Fermi gas of neutrons E UG ( n ) = ξE FG ( n ),where E FG ( n ) is the free neutron gas energy and theBertsch parameter is ξ ≈ .
376 [47]. In Ref. [39], it wasargued that E UG ( n ) can be used as a lower bound for thePNM energy. As seen in Fig. 1, the unitary gas boundreduces the uncertainties in the results by Tews et al. [12]and Lynn et al. [13], while the ones by Hebeler et al. [6]and Drischler et al. [14, 37] are unaffected.
3. Symmetry energy constraints
Equation 4 provides the basis for a benchmark anal-ysis of related aspects of nuclear physics and astro-physics [48, 49]. The symmetry energy coefficient E sym and slope parameter L are crucial for a variety of phe-nomena in nuclear physics and astrophysics, ranging fromnuclear masses [50], neutron skins and the dipole polar-izability [51–55], to heavy-ion collisions [56], and core-collapse supernovae [22, 27]. In particular, the neutronstar radius scales with the pressure of PNM at saturationdensity [49, 57], i.e., with the L parameter [see Eq. (5)].Many experimental and theoretical efforts have beenundertaken to constrain the symmetry energy. In Fig. 2,we show results for the correlation between E sym and L obtained from the microscopic PNM calculations byHebeler et al . (H) [6], Tews et al . (TK) [12], andDrischler et al . (GP-B (450), GP-B (500)) [37] discussedin Sec. II A 2 above. In addition, we show results fromauxiliary-field diffusion Monte Carlo calculations by Gan-dolfi et al . (G) [58]. The uncertainties in the H, TK, and FIG. 2. Theoretical constraints on the symmetry energy E sym and the slope parameter L from Refs. [6, 12, 37, 58], see textfor details. Also shown are constraints extracted from fits tonuclear masses (orange band) [50] and the conjectured uni-tary gas (UG) bound [39]. The corners (green dots) of thegreen shaded area (“This work”) correspond to the four rep-resentative ( E sym , L ) pairs adopted in this work. The figureis adapted from Refs. [37, 49] using the Jupyter notebookprovided in [37, 59]. G results were obtained by using various (chiral) two-and three-nucleon interactions. The TK results involvethe largest uncertainties, which stems in part from largervariations of the low-energy couplings in three-nucleonforces. In the GP-B case, we show results obtained fromchiral potentials with two different cutoffs, GP-B (450)and (500), where in each case the uncertainties were ob-tained from a Bayesian analysis using Gaussian processesof the fixed-cutoff EFT systematics. Note that in Fig. 1the two GP-B bands are combined in one single band.The region of ( E sym , L ) values spanned by these the-oretical results overlaps with the constraints extractedfrom various experiments; see, e.g., Refs. [37, 49, 60].As an example, in Fig. 2 we show the constraint ex-tracted from nuclear masses [50]. Also shown is the con-jectured unitary gas (UG) bound [39]. The theoreticalconstraints (H,K,TK,GP-B) are all consistent with theUG constraint.For the EOS constructed in this work we choose fourrepresentative ( E sym , L ) pairs that lie within the com-bined theoretical constraints (H,K,TK,GP-B):( E sym , L ) / MeV ∈ (cid:8) (30 , , (31 , , (33 , , (34 , (cid:9) . (9)These four pairs are shown as green dots in Fig. 2 andthe (green shaded) region spanned by them is labeled“This work”. Note that we exclude very large symmetryenergies and slope parameters in order to avoid havingtoo many EOSs exceed the PNM uncertainty band shownin Fig. 1. A detailed study of the EOS and neutronstar properties associated with our four ( E sym , L ) pairsis provided in Sec. IV A. B. Constraints from neutron star observations
Neutron star observations play a crucial role in con-straining the dense-matter EOS. In particular, mass mea-surements of two-solar-mass neutron stars [15, 16, 61]have narrowed the uncertainties in the neutron star mass-radius relation considerably. To support neutron starsof such mass the EOS cannot be too soft, which chal-lenges neutron star models that include substantial por-tions of exotic condensates or deconfined quark matter.The present lower bound for the maximal mass M max isgiven by the mass of the heaviest observed neutron stars:PSR J0740+6620 with a mass of M = 2 . +0 . − . M (cid:12) [16]at the 2 σ level measured using relativistic Shapiro delay.This is in line with the radio-timing observation of thepulsar J0348+0432 with M = 2 . ± . M (cid:12) [15]. In ourwork, we use the averaged lower bound of M = 1 . M (cid:12) as a constraint for the lower bound of the maximal mass.Since the observation of the first NSM GW170817 [1, 2]by LIGO/Virgo together with the observation of the cor-responding kilonova AT2017gfo and the short gamma-rayburst GRB170817A [62], there have been many efforts toinfer an upper bound on the maximum neutron star mass M max from the remnant behavior. The suggested limitsare generally in the range M max (cid:46) . − . M (cid:12) [63–69],which would rule out overly stiff EOS.Even with a relatively narrow range on the maximalmass, the radius of a typical neutron star with M =1 . M (cid:12) is uncertain, with a typical conservative range10 (cid:46) R/ km (cid:46)
14, see, e.g., Ref. [6]. Recently, a ma-jor step towards precise radius measurements was madeby the NICER collaboration [3, 4], which simultaneouslydetermined the radius and mass of PSR J0030+0451 viaX-ray pulse-profile modeling.Implications of this measurement on the EOS havebeen studied by Raaijmakers et al. [5] by applyingtwo parametrizations for the neutron star EOS (in β -equilibrium): a piecewise polytropic (PP) model [6] anda speed of sound (CS) parametrization [10]. Raaijmakers et al. [70] performed a joint analysis of these models toinfer implications on the EOS from the NICER measure-ment, GW170817, and the 2 . M (cid:12) mass measurement.Their results for the pressure as a function of density areshown in Fig. 3 (green bands). While the PP and CS Even with ( E sym , L ) values inside the GP-B region our EOS func-tional can still lead to PNM properties, which are incompatiblewith the theoretical PNM uncertainty band. This is because ofhigher-order terms in the density behavior. n/n P [ M e V f m − ] PP model , β -equil . CS model , β -equil . fRG , SNM
FIG. 3. Constraints on the pressure of neutron star mat-ter as a function of density n/n (green bands) from a jointanalysis [70] of the 2 . M (cid:12) mass constraint, GW170817, andthe NICER results, obtained using two different EOS models:piecewise polytropes (PP) and speed of sound model (CS),see Sec. II B for details. The bands are for the 95% credibleregions. Also shown are the results for the pressure of SNM(orange band) from the fRG study of Ref. [29], see Sec. II C. bands are consistent with each other, the PP model al-lows stiffer EOSs for densities n (cid:46) n and in generalsmaller maximal densities. Consequently, the M - R re-lation of the CS model involves somewhat smaller radiicompared to the PP model. In our work, we use the com-bined PP and CS bands by Raaijmakers et al. [70] as aconstraint for our EOS parametrization. C. Theoretical calculations at high densities
The ultra-high-density regime ( n (cid:38) n ) of the EOScorresponds to deconfined quark matter. PerturbativeQCD provides the expansion of the EOS about the high-density limit [71]. This expansion can be used to con-struct astrophysical EOSs from interpolating betweenthe chiral EFT band for the EOS at nuclear densitiesand the perturbative QCD region [72, 73]. As discussedin Sec. II B, in the present paper we incorporate high-density constraints from neutron star observations explic-itly. Moreover, we base the high-density extrapolation ofthe EOS of SNM on a recent fRG calculation at more rel-evant densities. Therefore, in our case the perturbativeQCD expansion does not provide significant additionalconstraints on the dense-matter EOS.At present no reliable and accurate method exists forcomputing the properties of strongly interacting matterat densities n (cid:38) n (apart from the perturbative QCDlimit). However, a notable step towards systematic high-density calculations was made by Leonhardt et al. inRef. [29]. Starting from the QCD action, they use thefRG to derive a low-energy quantum effective action witheffective four-quark interactions and diquark degrees offreedom. The uncertainties in the results for the (zero-temperature) EOS of SNM from this approach have beenestimated in terms of their RG scale dependence. Othersources of error are, e.g., due to neglected quark flavorsand higher-order interaction effects.The fRG results of Leonhardt et al. [29] for the pressureof SNM are shown in Fig. 3 (orange band). They span aband from n = 3 n to n = 10 n . The band lies mostlybelow the observational neutron star matter constraintsfrom Raaijmakers et al. [70], with small overlaps for den-sities near n = 3 n and near n = 8 n . Note that thefRG band is significantly smaller than the neutron starmatter bands of Raaijmakers et al. . In this paper, we willsee that neutron star constraints, in particular the maxi-mum mass constraint M max (cid:62) . M (cid:12) , tend to favor apressure of SNM that lies somewhat above the fRG bandeither near n = 3 n or near n = 8 n ; see Secs. III B, IV Cand V. Nevertheless, we find several EOS that are con-sistent with neutron star observations and for which thepressure of SNM lies within the fRG band for n (cid:38) n . III. NEW EQUATION OF STATE FUNCTIONAL
We now come to the construction of a new EOSfunctional that takes into account the constraints fromnuclear physics, neutron star observations and high-density QCD calculations summarized in Sec. II. First,in Sec. III A we examine recent microscopic calculationsof the neutron and proton effective mass m ∗ n , p ( n, β ) inSNM ( β = 0) and PNM ( β = 1), and introduce a con-venient parametrization of m ∗ n , p ( n, β ) to be implementedin our EOS functionals. The construction of the EOSfunctional is the subject of Sec. III B. A. Temperature dependence andnucleon effective mass
Recently, Carbone and Schwenk [26] computed thefinite-temperature EOS of PNM and SNM from chiralEFT interactions using the self-consistent Green’s func-tion method. In addition, they also calculated the ef-fective masses m ∗ n , p ( n, β = 0 , th obtainedfrom the pressure and the energy density, see Eq. (7),can be accurately parametrized in terms of the effectivemass, via the form given by Eq. (8). Recent microscopicneutron-matter calculations in many-body perturbationtheory have confirmed this result [43]. Therefore, a reli-able implementation of the effective masses of neutrons and protons m ∗ n , p ( n, β ) is crucial to capture thermal ef-fects in astrophysical applications.To this end, we introduce an effective massparametrization that fits the results for m ∗ n , p ( n, β = 0 , n (cid:46) n from Ref. [26] based on the N LONN potential from Ref. [74] and N LO 3N interactionsconstructed in Ref. [75]. The behavior of m ∗ n , p ( n, β ) athigher densities is uncertain. We explore different scenar-ios in this regime. Our effective mass parametrization asa function of density is given by m ∗ t m = 1 + (cid:16) α n t + β n − t + α n / t + β n / − t + α n / t + β n / − t (cid:17)
11 + e n + (cid:16) (cid:15) t n t n + (cid:15) − t n − t n − (cid:17) − e − n e − n − n off ) , (10)where the nucleon with opposite isospin is denoted by − t . The six parameters α i , β i with i ∈ (1 ,
3) are fit to theSNM and PNM results of Ref. [26]. The factor 1 / (1+ e n )(sigmoid function) has the effect that the fitted part goesto zero with increasing density. The high-density behav-ior of m ∗ n , p ( n, β ) is then fixed by the parameters (cid:15) t and (cid:15) − t as well as by the offset n off of the (modified) logis-tic function (1 − e − n ) / (1 + e − n − n off ) ). For instance,for β = 1 the high-density limit of the neutron effectivemass is given by (cid:15) n , and for β = 0 the nucleon effectivemass approaches the value ( (cid:15) n + (cid:15) p ) /
2. Note that thisimplies that the high-density limit of the thermal indexis Γ th → /
3, see Eq. (8). The correct ultra-relativisticlimit is Γ th → /
3, but this matters for the nucleonic partof the EOS only for densities far above those relevant forneutron stars.In Fig. 4, we show the results for the nucleon effec-tive mass in PNM and SNM from Ref. [26] and threerepresentative effective mass parametrizations based onEq. (10). (In the PNM case we show the results for theneutron effective mass m ∗ n ( n, β = 1).) At low densities,the effective mass is a decreasing function of density, butstarting at around nuclear saturation density it increaseswith density for both PNM and SNM mainly due to 3Nforces. The effective mass in PNM is larger than the onefor SNM (see also Ref. [76]) and for densities n (cid:38) . n it exceeds the bare nucleon mass.From Eq. (8) it follows that the thermal index isΓ th < / th < th < th > n off = 0 . − and restrict our-selves to cases where (cid:15) t = (cid:15) − t = (cid:15) in Eq. (10) suchthat the effective mass has the same high-density limit in n/n . . . . . . . . m ∗ / m PNM , Carbone , Schwenk , PRC (2019)SNM , Carbone , Schwenk , PRC (2019) m ∗ . m ∗ . m ∗ . PNMSNM
FIG. 4. Effective mass m ∗ /m as a function of density n/n forPNM and SNM. The blue (PNM) and red (SNM) solid linesup to n/n = 2 show the results of Carbone and Schwenk [26].The grey lines (connected by colored bands) correspond to thethree representative effective mass parametrizations employedin this work. The associated high-density limits are given by m ∗ /m → . m ∗ /m → . m ∗ /m → . SNM and PNM. We employ three representative valuesof the high-density limit, i.e., (cid:15) ∈ { . , . , . } , denotedby m ∗ . , m ∗ . , and m ∗ . in Fig. 4. As seen in Fig. 4, thesethree scenarios span a reasonably wide range for the be-havior of the effective mass at high densities. Note alsothat the fit below twice saturation density is unaffectedby the high-density behavior (with our parametrizationthis holds true even for more extreme values of (cid:15) ).Regarding isospin-asymmetric nuclear matter, our ef-fective mass parametrization ensures that the neutroneffective mass m ∗ n ( n, β ) increases with β and satisfies m ∗ n ( n, β ) > m ∗ p ( n, β ), in agreement with theoretical re-sults [34, 77–80]. In contrast to m ∗ n ( n, β ), which is con-strained by fits to microscopic calculations for both PNMand SNM, in our approach the β dependence (at finite n ) of the proton effective mass m ∗ p ( n, β ) is an outcomeof the fit of m ∗ p ( n,
0) = m ∗ n ( n,
0) to SNM (after m ∗ n ( n, β dependenceof both m ∗ n and m ∗ p is largest at n ≈ − n (thehigh-density limit (cid:15) is β independent). Moreover, ourparametrization leads to a β dependence of m ∗ p ( n, β ) thatis decreased compared to that of m ∗ n ( n, β ), which is con-sistent with the results from Refs. [78, 79]. For m ∗ . ,the proton effective mass decreases with β at low densi-ties n (cid:46) n and increases for n (cid:38) n . In the m ∗ . and m ∗ . case the proton effective mass decreases with β atall densities, with the decrease being significantly more pronounced for m ∗ . . Here, the behavior for m ∗ . and m ∗ . is more in line with nuclear theory results [78, 79].Future work may involve the construction of improvedEOS functionals that incorporate additional theoreticalconstraints on the β dependence of m ∗ p ( n, β ). B. Equation of state functional
We now introduce the new EOS functional that formsthe basis for the investigations carried out in the remain-der of this paper. The microscopic results for the rela-tion between the effective mass m ∗ t ( n, β ) and the thermalindex Γ th ( n, β ) discussed in Secs. II A and III A makeclear that a reasonable approach to the temperature de-pendence of an effective EOS functional is to use a T -dependent kinetic term with density-dependent effectivemass and a T -independent interaction part. This is alsosupported by the microscopic nuclear-matter calculationsof Refs. [24–26, 34, 43] where it was found that the T de-pendence of the interaction contribution in many-bodyperturbation theory is small compared to the one of thenoninteracting contribution.In the usual (Skyrme) energy density functionals, theinteraction part is modeled as a finite polynomial in frac-tional powers of density, see, e.g., Ref. [22, 81]. By con-struction, the high-density behavior of a polynomial EOSansatz involves a highly fine-tuned balance between dif-ferent density powers. In certain cases, i.e., for somejudicious choices of the density powers, a reasonablehigh-density extrapolation can result from fits to micro-scopic calculations at nuclear densities [8, 81]. However,for the systematic construction of EOS functionals con-strained by nuclear physics, neutron star observations,and high-density QCD calculations, a polynomial ansatzcan clearly encounter difficulties.We therefore choose a form of the interaction part thatameliorates this fine tuning. For the internal energy den-sity as a function of density n = n n + n p , proton fraction x = n p /n and temperature T we use the following form: EV ( n, x, T ) = (cid:88) t τ t ( n, x, T )2 m ∗ t ( n, x ) − xn ∆+ (cid:88) i (cid:20) a i d a + n ( δ i − / + 4 b i x (1 − x ) d b + n ( δ i − / (cid:21) n δ i / . (11)Here, the second term gives the rest mass contribution(modulo the neutron mass energy), with ∆ the neutron–proton mass difference. The first term corresponds to thekinetic part of the internal energy density; it is modeledas a noninteracting gas of neutrons and protons with ef-fective masses m ∗ n , p ( n, x ) given by Eq. (4). That is, theterm τ t is given by τ t ( n, x, T ) = 12 π (cid:90) ∞ dp p ×
11 + exp (cid:2) T (cid:0) p m ∗ t ( n,x ) − ˜ µ t ( n, x, T ) (cid:1)(cid:3) , (12)where the auxiliary chemical potential ˜ µ t ( n, x, T ) is de-fined via n t = 12 π (cid:90) ∞ dp p
11 + exp (cid:2) T (cid:0) p m ∗ t − ˜ µ t (cid:1)(cid:3) . (13)The T → τ t ( n, x,
0) = (3 π n t ) / / (5 π ). The role of the auxiliarychemical potential is similar to the one in many-bodyperturbation theory at finite temperature [82]. The truechemical potential is obtained from the thermodynamicpotential corresponding to the variables ( n, x, T ), i.e., thefree energy. Since the interaction part and the effectivemasses are T independent, the free energy is obtained bysubstituting the kinetic part of E/V with the free en-ergy density of a (non-relativistic) noninteracting gas ofneutrons and protons with effective masses m ∗ n , p ( n, x ).The third term in Eq. (11) is the interaction part. Thecrucial feature of the interaction part is that it is based onrational functions instead of density monomials. Whilethe parameters a i and b i with i ∈ (1 ,
4) are fit to low-and high-density results as specified below, the densityexponents δ i as well as d a and d b are not fit parametersbut set to specific values. We choose two different setsfor δ i : δ k F = (3 , , , , δ n = (3 , , , . (14)For the choice δ k F , the density exponents in the numer-ators of the interaction part are (1 , / , / , k F at zero temperature. The choice δ n corresponds to in-teger powers of n in the numerators. The density de-pendence of the denominators is chosen such that in thehigh-density limit the interaction part becomes propor-tional to n / . The purpose of the denominators is tomitigate the fine-tuning between the different parts of theinteraction term such that the EOS functional is stableunder variations of the fit input. For a given choice of δ i ,the fit performance of the EOS functional is controlled by We use the nonrelativistic quasiparticle dispersion relation forall densities. The high-density behavior of our EOSs is fit toobservational and fRG constraints, so only thermal effects atvery high densities are affected by this approximation, whichis however a minor effect in comparison to the effective-massuncertainties in that regime. The density dependence in the ultra-relativistic limit is ∼ n / ,but this matters only for densities far above those relevant forneutron stars, see Sec. III A. the two offset parameters d a and d b . We set d a = d b = d and use for d the following values: d k F ∈ { , , , } , d n ∈ { . , . , . , . } . (15)These choices provide a reasonably wide range of dif-ferent density behaviors, as examined in detail below.The smaller values of d for δ n are mandated by the largedensity exponents, i.e., the suppression of higher densitypowers must set in earlier there. We note that remov-ing the restriction d a = d b has no notable impact on ourresults.We fix the eight parameters a , , , and b , , , bymatching to the following input: • the energy per particle of PNM at n = 0 .
05 fm − ,determined by the QMC result from Ref. [38] as E/N (0 .
05 fm − ) = 2 . • the nuclear matter properties ( n , B, K, E sym , L ), • the pressure of PNM at n = 1 .
28 fm − ≈ n , • the pressure of SNM at n = 1 .
28 fm − ≈ n .Here, the six nuclear-density inputs (first two items) arevaried according to their uncertainties, as examined inSec. II A. The high-density input is taken such that theresulting EOS is consistent with constraints from neutronstar observations [Sec. II B] and the pressure of SNM isin reasonable agreement with the fRG results [Sec. II C](we allow a 10% deviation from the fRG band to accountfor further fRG uncertainties).The results for the energy per particle E/A of PNMand SNM obtained for one particular input set are shownin Fig. 5. Specifically, here we set the nuclear matterproperties to (
B, K, E sym , L ) = (15 . , , ,
35) MeVand n = 0 .
157 fm − , the pressure at n = 1 .
28 fm − forSNM to 600 MeV fm − and for PNM to 1000 MeV fm − ,and use the effective mass scenario m ∗ . . Fig. 6 shows thecorresponding results for the pressure P and the squareof the speed of sound c s (in units where c = 1), whichis given by the derivative of the pressure with respect toenergy density. Variations of the input are investigatedin Sec. IV. For each δ i set [Eq. (14)], we show the resultsfor the smallest and the largest d in the correspondingset of possible d values [Eq. (15)]. One sees that thenuclear-density input has the effect that for densities upto roughly twice saturation density the functional is notvery sensitive to the values of δ i and d . At intermediatedensities 2 n (cid:46) n < n , different choices of δ i and d result in the following systematics: • δ k F , d = 1 . c s ≈ .
5. In the case of PNM, the EOS isstiffer in order to support a 2 M (cid:12) neutron star and c s shows a broad peak around 5 n . n/n E / A [ M e V ] δ k F , d . δ k F , d . δ n , d . δ n , d . FIG. 5. Results for the energy per particle of PNM (blue)and SNM (red) as a function of density n/n obtained fromthe two δ i sets with the minimal and maximal value of thecorresponding d . We use the same set of low- and high-densityfit points for all depicted EOS. The light blue band in the insetcorresponds to the combined chiral EFT results from Fig. 1. • δ k F , d = 7 . d parameter results in no notable changesfor SNM. In contrast, increasing d softens the EOSof PNM for n (cid:46) n , while at high densities it issignificantly stiffer and c s almost reaches the speedof light at 8 n . Nevertheless, the energy per parti-cle as well as the pressure of PNM is smaller in thegiven density range. • δ n , d = 0 . δ k F , the larger density exponents in δ n lead to arapid increase of the energy per particle and thepressure of both SNM and PNM at comparativelylow densities. As the density increases the EOSbecomes softer again in order to match the pres-sure fit point at 8 n . As a consequence, we find apronounced peak for the speed of sound due to thestiffness of the EOS. • δ n , d = 0 . d param-eter softens the EOS. In comparison to δ k F , thesensitivity of the functional with respect to d ismore pronounced for PNM as well as for SNM. Inthis specific case, the speed of sound has a secondmaximum at high densities; this feature depends onthe chosen input values and is not present in mostEOS.These characteristics of the EOS functional are quiterobust throughout the input parameter space. In linewith this, our new approach (rational functions instead n/n P [ M e V f m − ] δ k F , d . δ k F , d . δ n , d . δ n , d . fRG , SNM1 2 3 4 5 6 7 8 n/n . . . . c s FIG. 6. Same as Fig. 5 but here we show the pressure P andspeed of sound c s of PNM (blue) and SNM (red). Note thepressure fit points at 8 n . The orange band corresponds tothe fRG results of Ref. [29]. of density monomials) ensures that the fit parameters( a i , b i ) remain of reasonable size (i.e., there is no “un-natural” fine-tuning) under comprehensive variations ofthe low- and high-density input. More precisely, forthe 16128 input sets considered in Sec. V A, those thatare consistent with the imposed constraints from nuclearphysics and observations with the density exponents δ n have absolute values of the dimensionless parameters ofat most 2.6, with mean values between 0.2 and 1.1. For δ k F the parameters are in general larger due to the smallerdensity exponents in the numerators of the interactionterms, i.e., δ k F involves more tuning than δ n . The meanvalues of the fit parameters for δ k F lie in the range 1.4to 15.5. For the combined ( δ k F and δ n ) parameter space,80% of the constrained EOS have a i and b i in a range0from −
12 (lower bound b ) to 20 (upper bound b ), withfive out of the eight fit parameters spanning only a rangeat most from − IV. EQUATION OF STATE VARIATIONS
With the energy density functional in place, we per-form variations of the input choices in order to span arange of EOS that covers the uncertainties of the con-straints discussed in Sec. II. First, we discuss the vari-ations of nuclear matter properties in Sec. IV A. Thisinvolves the four representative ( E sym , L ) pairs shown inFig. 2 as well as three choices for the saturation proper-ties of SNM. In Sec. IV B, we then analyze the behavior ofthe EOS parametrization for the three effective mass sce-narios introduced in Sec. III A. This is followed by high-density variations of the pressure for PNM and SNM inSec. IV C. These variations are performed for each set ofexpansion coefficients δ i and each of the corresponding d values given by Eqs. (14) and (15), respectively. A. Variations of nuclear matter properties
In order to cover the uncertainties of the energy ofPNM from many-body calculations based on chiral EFTas elaborated in Sec. II A 2 and Sec. II A 3, four com-binations of the symmetry energy E sym and the slopeparameter L were identified, namely( E sym , L ) / MeV ∈ (cid:8) (30 , , (31 , , (33 , , (34 , (cid:9) . (16)These values (green dots in Fig. 2) cover a reasonablerange of the combined theoretical results for the E sym − L correlation, see Fig. 2. In fact, we have investigated awhole grid of ( E sym , L ) pairs that encompasses and ex-ceeds the green-shaded region in Fig. 2: the grid rangesare 28–36 MeV for the symmetry energy and 30–75 MeVfor the slope parameter (with step sizes 1 MeV and5 MeV). We have examined the EOS and neutron starproperties obtained from each ( E sym , L ) pair in this gridfor each ( δ i , d ) choice, each effective mass scenario, andeach of the different ( n , B, K ) values and high-densityinput specified in Sec. IV C. For every ( E sym , L ) pair wethen counted the number of EOSs, which fulfill the con-straints from nuclear physics and neutron star observa-tions discussed in Secs. II A and II B. The nuclear-densityinputs are discussed further below; see Sec. IV C for de-tails on the high-density input range considered. Notethat the high-density fRG results for SNM from Sec. II Care not enforced as a strict constraint.
28 29 30 31 32 33 34 35 36 E sym [MeV]30354045505560657075 L [ M e V ] δ k F δ n UG “This work” E O S FIG. 7. Grid of ( E sym , L ) pairs used to determine the fourrepresentative pairs given by Eq. (16). The black line (labeled”UG”) corresponds to the constraint on ( E sym , L ) obtainedfrom the unitary gas boundary on the PNM energy [39], seeSec. II A 3. The color coding gives the number of EOS thatfulfill the theoretical and observational constraints discussedin Secs. II, see the text for details. The results from this study are analyzed in Fig. 7where one sees that larger slope parameters become dis-favored the smaller the symmetry energy is. This feature,which is more pronounced for δ k F (see the boundary inFig. 7), is reflected also in the microscopic constraints onthe E sym − L correlation, see Fig. 2. Our four choices for( E sym , L ) are based on the combined results from Figs. 2and 7, and on the observation that they lead to EOSsthat cover a broad range of neutron star properties.For the saturation properties of SNM, we use the em-pirical saturation point n = 0 . − and B =15 . K = 215(40) MeV deter-mined from microscopic nuclear-matter calculations [14,33, 40]. The uncertainties in ( n , B, K ) are coveredby three combinations. The triple ( K, n , B ) central usesthe central values. The two other triples combinethe minimal (maximal) values of n and B with thelargest (smallest) incompressibility: ( K max , ( n , B ) min )and ( K min , ( n , B ) max ). Overall:( K min , ( n , B ) max ) = (175 , . , . , (17)( K, n , B ) central ) = (215 , . , . , (18)( K max , ( n , B ) min ) = (255 , . , . , (19)in units MeV, fm − , and MeV, respectively. Thesecombinations have a physical motivation: First, they fol-low the Coester-band correlation between n and B val-ues [14]. Second, if SNM saturates at small densities andenergies one expects that the incompressibility increases,and vice versa.1 E / N [ M e V ] E sym = 30 MeV , L = 35 MeV δ k F δ n K min , ( n s , B ) max ( K, n s , B ) central K max , ( n s , B ) min Unitary gas ( ξ = 0 . E sym = 31 MeV , L = 55 MeV .
05 0 . . n [fm − ]05101520 E / N [ M e V ] E sym = 33 MeV , L = 65 MeV .
05 0 . . n [fm − ] E sym = 34 MeV , L = 55 MeV .
15 0 . − − δ k F δ n FIG. 8. Results for the energy per particle of PNM at nuclear densities. The four panels correspond to the four representative( E sym , L ) pairs given by Eq. (16). In each panel we show the results obtained from the three ( n , B, K ) triples given byEqs. (17)–(19). The light/dark blue lines are for the δ n /δ k F functional. All depicted EOSs employ the smallest d available forthe respective choice of δ i , see Eq. (15), the effective mass scenario m ∗ . , and the same high-density fits for the pressure (seeFig. 9, top panel). In the inset of the second panel, we show the corresponding SNM results. The light blue band in each panelcorresponds to the combined chiral EFT results from Fig. 1. The anthracite dash-dotted line in each panel is the conjecturedlower bound (unitary gas with ξ = 0 . The four ( E sym , L ) pairs and three ( n , B, K ) triplesamount to twelve possible combinations of nuclear matterproperties for each choice of ( δ i , d ). The correspondingresults for the PNM energy at nuclear densities are exam-ined in Fig. 8 where each of the four panels is for one ofthe four ( E sym , L ) pairs. One sees that the chosen varia-tions of the nuclear matter properties provide a thoroughrepresentation of the PNM uncertainty band from chiralEFT. The depicted EOS are for one particular choice of( δ i , d ), m ∗ t ( n, x ) and the high-density input, as specifiedin the caption of Fig. 8. Variations of these propertiesbroaden the covered area further. We use the unitary gas bound to rule out some of the EOSs (gray lines),in particular among those are K min , ( n , B ) max and δ k F .The EOS with intermediate symmetry energies and largeslope parameters ( E sym , L ) / MeV = (31 ,
55) and (33 , δ i and d .2 n/n P [ M e V f m − ] K min , ( n s , B ) max ( K, n s , B ) central K max , ( n s , B ) min PNM , δ k F PNM , δ n SNM , δ k F SNM , δ n fRG , SNM1 2 3 4 5 6 7 8 n/n . . . c s FIG. 9. High-density analog of Fig. 8 for the pressure and thespeed of sound of PNM (blue) and SNM (red) as a functionof density. The orange band corresponds to the fRG SNMresults from Ref. [29].
B. Effective mass variation
As discussed in Sec. III B, we employ three differ-ent parametrizations of the nucleon effective mass, m ∗ . , m ∗ . , and m ∗ . , where the subscript denotes the respec-tive high-density limit, see Fig. 4. Here, in Figs. 10 and11 we examine the impact of the high-density behaviorof the effective mass on the zero-temperature EOS. Theimpact on thermal effects is studied in Sec. V B.In Figs. 10 and 11, the nuclear matter properties arefixed as ( K, n , B ) = ( K max , ( n , B ) min ) , (20)( E sym , L ) / MeV = (30 , . (21)The high-density input is fixed as specified in Fig. 9. We n/n E / A [ M e V ] m ∗ . m ∗ . m ∗ . PNM , δ k F PNM , δ n SNM , δ k F SNM , δ n FIG. 10. Results for the energy per particle of PNM (blue)and SNM (red) as a function of density for the three effectivemass scenarios of Fig. 4. All depicted EOS employ the small-est d available for the respective choice of δ i . The nuclearmatter properties and high-density input are fixed, see textfor details. The light blue band in the inset corresponds tothe combined chiral EFT results from Fig. 1. see that the overall influence of the effective mass on theenergy and pressure at zero temperature is comparativelysmall by construction. The speed of sound, as a quantitythat is not directly constrained by the fit, is more sen-sitive to variations of the effective mass. The observedbehavior depends on the choice of δ i , and in each casethe results for m ∗ . and m ∗ . are more similar comparedto m ∗ . . In the case of δ n the maximum of c s increaseswith the high-density limit of the effective mass (i.e., theEOS becomes stiffer), while for δ k F the speed of soundpeak occurs at smaller densities for m ∗ . . C. High-density variations
With a careful implementation of the nuclear physicsconstraints at hand, see Sec. IV A, the objective now is tohave the EOS functional reproduce the high-density con-straints from neutron star observations. That is, the goalis to cover the band for neutron star matter obtained byRaaijmakers et al. [70] and have EOS that are consistentwith the mass measurements of Antoniadis et al . [15] andthe 2 σ confidence interval of the 2 . M (cid:12) measurementby Cromartie et al. [16], see Sec. II B. For this, we spana grid of fit points for the pressure of SNM and PNM at n = 1 .
28 fm − ≈ n . We fit the pressure of SNM tovalues { , , , , , , } MeV fm − andthe pressure difference between PNM and SNM to3 n/n P [ M e V f m − ] m ∗ . m ∗ . m ∗ . PNM , δ k F PNM , δ n SNM , δ k F SNM , δ n fRG , SNM1 2 3 4 5 6 7 8 n/n . . . c s FIG. 11. Analog of Fig. 10 for the pressure and speed ofsound of PNM (blue) and SNM (red) as a function of density.The orange band corresponds to the fRG SNM result fromRef. [29]. { , , , , , , , } MeV fm − , so thepressure of PNM ranges from 350 to 1300 MeV fm − .This results in 56 high-density fit combinations for eachlow-density and effective mass input. From these, we ex-clude all EOSs that, after including β -equilibrium andelectrons, are not consistent with the Raaijmakers etal. bands in Fig. 3.The high-density fRG calculations of SNM by Leon-hardt et al. [29] (see Sec. II C) lie on the lower endof the employed fit values for the SNM pressure:310 MeV fm − (cid:46) P fRG (8 n ) (cid:46)
410 MeV fm − . A lowerSNM pressure implies that the pressure of matter in β -equilibrium is small as well. More specifically, the pro-ton fraction increases with the SNM-PNM pressure dif-ference, leading to a decrease of the pressure of matter n/n E / A [ M e V ] SNM , δ k F SNM , δ n PNM , δ k F PNM , δ n SNM , δ k F SNM , δ n PNM , δ k F PNM , δ n SNM , δ k F SNM , δ n PNM , δ k F PNM , δ n δ k F δ n FIG. 12. High-density variations for the energy per particleof PNM (blue) and SNM (red) as a function of density, seetext for details. The bands for the two different δ i combina-tions include only EOS that are consistent with neutron starconstraints. The d parameter, the effective mass, and the nu-clear matter properties are kept fixed. The light blue bandin the inset corresponds to the combined chiral EFT resultsfrom Fig. 1. in β -equilibrium. As a result, enforcing consistency withthe fRG results reduces the range for neutron star mat-ter to a great extent, such that the uncertainty band byRaaijmakers et al. [70] cannot be fully covered. There-fore, we do not use the fRG band as a strict constraint.The subset of EOSs that are consistent with the fRGcalculations is studied further in Sec. V A.As an example, representative high-density variationsof the energy per particle of PNM and SNM are shownin Fig. 12. The corresponding results for the pressureand the speed of sound are displayed in Fig. 13. Here,the nuclear matter properties are set to K max , ( n , B ) min and E sym = 30 MeV, L = 35 MeV. The effective massis given by m ∗ . , and for each δ i combination we use thesmallest d value from Eq. (15). For each δ i we keep onlythe EOSs, which are consistent with the constraints fromnuclear physics and neutron star observations.As seen in Fig. 12, for SNM the EOSs with δ n spana much wider energy band that almost entirely enclosesthe δ k F energy band. At nuclear densities the pressuresof both the δ n and the δ k F EOSs lie mostly above the fRGpressure, see the top panel of Fig. 13. At high densities onthe other hand the δ n high-density variations encompassthe entire fRG band. In contrast, for δ k F the deviationsfrom the fRG band increase with density.Compared to SNM, the δ n and δ k F energy and pressurebands are of similar size for PNM. Regarding the speed ofsound of PNM and SNM, in the bottom panel of Fig. 134 n/n P [ M e V f m − ] SNM , δ k F SNM , δ n PNM , δ k F PNM , δ n SNM , δ k F SNM , δ n PNM , δ k F PNM , δ n SNM , δ k F SNM , δ n PNM , δ k F PNM , δ n fRG , SNM1 2 3 4 5 6 7 8 n/n . . . c s FIG. 13. Analog of Fig. 12 for the pressure and speed ofsound of PNM (blue) and SNM (red) as a function of density.The orange band corresponds to the fRG SNM results fromRef. [29]. one sees that the high-density variations do not lead tosignificant changes in the systematics for different ( δ i , d )choices, see Fig. 6. For both δ i sets, the EOS with largeststiffness regions involve a PNM speed of sound that at itspeak is close to c s = 0 .
8. In Fig. 6, for both the PNM andSNM speed of sound the two δ i sets give non-overlappingresults at high densities. However, one needs to keep inmind the displayed results are for one particular choice ofthe d parameter, the nuclear matter properties and theeffective mass. Varying these reduces the area that is notcovered with the specific input used. Plots that involvethe full range of the considered parameter variations areshown in Sec. V. V. ASTROPHYSICAL EQUATION OF STATE
Here, in Sec. V A we examine our results for cold mat-ter in β -equilibrium and study neutron star propertiessuch as the M - R relation and the electron fraction. Wetake into account the full set of parameter variations ofthe EOS functional, as discussed in Sec. IV. Thermaleffects, which are crucial for applications in CCSN andNSM simulations, are the subject of Sec. V B. A. Neutron star properties
The density of electrons and muons in neutron starmatter is equal to the proton density because of localcharge neutrality. For simplicity, we neglect muons asthis causes only a very small change in the neutron starEOS. The proton fraction x at a given baryon density n is fixed by the requirement of β -equilibrium, i.e., µ n ( n, x ) = µ p ( n, x ) + µ e ( n e = xn ) , (22)where µ n , p , e is the chemical potential of the respectiveparticle species. Electrons can be modeled as an ultra-relativistic degenerate Fermi gas, so the electron pres-sure is P e = E e / (3 V ), with the electron energy den-sity given by E e /V = (3 π n e ) / / (4 π ). The elec-tron chemical potential reads µ e = (3 π n e ) / and thechemical potentials of neutrons and protons are µ n , p = ∂ n n , p E ( n, x ) /V + m n , p .With the variations of the EOS input and the choicesfor the functional parameters δ i and d of Eq. (14) andEq.(15) in place, we perform all possible fit combinationsto obtain bands for the EOS of matter in β -equilibrium.The neutron star mass-radius ( M - R ) relation is then ob-tained by solving the Tolman-Oppenheimer-Volkoff equa-tions [83, 84]. To this end, we implement the Baym,Pethick, Sutherland (BPS) crust from Ref. [85] for den-sities below n crust = 0 .
08 fm − , as in Ref. [70]. As dis-cussed in the two previous sections, for each of the δ i setswe consider four d values, twelve variations of nuclearmatter properties, three effective mass scenarios, and 56high-density fits for the pressure of SNM and PNM, re-sulting in a set of 16128 EOS. Among these, we keep onlyEOSs that • are consistent with the theoretical PNM uncer-tainty band and the unitary gas bound for the en-ergy per particle up to 0 . − , • provide masses of neutron stars of at least 1 . M (cid:12) (combined lower bound of the measurements fromRef. [15] and the 2 σ interval of Ref. [16]), • and lie within the 95% credible regions based onthe joint analysis of GW170817 and NICER fromRaaijmakers et al. [70].The results for the pressure and the speed of sound ofneutron star matter are shown in Fig. 14. We see that5 FIG. 14. Results for the pressure and the speed of sound ofmatter in β -equilibrium as a function of density. As describedin the text, we show all EOSs that fulfill the constraints fromnuclear physics and neutron star observations. The colorcoding indicates the mass of the corresponding neutron star,where dark green corresponds to masses up to 1 . M (cid:12) andlight green to higher masses up to the respective maximummass M max . Gray lines correspond to the continuation of theEOSs to densities above the central densities n c , max of therespective heaviest neutron star. The light-gray band depictsthe 95% credible region of the neutron star constraints fromRaaijmakers et al. [70]. our EOS functional covers almost the entire band for thepressure by Raaijmakers et al. [70]. At the high pressureboundary the agreement is very close, but some of thesofter EOSs within the Raaijmakers et al. band are notreproduced. This feature can be largely attributed tothe fact that we use a strict lower bound for the minimalvalue of M max whereas Raaijmakers et al. have modeledthe mass likelihood function for the 2 . M (cid:12) pulsar [16]. In Fig. 14, the parts of the different EOSs that corre-spond to neutron stars with masses below the canonical1 . M (cid:12) are highlighted in dark green. Their continua-tion up to the respective maximum mass M max is coloredin light green. The central density n . for a neutronstar with 1 . M (cid:12) lies approximately within 2 − . n .The smallest and largest n max are roughly 4 . n and7 . n , respectively, where one particular EOS reaches n max ≈ . n , very similar to Ref. [6]. Up to n . , thespeed of sound is relatively strongly constrained, but athigher densities a large variety of speed of sound curvesis present that covers a range from c s ≈ . M - R diagram. Regarding the M - R relation, compared toour results the band of Raaijmakers et al. [70] includes M (cid:38) . M (cid:12) neutron stars with slightly smaller radii.This is a direct consequence of the softer EOSs includedthere, as discussed above. For a 1 . M (cid:12) neutron star wefind a radius range of R . = 11 . − . (cid:46) R/ km (cid:46)
13. Further, in Fig. 15 wealso show the mass-radius band obtained by Hebeler etal. [6] using polytropic extrapolations of chiral EFT re-sults. Compared to the other bands, the band by Hebeler et al. [6] allows for neutron stars with smaller radii andlarger maximum masses. This is mainly because it showsthe entire region (100% instead of 95% credible) compat-ible with the maximum mass constraint.The density exponents of the functional δ i mostly in-fluence the stiffness of the EOS. In particular, softerEOSs corresponding to neutron stars with smaller radiimainly involve δ k F , while δ n yields stiffer EOSs and largerradii. Moreover, the back-bending of the M - R lines at M ≈ . M (cid:12) is more pronounced for EOSs that use δ n .The EOSs for which SNM is consistent with the fRGband by Leonhardt et al. [29] are highlighted in orange inFig. 15. More specifically, for the orange EOSs the pres-sure of SNM starting at 5 n deviates from the fRG bandby at most 10%. Compared to the full band, the fRG-consistent neutron star EOSs have lower pressures at en-ergy densities ε (cid:38)
800 MeV fm − , which is a consequenceof the relatively low SNM pressures obtained by the fRGcalculation. This translates into comparatively largerneutron star radii and smaller maximum masses. Nev-ertheless, the fRG-consistent EOSs cover a broad rangeof the pressure-energy density and M - R bands of Raaij-makers et al. [70]. Overall, the fRG results for SNM pro-vide viable additional constraints for astrophysical EOSconstructions, and incorporating them leads to a consid-erable narrowing of the uncertainty band for the EOSand the M - R relation. In the future, improved fRG cal-culations will enable further advancements along theselines.Our results for the neutron star maximum mass are6 FIG. 15. Analog of Fig. 14 for the pressure-energy density relation (left) and the mass-radius relation (right) of cold neutronstars. The orange lines correspond to EOSs that are consistent with the fRG SNM results from [29], see text for details. Thegray band depicts the 95% credible region of the neutron star constraints from Raaijmakers et al. [70]. For comparison, inthe M - R plot we show also the uncertainty band obtained by Hebeler et al. [6] using piecewise polytrope extensions to highdensities (light gray band). examined further in Fig. 16 where we plot the numberof EOSs that yield a given value of M max . The distribu-tion in Fig. 16 shows a broad peak, which falls off steeplyfor M max (cid:38) . M (cid:12) , reaching zero at M max ≈ . M (cid:12) .Our largest maximum masses are only slightly above themodel-dependent bound M max (cid:46) . − . M (cid:12) inferredfrom GW170817 [63–69]. Again, in Fig. 16 the EOSsthat are consistent with the fRG band are highlightedin orange. As discussed above, the fRG constraint im-plies softer EOSs and thus leads to smaller values of themaximum mass with M max (cid:46) . M (cid:12) .Finally, in Fig. 17 we show the electron fraction Y e ( n )in neutron stars as obtained from the different EOS. Thedensity dependence of the electron fraction is stronglyrelated to that of the symmetry energy S ( n ). Similarto the results for the speed of sound, our Y e ( n ) band isfairly narrow up to around saturation density where theelectron fraction is given by Y e ( n ) ≈ (4 E sym ) / π n ≈ . − .
055 [86]. Above (2 − n , the Y e ( n ) band widensconsiderably, with different EOS given electron fractionsbetween 0 to about 30% in the core of heavy neutronstars. B. Thermal effects
Some EOSs used in NSM simulations start from acold neutron star EOS and add a thermal part that isparametrized in terms of a density-independent thermalindex Γ th = const., e.g., Γ th = 1 −
2. The validity of
FIG. 16. Number of EOSs per maximum mass M max corre-sponding to the mass-radius relation of Fig. 15. The orangedotted line corresponds to EOSs that are consistent with thefRG SNM results from Ref. [29]. this approximation was studied by Bauswein et al. inRef. [41] who concluded that a consistent treatment ofthermal effects beyond the Γ th = const. approximation7 FIG. 17. Analog of Fig. 14 for the electron fraction Y e as afunction of density. is important. As discussed in Sec. II A 1, the thermal con-tribution to the EOS is governed to a large extent by thenucleon effective mass. Full EOS tables for astrophysi-cal simulations mostly apply a mean-field effective massthat monotonically decreases with density. This howeveris not consistent with microscopic nuclear-matter calcu-lations [26, 43], which show that interaction effects be-yond the mean-field approximation are significant. Asdiscussed in Sec. III A, our EOS functional incorporatessuch microscopic effective-mass results explicitly.The temperature dependence of our EOS functional isdescribed solely by the kinetic term, which is modeledas a noninteracting nucleon gas with neutron and protoneffective mass m ∗ n ( n, x ) and m ∗ p ( n, x ), see Sec. III B. Fromthis one arrives at the following formula for the thermalindex Γ th of isospin-asymmetric nuclear matter (ANM)with proton fraction x :Γ th ( n, x, T ) = 53 − n (cid:80) t =n , p ε t, th ( n,x,T ) m ∗ t ( n,x ) ∂m ∗ t ( n,x ) ∂n (cid:80) t =n , p ε t, th ( n, x, T ) , (23)where ε t, th ( n, x, T ) = ε t ( n, x, T ) − ε t ( n, x,
0) is the ther-mal part of the kinetic energy density of neutrons andprotons, respectively. In the context of our EOS func-tional, Eq. (23) is an exact representation of Γ th ( n, x, T ),i.e., it is equivalent to Eq. (7). The T dependenceof Γ th ( n, x, T ) vanishes for x = 0 and x = 1 /
2, i.e.,for PNM and SNM one obtains the familiar expres-sion given by Eq. (8). However, for ANM the ther-mal index is a temperature-dependent quantity. (Equa-tion (23) would reduce to Eq. (8) also for ANM if theeffective masses of protons and neutron were identified, m ∗ n ( n, x ) = m ∗ p ( n, x ), but this would be in conflict with nuclear theory results, see Sec. III A.)Consequently, we explore whether the T dependenceof Γ th for ANM is a significant effect, and, if the T de-pendence is small, what is the appropriate temperature-independent approximative expression for the thermalindex. For a classical free nucleon gas with m ∗ n ( n, x )and m ∗ p ( n, x ) one obtains by substituting 3 T n t / ε t, th ( n, x, T ) in Eq. (23) the expressionΓ th,classical ( n, x ) = 53 − (cid:88) t =n , p n t ( n, x ) m ∗ t ( n, x ) ∂m ∗ t ( n, x ) ∂n , (24)with n n = (1 − x ) n and n p = xn . Note that for x = 0and x = 1 /
2, respectively, both Eq. (23) and Eq. (24)yield the correct formula for PNM and SNM, Eq. (8).Comparing the results obtained from Eq. (23) andEq. (24) one finds that the classical expression pro-vides a very good approximation, with relative errorswell below the 1% level (except at very low T (cid:46) n/n ∈ [0 ,
8] and the results obtained from the three effectivemass scenarios ( m ∗ . , m ∗ . , and m ∗ . ), the mean rela-tive errors at T = 10 MeV are (0 . , . , . x = (0 . , . , . T ; at T = 1 MeV and T = 50 MeV theyare (0 . , . , . . , . , . x = (0 . , . , . T and x andeach effective mass scenario the deviations first increasewith density up to n/n ≈ −
7, and then decrease again(since the high-density limit of the effective mass is x in-dependent). Overall, we conclude that Γ th,classical ( n, x )provides a very good representation of the temperaturedependence of the EOS of ANM.Our results for the thermal index Γ th of PNM, SNM,and ANM with x = 0 . th isthen for each x determined entirely by that of m ∗ n ( n, x )and m ∗ p ( n, x ). For PNM and SNM, an increasing (de-creasing) effective mass implies that Γ th is below (above)the free or unitary Fermi gas value Γ th = 5 /
3, and thethermal index of PNM is larger than that of SNM. Foreach x , at low densities Γ th first increases with n andthen decreases again such that Γ th = 5 / n ≈ n , corresponding to the minimum of m ∗ t ( n, x )at around saturation density (see Sec. III A). The high-density behavior is fixed by the respective effective massscenario, where the largest deviations from Γ th = 5 / m ∗ . at n ≈ . n (the n → th → / th ≈ . n ≈ . n for SNM, are obtained for m ∗ . .The detailed description of thermal effects within ourEOS functional may have interesting effects in astro-physical applications. In particular, the proto-neutronstar contraction in CCSN simulations is largely governedby the T dependence of the EOS [22, 27]. Lower effec-tive masses lead to larger thermal contributions and thus8 FIG. 18. Results for the thermal index of PNM (blue), SNM(red), and ANM with x = 0 . th = 5 / to a larger PNS radius. A faster contraction increasesthe temperature at the surface of the PNS. As a con-sequence, neutrinos emitted from the PNS have largerenergies, which aids the shock evolution towards a fasterexplosion. All effective mass scenarios result in a ther-mal index that is mostly well below 5 / − n ,which may lead to a faster PNS contraction and explo-sion compared to commonly used astrophysical EOS suchas the Lattimer-Swesty or Shen EOS [17, 18, 87]. Thevery high-density regime of the EOS is more importantin NSM than in CCSN applications. Investigating theeffects of our different high-density effective mass scenar-ios in NSM simulations may be an interesting subject forfuture research. VI. SUMMARY AND OUTLOOK
In this paper, we have developed a new EOS functionalfor application in CCSN and NSM simulations. The EOSfunctional is constrained by various chiral-EFT based cal-culations of neutron matter and nuclear matter, by as-trophysical observations, and (to a lesser extent) by re-sults from a recent QCD-based fRG study of high-densitymatter. In particular, the EOSs obtained from our func-tional are consistent with recent mass measurements ofheavy neutron stars from [15, 16] and the joint analysisof observational data from GW170817 and NICER fromRaaijmakers et al. [70]. Using as input the recent microscopic nuclear-matterresults from Ref. [26], we have modeled the kinetic partof the EOS as a noninteracting nucleon gas with density-dependent effective mass. The careful implementationof microscopic results for the nucleon effective mass inour EOS functional is an essential novelty compared toprevious constructions of astrophysical EOS [17, 18, 20–22, 88], in particular the commonly used Lattimer-SwestyEOS [17] and Shen EOS [18].The effective mass is a key quantity that determinesthe temperature dependence of the EOS. It was demon-strated in various studies [27, 28, 41] that a thoroughdescription of the effective mass is crucial for CCSN andNSM simulations. In this work, we have studied differenthigh-density extrapolations of the microscopic nuclear-density results for the effective mass. It will be interest-ing to investigate the impact of different effective massscenarios on thermal effects in astrophysical applications.The description of the interaction part in our approachrepresents an improvement over traditional EOS func-tionals [17, 18, 20–22, 88] as well. That is, we have mod-eled the interaction part not as a sum of density monomi-als but as a sum of density-dependent rational functions(the temperature-independence of the interaction part issmall [43], and thus was neglected in this work). Thisensures that the EOS functional is appropriately stableunder variations of the low- and high-density input. Wehave fitted the parameters of the interaction part to thecombination of state-of-the-art neutron matter calcula-tions and observational constraints. From this, we havederived a comprehensive uncertainty band for neutronstar matter. Our EOSs predict that the radius of acanonical 1.4 solar mass neutron star lies in the range R . = 11 . − . ACKNOWLEDGMENTS
We thank Almudena Arcones, Arianna Carbone,Svenja Greif, Kai Hebeler, Jonas Keller, Yeunhwan Lim,Mirko Pl¨oßer, Geert Raaijmakers, and Ingo Tews for use-ful discussions. This work was supported by the DeutscheForschungsgemeinschaft (DFG, German Research Foun-dation) – Project-ID 279384907 – SFB 1245.9 [1] B. P. Abbott et al. (LIGO Scientific Collaboration andVirgo Collaboration), Phys. Rev. Lett. , 161101(2018).[2] B. P. Abbott et al. (LIGO Scientific Collaboration andVirgo Collaboration), Phys. Rev. X , 011001 (2019).[3] T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M.Ludlam, S. Guillot, Z. Arzoumanian, C. L. Baker, A. V.Bilous, D. Chakrabarty, et al. , Astrophys. J. Lett. ,L21 (2019).[4] M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov,Z. Arzoumanian, K. C. Gendreau, S. Guillot, A. K. Hard-ing, W. C. G. Ho, J. M. Lattimer, et al. , Astrophys. J.Lett. , L24 (2019).[5] G. Raaijmakers, T. E. Riley, A. L. Watts, S. K. Greif,S. M. Morsink, K. Hebeler, A. Schwenk, T. Hinderer,S. Nissanke, S. Guillot, et al. , Astrophys. J. Lett. ,L22 (2019).[6] K. Hebeler, J. M. Lattimer, C. J. Pethick, andA. Schwenk, Astrophys. J. , 11 (2013).[7] E. Rrapaj, A. Roggero, and J. W. Holt, Phys. Rev. C , 065801 (2016).[8] Y. Lim and J. W. Holt, Phys. Rev. Lett. , 062701(2018).[9] I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astro-phys. J. , 149 (2018).[10] S. K. Greif, G. Raaijmakers, K. Hebeler, A. Schwenk,and A. L. Watts, Mon. Not. Roy. Astron. Soc. , 5363(2019).[11] K. Hebeler and A. Schwenk, Phys. Rev. C , 014314(2010).[12] I. Tews, T. Kr¨uger, K. Hebeler, and A. Schwenk, Phys.Rev. Lett. , 032504 (2013).[13] J. E. Lynn, I. Tews, J. Carlson, S. Gandolfi, A. Gezerlis,K. E. Schmidt, and A. Schwenk, Phys. Rev. Lett. ,062501 (2016).[14] C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev.Lett. , 042501 (2019).[15] J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris,R. S. Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa,V. S. Dhillon, T. Driebe, et al. , Science , 1233232(2013).[16] H. T. Cromartie, E. Fonseca, S. M. Ransom, P. B. De-morest, Z. Arzoumanian, H. Blumer, P. R. Brook, M. E.DeCesar, T. Dolch, J. A. Ellis, et al. , Nature Astron. ,72 (2019).[17] J. M. Lattimer and F. D. Swesty, Nucl. Phys. A , 331(1991).[18] H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi, Nucl.Phys. A 637 , 435 (1998).[19] G. Shen, C. J. Horowitz, and S. Teige, Phys. Rev. C ,035802 (2011).[20] M. Hempel, T. Fischer, J. Schaffner-Bielich, andM. Liebendrfer, Astrophys. J. , 70 (2012).[21] A. W. Steiner, M. Hempel, and T. Fischer, Astrophys.J. , 17 (2013).[22] A. S. Schneider, L. F. Roberts, and C. D. Ott, Phys.Rev. C , 065802 (2017).[23] T. Kr¨uger, I. Tews, K. Hebeler, and A. Schwenk, Phys.Rev. C , 025802 (2013).[24] C. Wellenhofer, J. W. Holt, N. Kaiser, and W. Weise,Phys. Rev. C , 064009 (2014). [25] C. Wellenhofer, J. W. Holt, and N. Kaiser, Phys. Rev.C , 015801 (2015).[26] A. Carbone and A. Schwenk, Phys. Rev. C , 025805(2019).[27] H. Yasin, S. Schfer, A. Arcones, and A. Schwenk, Phys.Rev. Lett. , 092701 (2020).[28] A. S. Schneider, L. F. Roberts, C. D. Ott, andE. O’Connor, Phys. Rev. C , 055802 (2019).[29] M. Leonhardt, M. Pospiech, B. Schallmo, J. Braun,C. Drischler, K. Hebeler, and A. Schwenk,arXiv:1907.05814.[30] I. Bombaci and U. Lombardo, Phys. Rev. C , 1892(1991).[31] A. Carbone, A. Polls, C. Providˆencia, A. Rios, andI. Vida˜na, Eur. Phys. J. A , 13 (2014).[32] C. Drischler, V. Som`a, and A. Schwenk, Phys. Rev. C , 025806 (2014).[33] C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev. C , 054314 (2016).[34] C. Wellenhofer, J. W. Holt, and N. Kaiser, Phys. Rev.C , 055802 (2016).[35] R. Somasundaram, C. Drischler, I. Tews, and J. Mar-gueron, (2020), arXiv:2009.04737.[36] N. Kaiser, Phys. Rev. C , 065201 (2015).[37] C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R.Phillips, arXiv:2004.07232.[38] A. Gezerlis and J. Carlson, Phys. Rev. C , 025803(2010).[39] I. Tews, J. M. Lattimer, A. Ohnishi, and E. E. Kolomeit-sev, Astrophys. J. , 105 (2017).[40] K. Hebeler, S. K. Bogner, R. J. Furnstahl, A. Nogga, andA. Schwenk, Phys. Rev. C , 031301(R) (2011).[41] A. Bauswein, H.-T. Janka, and R. Oechslin, Phys. Rev.D , 084043 (2010).[42] C. Constantinou, B. Muccioli, M. Prakash, and J. M.Lattimer, Phys. Rev. C , 025801 (2015).[43] J. Keller, C. Wellenhofer, K. Hebeler, and A. Schwenk,in preparation (2020).[44] E. Epelbaum, H.-W. Hammer, and U.-G. Meißner, Rev.Mod. Phys. , 1773 (2009).[45] S. K. Bogner, R. J. Furnstahl, and A. Schwenk, Prog.Part. Nucl. Phys. , 94 (2010).[46] R. Machleidt and D. R. Entem, Phys. Rep. , 1 (2011).[47] M. J. H. Ku, A. T. Sommer, L. W. Cheuk, and M. W.Zwierlein, Science , 563 (2012).[48] M. B. Tsang, J. R. Stone, F. Camera, P. Danielewicz,S. Gandolfi, K. Hebeler, C. J. Horowitz, J. Lee, W. G.Lynch, Z. Kohley, et al. , Phys. Rev. C , 015803 (2012).[49] J. M. Lattimer and Y. Lim, Astrophys. J. , 51 (2013).[50] T. Kortelainen, M. Lesinski, J. Mor´e, W. Nazarewicz,J. Sarich, N. Schunck, M. V. Stoitsov, and S. Wild,Phys. Rev. C , 024313 (2010).[51] A. Tamii, I. Poltoratska, P. von Neumann-Cosel, Y. Fu-jita, T. Adachi, C. A. Bertulani, J. Carter, M. Dozono,H. Fujita, K. Fujita, et al. , Phys. Rev. Lett. , 062502(2011).[52] X. Roca-Maza, X. Vi˜nas, M. Centelles, B. K. Agrawal,G. Col`o, N. Paar, J. Piekarewicz, and D. Vretenar, Phys.Rev. C , 064304 (2015).[53] G. Hagen, A. Ekstr¨om, C. Forss´en, G. R. Jansen,W. Nazarewicz, T. Papenbrock, K. A. Wendt, S. Bacca, N. Barnea, B. Carlsson, et al. , Nat. Phys. , 186 (2016).[54] J. Birkhan, M. Miorelli, S. Bacca, S. Bassauer, C. A.Bertulani, G. Hagen, H. Matsubara, P. von Neumann-Cosel, T. Papenbrock, N. Pietralla, et al. , Phys. Rev.Lett. , 252501 (2017).[55] S. Kaufmann, J. Simonis, S. Bacca, J. Billowes, M. L.Bissell, K. Blaum, B. Cheal, R. F. Garcia Ruiz, W. Gins,C. Gorges, et al. , Phys. Rev. Lett. , 132502 (2020).[56] M. B. Tsang, Y. Zhang, P. Danielewicz, M. Famiano,Z. Li, W. G. Lynch, and A. W. Steiner, Phys. Rev. Lett. , 122701 (2009).[57] J. M. Lattimer and M. Prakash, Phys. Rep. , 109(2007).[58] S. Gandolfi, J. Carlson, and S. Reddy, Phys. Rev. C ,032801(R) (2012).[59] BUQEYE collaboration, https://buqeye.github.io/software/ (2020).[60] J. M. Lattimer, Nucl. Phys. A , 276 (2014).[61] P. Demorest, T. Pennucci, S. Ransom, M. Roberts, andJ. Hessels, Nature , 1081 (2010).[62] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese,K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Ad-hikari, V. B. Adya, et al. , Astrophys. J. Lett. , L12(2017).[63] B. Margalit and B. D. Metzger, Astrophys. J. , L19(2017).[64] M. Shibata, S. Fujibayashi, K. Hotokezaka, K. Kiuchi,K. Kyutoku, Y. Sekiguchi, and M. Tanaka, Phys. Rev.D (2017).[65] M. Ruiz, S. L. Shapiro, and A. Tsokaros, Phys. Rev. D (2018).[66] L. Rezzolla, E. R. Most, and L. R. Weih, Astrophys. J. , L25 (2018).[67] M. Shibata, E. Zhou, K. Kiuchi, and S. Fujibayashi,Phys. Rev. D (2019).[68] S. Ai, H. Gao, and B. Zhang, Astrophys. J. , 146(2020).[69] B. P. Abbott, R. Abbott, T. D. Abbott, S. Abraham,F. Acernese, K. Ackley, C. Adams, V. B. Adya, C. Af- feldt, M. Agathos, and et al., Class. Quant. Grav. ,045006 (2020).[70] G. Raaijmakers, S. K. Greif, T. E. Riley, T. Hinderer,K. Hebeler, A. Schwenk, A. L. Watts, S. Nissanke,S. Guillot, J. M. Lattimer, and et al., Astrophys. J. , L21 (2020).[71] T. Gorda, A. Kurkela, P. Romatschke, M. S¨appi, andA. Vuorinen, Phys. Rev. Lett. , 202701 (2018).[72] A. Kurkela, E. S. Fraga, J. Schaffner-Bielich, andA. Vuorinen, Astrophys. J. , 127 (2014).[73] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen,Phys. Rev. Lett. , 172703 (2018).[74] D. R. Entem and R. Machleidt, Phys. Rev. C ,041001(R) (2003).[75] P. Klos, A. Carbone, K. Hebeler, J. Men´endez, andA. Schwenk, Eur. Phys. J. A , 168 (2017), [Erratum:Eur. Phys. J. A , 76 (2018)].[76] K. Hebeler, T. Duguet, T. Lesinski, and A. Schwenk,Phys. Rev. C , 044321 (2009).[77] O. Sj¨oberg, Nucl. Phys. A , 511 (1976).[78] E. N. E. van Dalen, C. Fuchs, and A. Faessler, Phys.Rev. Lett. , 022302 (2005).[79] A. Li, J. N. Hu, X. L. Shang, and W. Zuo, Phys. Rev.C , 015803 (2016).[80] B.-A. Li, B.-J. Cai, L.-W. Chen, and J. Xu, Prog. Part.Nucl. Phys. , 29 (2018).[81] Y. Lim and J. W. Holt, Phys. Rev. C , 065805 (2017).[82] C. Wellenhofer, Phys. Rev. C , 065811 (2019).[83] R. C. Tolman, Phys. Rev. , 364 (1939).[84] J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. ,374 (1939).[85] G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. , 299 (1971).[86] K. Hebeler, J. M. Lattimer, C. J. Pethick, andA. Schwenk, Phys. Rev. Lett. , 161102 (2010).[87] J. Lattimer, C. Pethick, D. Ravenhall, and D. Lamb,Nucl. Phys. A , 646 (1985).[88] H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi, As-trophys. J. Suppl.197