Constraints on the nuclear symmetry energy from asymmetric-matter calculations with chiral NN and 3N interactions
LLA-UR-20-26894
Constraints on the nuclear symmetry energy from asymmetric-matter calculationswith chiral NN and 3N interactions
R. Somasundaram, ∗ C. Drischler,
2, 3, † I. Tews, ‡ and J. Margueron § Univ Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3,IP2I Lyon, UMR 5822, F-69622, Villeurbanne, France Department of Physics, University of California, Berkeley, CA 94720, USA Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA (Dated: September 11, 2020)The nuclear symmetry energy is a key quantity in nuclear (astro)physics. It describes the isospindependence of the nuclear equation of state (EOS), which is commonly assumed to be almostquadratic. In this work, we confront this standard quadratic expansion of the EOS with explicitasymmetric nuclear-matter calculations based on a set of commonly used Hamiltonians includingtwo- and three-nucleon forces derived from chiral effective field theory. We study, in particular,the importance of non-quadratic contributions to the symmetry energy, including the non-analyticlogarithmic term introduced by Kaiser [Phys. Rev. C , 065201 (2015)]. Our results suggestthat the quartic contribution to the symmetry energy can be robustly determined from the variousHamiltonians employed, and we obtain 1.00(8) MeV (or 0.55(8) MeV for the potential part) atsaturation density, while the logarithmic contribution to the symmetry energy is relatively smalland model-dependent. We finally employ the meta-model approach to study the impact of thehigher-order contributions on the neutron-star crust-core transition density, and find a small 5%correction. I. INTRODUCTION
The nuclear-matter equation of state (EOS) is of greatinterest for nuclear physics, see recent reviews [1, 2] andreferences therein. It connects bulk properties of atomicnuclei, with small isospin asymmetry, with neutron-richmatter inside neutron stars (NSs) [3, 4]. The isospindependence of the nuclear-matter EOS is described bythe nuclear symmetry energy which, for example, gov-erns the proton fraction in beta-equilibrium, determinesthe pressure in the core of NSs, and hence, the NS mass-radius relation [5–7], or cooling via the direct URCA pro-cess [8]. Due to its importance for many physical systems,the symmetry energy and its density dependence wereidentified as key quantities for nuclear (astro)physics inthe 2015 DOE/NSF Nuclear Science Advisory Commit-tee Long Range Plan for Nuclear Science [9], and areactively investigated by combining information from nu-clear theory, astrophysics, and experiments.Because NS observations still come with sizable un-certainties, the symmetry energy and its density depen-dence can not be inferred from NS properties alone [8].Hence, various constraints on the symmetry energy havebe been inferred from experimental data, e.g., precisedeterminations of neutron skins in lead (PREX) and cal-cium (CREX) [10, 11], collective modes such as giantdipole resonances [12], and heavy-ion collisions [13, 14].Typically, experimental constraints are in the range ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] e sym ( n sat ) ∼ −
35 MeV [8, 15, 16] at nuclear satu-ration density, n sat = 0 .
16 fm − ≡ n empsat (for a reviewsee, e.g., Ref. [13]). The determination of the symmetryenergy is on the road-map for several future experimentsconducted at rare-isotope beam facilities such as FRIBat MSU, SPIRAL2 at GANIL, and FAIR at GSI. Theo-retical calculations provide additional information on thesymmetry energy. While there are nuclear EOS modelsfor a wide range of values for the symmetry energy andits density dependence [17, 18], microscopic EOS calcula-tions have improved theoretical constraints over the lastyears [5–7, 19–22].A consistent extraction of the nuclear symmetry en-ergy from nuclear theory as well as experimental and as-trophysical programs requires that the measured quan-tities in these different approaches, as well as their rela-tions, are well defined. Different approximations for thesymmetry energy are commonly used. It is, therefore,important to clarify whether the symmetry energy mea-sured in laboratory experiments is the same quantity asthe one inferred from NS properties. For example, theenergy per particle of nuclear matter at zero tempera-ture is a function of the baryon density n = n n + n p and isospin asymmetry δ = ( n n − n p ) /n , where n n ( n p )denotes the neutron (proton) number density. The fol-lowing isospin-asymmetry expansion from symmetric nu-clear matter (SNM, δ = 0) to pure neutron matter (PNM, δ = 1) is often employed, e ( n, δ ) ≈ e ( n, δ = 0) + δ e sym,2 ( n )+ δ e sym,4 ( n ) + O ( δ ) . (1)Here e sym,2 ( n ) and e sym,4 ( n ) are the quadratic and quar-tic contributions to the symmetry energy, respectively.Given the expansion (1), the quadratic contribution to a r X i v : . [ nu c l - t h ] S e p the symmetry energy is defined by the second derivative e sym , ( n ) = 12 ∂ e ( n, δ ) ∂δ (cid:12)(cid:12)(cid:12)(cid:12) δ =0 , (2)similar to the empirical Bethe-Weizs¨acker mass formulafor finite nuclei. Hence, e sym,2 ( n ) is often referred to asthe symmetry energy , and used in nuclear experiments.In practice, however, the more commonly used definitionof the symmetry energy is given by the difference betweenthe energy per particle in pure neutron matter (PNM)and SNM, e sym ( n ) = e PNM ( n ) − e SNM ( n ) . (3)While this definition needs information on the EOS onlyin the limits of PNM and SNM, Eq. (2) necessitatesexplicit calculations of isospin-asymmetric nuclear mat-ter (ANM). Both e sym , ( n ) and e sym ( n ) are equal if theisospin dependence of the energy per particle is exactly quadratic, i.e., non-quadratic terms in the expansion (1)vanish. However, there is no a priori argument why theseterms should vanish. In fact, they have been found to berelevant for, e.g., accurate studies of nuclear matter inbeta-equilibrium at supra-saturation density [23–26] andthe crust-core transition density in NSs [25, 27].In this work, we study the expansion (1) and quantifythe impact of non-quadratic contributions to the sym-metry energy. We investigate to which extent uncertain-ties in the microscopical approach affect the extractionof non-quadratic contributions to the symmetry energy.The paper is organized as follows. In Sec. II, we discusssome previous studies of non-quadratic contributions tothe symmetry energy. In Sec. III, we present our com-putational setup and summarize the results calculatedin Ref. [20]. General expressions for the energy expan-sion in terms of the isospin-asymmetry parameter δ aregiven in Sec. IV, in particular, for the total energy perparticle as well as for the contributions of the potentialenergy terms. We then discuss the EOS in the limits ofSNM and PNM in Sec. V, followed by the symmetry en-ergy in Sec. VI. In Sec. VII, we study the impact of thenon-quadratic contributions to the symmetry energy indeterminations of the core-crust transition in NSs. Fi-nally, we conclude in Sec. VIII. The Python codes usedto perform the analysis and generate the figures in thispaper are publicly available on GitHub [28] and brieflydescribed in the Supplemental Material associated to thispublication. II. PREVIOUS STUDIES OF NON-QUADRATICCONTRIBUTIONS
As stated above, there is no a priori reason for theisospin-asymmetry expansion to be purely quadratic. Ingeneral, even the free Fermi gas (FFG) energy per parti-cle, given by e FFG ( n ) = t satSNM (cid:18) nn sat (cid:19) / h (1+ δ ) / +(1 − δ ) / i , (4) with t satSNM = m N (cid:16) π n sat (cid:17) / ≈ . e FFGsym , ( n ) ’ .
45 MeV × (cid:18) nn sat (cid:19) / , (5)represents a ∼ .
5% correction to the FFG symmetryenergy at n sat . Nuclear interactions also contribute tonon-quadratic terms. For example, the phenomenolog-ical Skyrme interaction [29] gives the following quarticcontribution to the symmetry energy: e Skyrmesym , ( n ) ’ e FFGsym , ( n )+ k π [3 t (1 + x ) + t (1 − x )] , (6)where the Skyrme parameters ( t , t ) represent the cor-rection to the bare nucleon mass generated by in-mediumeffects. Since the Skyrme in-medium mass is generally ∼ −
40% lower than the bare mass [29], these terms con-tribute to an increase of the FFG e sym , of ∼ −
40% to ∼ (0 . − .
8) MeV. In a recent work, Cai and Li [26] found e sym , ( n sat ) = (7 . ± .
5) MeV, which indicates a rathersignificant difference between e sym and e sym , . Theyemployed an empirically constrained isospin-dependentsingle-nucleon momentum distribution and the EOS ofPNM near the unitary limit. Subsequently, Bulgac et al. found that e sym , ( n = 0 . − ) = 2 .
635 MeV is neces-sary in order to reproduce properties of both finite nucleiand the PNM EOS as calculated in Ref. [30]. In con-trast, previous works, e.g., based on Brueckner-Hartree-Fock (BHF) approaches and hard-core interactions [31–34] obtained only small non-quadratic contributions tothe symmetry energy.In a recent study of nuclear matter in many-body per-turbation theory (MBPT) with two-nucleon (NN) andthree-nucleon (3N) interactions derived from chiral effec-tive field theory (EFT), Kaiser [35] did not confirm suchlarge values for e sym , . Instead, Kaiser found e sym , ’ . n sat , which is still about three times largerthan the FFG contribution. Moreover, Kaiser found con-tributions to the energy per particle whose fourth deriva-tive with respect to δ was singular at δ = 0. This was fur-ther substantiated by analytic MBPT calculations basedon an S -wave contact interaction, which gave rise to asingular term ∝ δ log | δ | —a term that only contributesto ANM when δ = 0 and δ = 1, and which will be referredto as the leading-order logarithmic term in the following.Subsequently, Wellenhofer et al. performed a more de-tailed analysis of such divergences by examining the δ dependence of the nuclear EOS as a function of densityand temperature [36]. They found that the asymmetryexpansion is hierarchically ordered, i.e., the lower-ordercoefficients are dominant at finite temperature and lowdensity, but the expansion diverges with alternating signin the zero-temperature limit. Around saturation den-sity, their results indicate that the convergence of the se-ries expansion is restored for T & TABLE I. Nonlocal N LO NN and N LO 3N interactionsused in the MBPT calculations of Ref. [20]. Each Hamil-tonian in the Table is based on the N LO NN potentialEM 500 MeV [37] evolved to the SRG resolution scale λ . Thelow-energy couplings c D and c E were subsequently fit to thetriton binding energy and the charge radius of He in Ref. [38]for different combinations of λ and the 3N cutoff Λ . The3N two-pion exchange is governed by the π N low-energy cou-plings c , c , and c , which were taken from the NN potential,except for H7 which uses the values obtained from the NNpartial-wave analysis (PWA) of Ref. [39]. Hamiltonian H6has been excluded as discussed in Section IV B of Ref. [20].label λ [ fm − ] Λ [ fm − ] 3N c , , c D c E H1 1 . . . − . . . . − . . . − . − . . . . − . . . . − . . . − . − . have argued that the logarithmic term at leading orderconsiderably improves the isospin-asymmetry expansionat zero temperature and suggested to include this termin future fits of the EOS.While mathematically well-defined, it is not clear thatthe aforementioned divergence of the series expansion inthe isospin asymmetry parameter substantially impactsthe practical usability of the expansion (1), because cor-rections remain small at nuclear densities. As remarkedby Wellenhofer et al. , the lack of precision in nonperturba-tive approaches to the nuclear many-body problem makesit practically impossible to precisely extract the high-order isospin-asymmetry derivatives. Furthermore, ourknowledge of the symmetry energy, and even more fun-damentally of the nuclear interaction itself, is limited byexperimental precision and by a limited theoretical un-derstanding of strongly-interacting systems. As a conse-quence, while the series expansion in the isospin asymme-try can be determined with high accuracy when the nu-clear interaction and the many-body treatment are fixed(with numerical limitations as discussed in Ref. [36]), cur-rent theoretical uncertainties reduce our ability to accu-rately determine high-order contributions in general. Inthis paper, we analyze the impact of these uncertaintieson the determination of the symmetry energy. III. NUCLEAR-MATTER EQUATION OFSTATE
The nuclear-matter EOS has been investigated sincethe nineteen-fifties employing various theoretical ap-proaches, see e.g. Refs. [1, 2] and references therein.Recently, improved studies of nuclear matter and finitenuclei [40–42] have become available with the adventof chiral EFT [43, 44] and renormalization group meth-ods [45, 46]. Importantly, chiral EFT interactions have also enabled theoretical uncertainty estimates [22, 47,48]. As for the various many-body approaches, signifi-cant progress has been made in predicting the EOS inthe limits of PNM and SNM using the self-consistentGreen’s function method [49, 50], coupled cluster the-ory [51, 52], MBPT [53–57], in-medium chiral perturba-tion theory [58], and Quantum Monte Carlo (QMC) stud-ies [21, 59–62]. On the other hand, only a few explicitcalculations of ANM exist to date [20, 35, 63, 64], becausetypically ANM is computationally more involved.In this work, we use the explicit ANM calculations ofRef. [20] to study the importance of non-quadratic con-tributions to the symmetry energy. We focus our studyon the zero-temperature limit in which the dominant ef-fects are expected. Drischler et al. evaluated the en-ergy per particle
E/A ( n, δ ) at 11 isospin asymmetries( δ = 0 . , . , . . . , .
0) using MBPT up to second order,and estimated the neglected third-order ladder contribu-tions to be small. For each proton fraction, the energyper particle is sampled on an equidistant grid in the neu-tron Fermi momentum (not the density), and availableup to 2 fm − leading to 385 points in total for each Hamil-tonian. The MBPT calculations of Ref. [20] are based onthe set of six Hamiltonians with the chiral NN and 3Ninteractions summarized in Table I. These chiral Hamil-tonians are also commonly used in nuclear-structure cal-culations [65–74]. They combine the N LO NN potentialEM 500 MeV [37] evolved to lower momentum scales us-ing the similarity renormalization group (SRG) with bareN LO 3N forces regularized by a nonlocal regulator withmomentum cutoff Λ . Hebeler et al. then refit the two3N low-energy couplings c D and c E for different combina-tions of λ and Λ shown in Table I to the triton bindingenergy as well as the charge radius of He [38]. AssumingN LO 3N forces provide a sufficiently complete operatorbasis, and the long-range low-energy couplings c , c , and c are SRG-invariant, this approach captures dominantcontributions from induced three- and higher-body forcesdue to the SRG transformation. Note that the c i ’s ap-pear both, in the NN and 3N interactions at N LO. Asdiscussed in Ref. [20], the spread of energies per particleobtained from these nuclear interactions can serve as asimple uncertainty estimate—although not allowing fora statistical interpretation.Normal-ordering is the standard approach for includ-ing (dominant) 3N contributions in many-body calcula-tions in terms of density-dependent effective two-bodypotentials. In infinite nuclear matter, normal orderingcorresponds to summing one particle in the 3N forcesover the occupied states of the (e.g., Hartree-Fock) ref-erence state. This induces a dependence on the totalmomentum P of the two remaining particles, in contrastto the Galilean invariant NN potentials. To straightfor-wardly combine NN and effective potentials, the approx-imation P = 0 had previously been imposed [53, 75].Reference [20] relaxed this approximation by averagingover all directions of P using a novel partial-wave-basedmethod that generalizes normal-ordering of arbitrary 3N n (fm ) e P N M ( M e V ) APR (1998)Wlazlowski et al . (2014)Tews et al . (2016)Drischler et al . (2016)fit (68% CL)fit (95% CL) k F (fm ) e P N M / e FF G P N M FIG. 1. Comparison of the MBPT predictions for the en-ergy per particle in PNM [20] (blue points) with the APREOS [76] (green squares), QMC calculations Wlazlowski etal. (2014) [30] (cyan triangles), and Tews et al. (2016) [61](black dots) using different chiral EFT Hamiltonians with NNand 3N forces. The latter points include simple estimates forthe EFT truncation error of the chiral expansion. We alsoshow our fit posterior at 68% (95%) confidence level as dark(light) red bands. See the main text for details. interactions. We focus the discussion here on the resultsobtained with this improved (angle averaging) approxi-mation in a Hartree-Fock single-particle spectrum.
A. Energy per particle
The energies per particle in PNM predicted by theMBPT approach [20] are shown in the upper panel ofFig. 1 as blue dots (the lower panel represents the ratioof the energy per particle over the free Fermi gas energy).For the fits to the MBPT data that we use later in ouranalysis, we also show 68% (95%) confidence intervals ofthe fit posteriors as dark (light) red bands. In addition,we compare with the variational calculation of Ref. [76](APR), Fock-space formulated Quantum Monte Carlo(QMC) calculations of Ref. [30] (Wlazlowski et al. et al. LO NN interactions of Ref. [79] combined with N LO3N forces as specified in Ref. [80], and Ref. [61] useslocal coordinate-space chiral interactions constructed inRefs. [59, 60, 62]. The first two calculations do not pro-vide theoretical uncertainties, while the latter estimatesthe standard EFT uncertainty [47]. Note that, in general,order-by-order calculations are required for estimatingEFT truncation errors. Such calculations are not pos-sible with the chiral Hamiltonians given in Table I.When comparing the approaches using chiral EFT in-teractions, the QMC calculations of Ref. [61] agree withthe MBPT approach employed in this work within uncer-tainties above n ∼ .
08 fm − , while QMC finds slightlyhigher energies at lower densities. In contrast, the QMCcalculations of Ref. [30] find a higher PNM energy perparticle at all densities, by about ∼ e PNM /e FFGPNM as a function of neutronFermi momentum k F for the various calculations in thebottom panel of Fig. 1. In the figure, we can identify thedensity region where the ratio exhibits a plateau, indi-cating a similar scaling of e PNM and e FFG with k F . Forthe MBPT calculation, we find the ratio at the plateauto be ≈ . k F ≈ . − ,which describes densities of ≈ n sat / B. Single-particle energies
The single particle energy (cid:15) τ ( k ) is determined by theself-consistent solution to the Dyson equation. In aHartree-Fock spectrum it is given by (cid:15) τ ( k, n, δ ) ≈ k m τ + Σ (1) ( k, n, δ ) . (7) k (fm ) n ( M e V ) SNMPNMNN-only n = n empsat k F , SNM k F , PNM k (fm ) n ( M e V ) SNMPNM k F , SNM k F , PNM NN+3N n = n empsat H1H2H3H4H5H7
FIG. 2. The neutron single-particle energies (cid:15) n ( k ) as a function of the momentum k calculated at n empsat , and extracted fromthe MBPT calculations of Ref. [20]. The different colors correspond to the six Hamiltonians as labeled in the legend. We showthe single particle energies obtained from only NN forces (left panel) and when including 3N contributions (right panel). Ineach panel, we present results for both SNM and PNM. The first term in Eq. (7) is the single-particle kinetic en-ergy, while the second term Σ (1) denotes the spin-isospin-averaged first-order self-energy. We refer the reader to,e.g., Refs. [20, 53] for more details.Figure 2 shows the single-particle energy (cid:15) n ( k ) in SNMand PNM evaluated at n empsat . The left (right) panel de-picts the NN-only (NN+3N) results, and the verticallines mark the position of the neutron Fermi momen-tum in SNM ( k F, SNM = 1 .
33 fm − ) and PNM ( k F, PNM =1 .
68 fm − ) associated with the nuclear saturation density, n empsat . The different curves show the results for the sixHamiltonians H1 to H7 specified in Table I. The spread islarger in SNM (about 15 MeV) compared to PNM (about5 MeV) because the 3N short- and intermediate-rangecontributions governed by c D and c E do not contributeto the PNM EOS for nonlocal regulator functions. Asexpected, SNM is more attractive than PNM, as a resultof the attractive contributions from the T = 0 channels,which are absent in PNM. C. Landau mass
The momentum dependence of the nuclear interactionscan be absorbed by a modification of the mass of the nu-cleons, which gives rise to the so-called in-medium effec-tive mass and the Landau mass. Specifically, the single-particle energy can be approximated as, (cid:15) τ ( k, n, δ ) ≈ k m ∗ τ ( k, n, δ ) + Σ (1) ( k = 0 , n, δ ) , (8)where the in-medium effective mass is defined as [81], m ∗ τ ( k, n, δ ) m τ = km τ (cid:18) d(cid:15) τ ( k, n, δ ) dk (cid:19) − . (9) Finally, the Landau mass is defined as the effectivemass (9) taken at k = k F .The effective masses in SNM and PNM are shown inFig. 3 as functions of the momentum k at a fixed density n empsat . The effective masses are lower in SNM comparedto PNM, in agreement with BHF calculations [31, 82, 83].We find that the inclusion of 3N forces leads to several in-teresting effects on the effective mass: (a) 3N forces gen-erate a stronger momentum dependence compared withNN-only calculations, and (b) 3N forces have a largerimpact on PNM than on SNM. Furthermore, the disper-sion among the different Hamiltonians is slightly largerwhen 3N forces are included. From Fig. 3, we find forthe Landau mass m ∗ n /m ( δ = 0) = 0 . m ∗ n /m ( δ =1) = 0 . Dm ∗ n, sat = m ∗ n ( n sat , δ = 1) − m ∗ n ( n sat , δ = 0) , (10)is about Dm ∗ n, sat = 0 . n sat .In Fig. 4, we show the Landau mass (left) and its in-verse (right) considering NN-only forces (dashed lines)and NN and 3N forces (gray bands) in SNM and PNMas a function of the density, n . The difference of the Lan-dau masses in PNM and SNM, Dm ∗ n ( n ), increases withdensity, and is found to be about 0 .
24 at saturation den-sity, see also Fig. 3. While it is usually found that theLandau mass decreases with density [31, 82, 83], we findthat in PNM the Landau mass first decreases at lowerdensity, but increases again for n > . − (except forHamiltonian H3, which has a higher momentum cutoffapplied to the 3N forces). This effect is due to the inclu-sion of 3N interactions in the Hamiltonian.Because many energy density functional (EDF) ap-proaches approximate the inverse of the Landau massby a linear function in density [29], we show the inverse k (fm ) m * n / m SNMPNM k F , SNM k F , PNM NN-only n = n empsat k (fm ) m * n / m SNMPNM k F , SNM k F , PNM NN+3N n = n empsat H1H2H3H4H5H7
FIG. 3. Same as Fig. 2 but for the neutron effective mass as function of the momentum k . n (fm ) m * n / m SNMPNM k = k F NN-onlyNN+3N n (fm ) m / m * n PNMSNM k = k F NN-onlyNN+3N
H1H2H3H4H5H7
FIG. 4. Landau effective mass (left) and its inverse (right) in SNM and PNM as a function of the density. The black-dashedlines represent the upper and lower limits when only NN forces are considered, while the grey-shaded regions show the resultswith 3N forces included. The different colors correspond to the six Hamiltonians as labeled in the legend.
Landau mass in the right panel of Figure 4. In contrastto the EDFs approaches, we find that the density depen-dence of the inverse Landau mass is not linear, and that3N forces enhance the nonlinear behavior. To describethe inverse of the Landau mass as function of the density n and asymmetry parameter δ , we consider the followingfunctional form: (cid:18) m ∗ τ m ( n, δ ) (cid:19) − = 1 + (cid:18) κ sat n sat + τ δ κ sym n sat (cid:19) n + (cid:18) κ sat , n + τ δ κ sym , n (cid:19) n , (11) where τ = 1( −
1) for neutrons (protons). Using Eq. (11),we obtain in SNM and PNM, respectively, (cid:18) m ∗ SNM m ( n ) (cid:19) − = 1 + κ sat n sat n + κ sat , n n , (12) (cid:18) m ∗ PNM m ( n ) (cid:19) − = 1 + κ PNM n sat n + κ PNM , n n , (13)with κ PNM = κ sat + κ sym and κ PNM , = κ sat , + κ sym , .The parameters κ sat , κ sat , , κ sym , and κ sym , are ob-tained from fitting the expressions (12) and (13) to theresults shown in the right panel of Figure 4. The detailsof our parametric fits is discussed in Appendix A. Therelevant fit parameters, p α , are p α = { κ sat /n sat , κ sat , /n , κ PNM /n sat , κ PNM , /n } . (14)These fit parameters are determined from the predictedLandau effective masses for each of the six Hamiltoni- TABLE II. Fit parameters of the inverse Landau mass con-sidering linear and quadratic density expansions. The fitsare compared to three Skyrme-type interactions: NRAPR [3],LNS5 [84], and SAMI [85]. κ sat /n sat κ sat , /n κ PNM /n sat κ PNM , /n [fm ] [fm ] [fm ] [fm ]linear 3 . . . − . . − . .
75 - 1 .
40 -LNS5 [84] 4 .
12 - 2 .
19 -SAMI [85] 3 .
03 - 2 .
87 - ans. The results of the fits for the inverse of the Landaumass are given in Table II, where we have consideredboth, a linear and a quadratic fit function. The priordistribution for each of the fit parameters is given by anormal distribution with mean 0 and standard deviation100, providing an uninformative prior. The fits are com-pared to three Skyrme-type interactions: NRAPR [3],LNS5 [84], and SAMI [85] that satisfy the following con-ditions: 0 . m ∗ /m (SNM) .
7, ∆ m ∗ /m > < L sym <
60 MeV.In Fig. 5, we compare the posterior distribution func-tions for the Landau mass in SNM (top panel) and PNM(bottom panel), and the input data. The predictionsfrom the six Hamiltonians are plotted as solid lines, and,at each density, we calculate the centroid and 1 σ intervalgiven the six Hamiltonians (black points with error bars).The ± σ ( ± σ ) contours of the posterior, correspond-ing to the 68% (95%) confidence region, are depicted inred (light red) for the quadratic fit and dark blue (lightblue) for the linear fit. We fit the models to the data inthe range n = 0 . − .
17 fm − for the linear fit (3 datapoints) and from n = 0 . − .
20 fm − for the quadraticfit (14 data points). These values are chosen to allow forthe ranges to be as large as possible while, at the sametime, ensuring that the fits reproduce the data aroundsaturation density. While the quadratic fit performs welleven outside the fit interval, down to n ∼ .
05 fm − inSNM and PNM, the linear fit does not because of thestrong curvature of the Landau mass. The differencesbetween the linear and quadratic fits are further anal-ysed in Sec. VI.Finally, we study the splitting of the neutron and pro-ton Landau masses in ANM, defined as∆ m ∗ sat ( δ ) = m ∗ n ( n sat , δ ) − m ∗ p ( n sat , δ ) . (15) n (fm ) m * n / m H1H2H3H4H5H7
Symmetric matter quadratic fit (68% CL)quadratic fit (95% CL)linear fit (68% CL)linear fit (95% CL)data (68% CL)0.00 0.05 0.10 0.15 0.20 n (fm ) m * n / m H1H2H3H4H5H7
Neutron matter
FIG. 5. Results of the Bayesian parametric fits of the in-verse Landau mass. The 68% (95%) confidence levels forthe posterior distribution functions are shown as dark-shaded(light-shaded) bands. The black lines represent the individualHamiltonians, and the black points show the average over thesix Hamiltonians with ± σ uncertainty bands. In PNM ( δ = 1), this splitting can be expressed in termsof the difference Dm ∗ sat , see Eq. (10), as∆ m ∗ sat m (PNM) ≈ Dm ∗ n, sat m + O (cid:18) κ sym + κ sym , κ sat + κ sat, (cid:19) ! . (16)From our fits, we can estimate that the neglected termsaccount for about 5% of the splitting (more precisely,7% for the linear fit of the effective mass, and 3% forthe quadratic fit), which is small considering the presentuncertainty of this quantity. The splitting of the Landaumass is, thus, approximately given by the difference ofthe Landau mass between PNM and SNM. The splittingof the Landau mass obtained here is compatible with theone obtained in the literature for BHF [31, 86, 87] andDirac-BHF [82, 83] approaches. IV. ENERGY EXPANSION IN THE ISOSPINASYMMETRY PARAMETER δ A general expression for the expansion (1) of energyobservables in nuclear matter was suggested in Eq. (27)of Ref. [36]. From this expression, we consider all con-tributions up to δ , including the logarithmic term, andrewrite it as y ( n, δ ) ≈ y SNM ( n ) + y sym , ( n ) δ + y sym , ( n ) δ + y sym,log ( n ) δ log | δ | . (17)The term y sym , log , δ log | δ | originates from the second-order contribution in the many-body expansion in theneutron-proton channel, as explained in Ref [36]. Thecorresponding contribution to the symmetry energy y sym is defined as y sym ( n ) = y PNM ( n ) − y SNM ( n ) , (18)where y SNM ( n ) = y ( n, δ = 0) and y PNM ( n ) = y ( n, δ = 1).The non-quadratic contribution to the symmetry energyis defined as y sym,nq ( n ) = y sym ( n ) − y sym,2 ( n ) . (19)Note, that since the logarithmic term vanishes in SNMand PNM, it also does not contribute to the non-quadratic term (19).The quantity y in Eq. (17) can be the energy per parti-cle e , as originally suggested by Kaiser [35], or any otherenergy contribution can be expanded as well. For in-stance, as discussed in Sec. I, one can expand the kineticenergy t , which coincides with the FFG energy given inEq. (4), t = e FFG . The δ -expansion of the kinetic energyup to fourth order is given by t ( n, δ ) ≈ t SNM ( n ) + t sym , ( n ) δ + t sym , ( n ) δ , (20)with t SNM ( n ) = t satSNM (cid:18) nn sat (cid:19) / , (21) t sym , ( n ) = 59 t satSNM (cid:18) nn sat (cid:19) / , (22) t sym , ( n ) = 5243 t satSNM (cid:18) nn sat (cid:19) / . (23)We note that t sym, log = 0 for the kinetic energy becauseit is a one-body operator. At saturation density, one finds t sym , ( n sat ) ≈ .
45 MeV, while the contribution of allhigher-order terms starting with t sym,6 sum up to about0.25 MeV.When including nuclear interactions, the total energyis given by the sum of the kinetic energy t ( n, δ ) and thepotential energy e pot ( n, δ ). The potential energy can,thus, be defined as e pot ( n, δ ) = e ( n, δ ) − t ( n, δ ) , (24) where the total energy e ( n, δ ) is obtained from theMBPT calculations [20] and the kinetic energy is ana-lytic. Again, the potential energy can be expanded usingEq. (17), e pot ( n, δ ) ≈ e potSNM ( n ) + e potsym , ( n ) δ + e potsym , ( n ) δ + e potsym,log , ( n ) δ log | δ | . (25)Another way to express the total energy is to explicitlyinclude the in-medium correction to the nucleon mass,modifying the kinetic and potential energies but not theirsum, e ( n, δ ) = t ∗ ( n, δ ) + e pot ∗ ( n, δ ) , (26)where t ∗ ( n, δ ) is the effective kinetic energy calculatedwith the Landau mass m ∗ τ ( δ ). Here, τ = n or p for neu-trons and protons, respectively, and e pot ∗ ( n, δ ) is the ef-fective potential energy. We refer to Sec. III C for moredetails on the Landau mass. The effective kinetic energyis defined as t ∗ ( n, δ ) = t satSNM (cid:18) nn sat (cid:19) / h mm ∗ n ( δ ) (1 + δ ) / + mm ∗ p ( δ ) (1 − δ ) / i , (27)and the effective potential energy is given by e pot ∗ ( n, δ ) = e ( n, δ ) − t ∗ ( n, δ ) (28)= t ( n, δ ) − t ∗ ( n, δ ) + e pot ( n, δ ) . Similarly to e and e pot , the effective potential energy e pot ∗ can be expanded using Eq. (17), e pot ∗ ( n, δ ) ≈ e pot ∗ SNM ( n ) + e pot ∗ sym , ( n ) δ + e pot ∗ sym , ( n ) δ + e pot ∗ sym,log ( n ) δ log δ . (29)In the following, we use these notations for analyzingthe δ -dependence of the total, potential, and effectivepotential energies. V. META-MODEL FOR SYMMETRIC ANDNEUTRON MATTER
To describe the MBPT data for the energy per parti-cle in SNM and PNM, we use in this work a functionalform described by a meta-model (MM) for nuclear mattersimilar to the one suggested in Ref. [16], but generalizedto a potential energy with non-quadratic δ dependence.The MM is adjusted to MBPT data sampled on a givengrid in the asymmetry parameter δ [20]. This is in con-trast to Wellenhofer et al. [36], who used a finite dif-ference method [88] on an adjustable grid to determineall derivatives with respect to δ of interest. The MM,instead, provides a flexible polynomial-type approach tonuclear matter, which allows us to accurately determinethe higher-order coefficients in the δ expansion, even forthe fixed grid considered here.For SNM and PNM, the energy per particle in the MMreads e MMSNM ( n ) = t ∗ SNM ( n ) + e pot ∗ SNM ( n ) , (30) e MMPNM ( n ) = t ∗ PNM ( n ) + e pot ∗ PNM ( n ) . (31)The kinetic energy is determined from Eq. (27) with theLandau mass, see Sec. III C. The potential energies areexpanded about n sat in terms of the parameter x ≡ n − n sat n sat as follows e pot ∗ SNM ( n ) = N X j =0 j ! v SNM ,j x j + v low − nSNM x N +1 e − b sat nn empsat ,e pot ∗ PNM ( n ) = N X j =0 j ! v PNM ,j x j + v low − nPNM x N +1 e − b PNM nn empsat , where b PNM = b sat + b sym , and the second term on theright is a low-density correction. This correction rep-resents the low-density contribution of all higher-orderterms neglected in the summation, and scales like x N +1 at leading order, where N is the upper limit of the powerin the density expansion. In the original nucleonic MM ofRef. [16], the low-density EOS correction was simply pa-rameterized by a fixed coefficient b = b sat = b PNM ≈ . b sat and b sym ) controlling the density depen-dence of the low-density corrections in PNM and SNMseparately. It was suggested that considering an expan-sion up to N = 4 allows the reproduction of the pressureand sound speed of about 50 known EDF up to about4 n sat , see Ref. [16] for more details. In principle, it isnot necessary to consider such a high N in the presentanalysis. The inclusion of high-order contributions, how-ever, affects the determination of the low-order ones, asdiscussed in Ref. [89], even if the data does not constrainthe high-order terms themselves.Imposing that the energies per particle vanish at n = 0 fm − , we obtain the following relations e pot ∗ SNM ( n ) = N X j =0 j ! v SNM ,j x j u SNM ,j ( x ) , (32) e pot ∗ PNM ( n ) = N X j =0 j ! v PNM ,j x j u PNM ,j ( x ) , (33)where u α,j ( x ) = 1 − ( − x ) N +1 − j e − b α n/n empsat , (34)and α indicates either SNM or PNM and the correspond-ing b sat or b PNM .In the MM, the coefficients v α, to v α,N are related tothe nuclear empirical parameters (NEPs), such as E sat , K sat , E sym , L sym , etc. The NEPs for SNM are definedby the density expansion e SNM ( n ) = E sat + K sat x + Q sat x + Z sat x + . . . , (35)whereas the NEPs for PNM are defined by e PNM ( n ) = E PNM + L PNM x + K PNM x + Q PNM x + Z PNM x + . . . . (36)We use the following relations between the MM param-eters and the NEPs for the isoscalar parameters control-ling the SNM EOS: v SNM , = E sat − t SNM (1 + κ sat + κ sat , ) , (37) v SNM , = − t SNM (2 + 5 κ sat + 8 κ sat , ) ,v SNM , = K sat − t SNM ( − κ sat + 20 κ sat , ) ,v SNM , = Q sat − t SNM (4 − κ sat + 40 κ sat , ) ,v SNM , = Z sat − t SNM ( − κ sat − κ sat , ) , while the isovector parameters describing the PNM EOSare v PNM , = E PNM − t SNM (1 + κ PNM + κ PNM , ) , (38) v PNM , = L PNM − t SNM (2 + 5 κ PNM + 8 κ PNM , ) ,v PNM , = K PNM − t SNM ( − κ PNM + 20 κ PNM , ) ,v PNM , = Q PNM − t SNM (4 − κ PNM + 40 κ PNM , ) ,v PNM , = Z PNM − t SNM ( − κ PNM − κ PNM , ) . These relations represent another difference to the orig-inal nucleonic MM of Ref. [16], where the isovector co-efficients were determined assuming a quadratic isospin-asymmetry dependence of the symmetry energy. Theisovector contribution of the present MM is built onthe global symmetry energy (3), which allows for possi-ble non-quadratic contributions to the symmetry energy.These contributions will be estimated from the differencebetween the global symmetry energy and its quadraticcontribution, as detailed in Sec. VI.In the MM used here, there are five NEP in SNM,including n sat , and five additional NEP in PNM. Con-sidering the two parameters controlling the low-densityEOS, b sat and b sym , there is a total of 12 parametersthat need to be determined. Note, that these parame-ters carry uncertainties that reflect the current lack ofknowledge of the nuclear EOS. In our Bayesian fits, weuse the priors for the NEP from the analysis presentedin Refs. [16, 90] and summarized in Table III. Here, weadditionally vary the parameters b sat and b PNM to repro-duce the low-density behavior of the energy per particlein SNM and PNM.We show the fitted parameters in Table III. Note, thatboth the posteriors and priors given in the table repre-sent the means of normal distributions with the standarddeviations given in parenthesis. The posteriors we obtainfor the NEPs may depend on the exact representation of0
TABLE III. Priors and posteriors of the NEP from analyses of SNM and PNM. We report results for the different scalingsdescribed in the text. Values within parentheses represent the error bars at the ± σ level. NEPs for the following threeSkyrme-type interactions are given: NRAPR [3], LNS5 [84] and SAMI [85].Scaling n sat E sat K sat Q sat Z sat b sat (fm − ) (MeV) (MeV) (MeV) (MeV)prior 0 . − . − . − . − . − . − . − . − . − . − † NRAPR [3] 0 . − .
85 226 −
363 1611LNS5 [84] 0 . − .
57 240 −
316 1255SAMI [85] 0 . − .
93 245 −
339 1330Scaling n sat E PNM L PNM K PNM Q PNM Z PNM b PNM (fm − ) (MeV) (MeV) (MeV) (MeV) (MeV)prior - 16 . − . †† . − − . †† . − . †† . − . †† . − − † NRAPR [3] 0 .
161 18 .
33 65 108 − − .
160 15 .
29 57 130 − − .
159 13 .
32 47 127 35 − † Fixed parameter. †† Quoted values are the n sat priors considered in PNM and obtained from SNM posteriors. the data points, i.e., if the data is equidistant in n or k F . To gauge the sensitivity to this choice, we investi-gate in the following three possible data representations,here called scalings. We show the results for all scalingsin Fig. 6.Since we require our fit to provide a faithful represen-tation of low-density nuclear matter (with a fair weightfor the low-density points) in order to fix b sat and b PNM ,we adopt for scaling 1 the representation of e/e
FFG as afunction of the Fermi momentum k F . Scaling 1 providesthe best representation to analyze the low-density prop-erties of the energy per particle because an equidistantgrid in k F leads to a very dense data set at low densities.Note that, as mentioned in Sec. III, the original MBPTdata [20] is provided on an equidistant grid in k F . Thescaling of the y -axis normalizes the energies to the sameorder of magnitude at all k F . For scaling 2, we choose the representation of e/e FFG on an equidistant grid indensity. We use a cubic spline to interpolate the energiesper particle from the original Fermi momentum grid toa equally spaced density grid. By switching the equidis-tant grid in momenta to densities, scaling 2 reduces theweight for the low-density data points and, therefore, ismore appropriate to fit the NEPs, which are determinedaround saturation density. Finally, scaling 3 representsthe energy per particle on an equidistant grid in density,as it is more often presented in the literature. Hence, theonly difference to scaling 2 is the scaling of the y -axis.The results for each of the three scalings are shown inFig. 6 for SNM and PNM, while the posteriors for theNEPs are given in Table III. Note, that the NEP n sat isonly meaningful in SNM, while its uncertainty influencesthe determination of the NEPs in PNM. In our approach,we therefore vary n sat in PNM within the posterior un-1 e S N M / e FF G S N M Scaling 1, SNM e S N M / e FF G S N M Scaling 2, SNM e S N M ( M e V ) Scaling 3, SNM k F (fm ) e P N M / e FF G P N M Scaling 1, PNM fit (68% CL)fit (95% CL)data (68% CL) n (fm ) e P N M / e FF G P N M Scaling 2, PNM n (fm ) e P N M ( M e V ) Scaling 3, PNM
FIG. 6. Comparison of the Bayesian inference results for the MM of this work (red bands) with the MBPT data (blue points) forSNM (top panels) and PNM (bottom panels) and the three different scalings described in the main text. The bands are givenat 65% (dark-shaded) and 95% confidence level (light-shaded), whereas the data points are shown with the ± σ uncertaintyestimate. certainty obtained from the fit in SNM. In this way, theNEP in PNM naturally include the uncertainty in n sat .Scaling 1 is ideal for analyzing the low-density behaviorof the energy per particle. When simultaneously vary-ing the 12 MM parameters, we find b SNM = 17(5) and b PNM = 42(4) as well as the values for the 10 NEP givenin Table III. The density dependence of the low-densitycorrection is, thus, very different in SNM and PNM, incontrast to the original MM of Ref. [16]. Illustratingthe importance of the scaling choice, we find some dif-ferences between the NEPs obtained from scaling 1 andscalings 2 and 3. These differences are usually small com-pared to the uncertainties, except for Q sat in SNM, aswell as Q PNM in PNM. We note that the fits from scal-ings 2 and 3 are identical and, hence, the scaling of theenergy with respect to e FFG has a negligible effect.Finally, we fix the values for b SNM and b PNM from scal-ing 1, and re-fit all remaining NEP considering scaling 3.The results are referred to as scaling 3 ∗ . Fixing b SNM and b PNM has the largest impact on Q sat and Q PNM , asexpected, but differences between scaling 3 and 3 ∗ aresmall compared to the overall uncertainties. We, there-fore, conclude that the parameters b SNM and b PNM donot have a significant impact on the determination of theNEPs, and can be fixed from the fit to low-density mat-ter (scaling 1). We stress that the higher-order NEPs Z sat and Z PNM are not constrained by our data, becausethe density range of the MBPT data is limited to den- sities below 0 . − . However, they contribute to theuncertainty of the other NEPs [89].For the NEPs describing nuclear saturation we ob-tain n sat = 0 . − , E sat = − . K sat = 226(18) MeV. The results are consistentwith the original analysis in Ref. [20], which obtained n sat = 0 . − .
190 fm − , E sat = − (15 . − .
3) MeV,and K sat = 223 −
254 MeV using a Hartree-Fock spec-trum. However, our uncertainties are generally smallerbecause we explicitly guide our fits in Fig. 6 by empirical(or “expert”) knowledge [16] through prior distributionsof the fit parameters. In PNM, where empirical con-straints are lacking, the fits are therefore closer to theMBPT data.
VI. SYMMETRY ENERGY
Using the results obtained in Secs. III C and V, we nowdetermine the properties of the symmetry energy and therelative contributions of the quadratic and non-quadraticterms.
A. Global symmetry energy e sym The global symmetry energy e sym is simply determinedfrom the fits of PNM and SNM, see Eq. (3) and Sec. V,2 n (fm )05101520253035 e s y m ( M e V ) n (fm )05101520253035 e p o t s y m ( M e V ) model (68% CL)model (95% CL)data (68% CL) 0.00 0.05 0.10 0.15 0.20 n (fm )05101520253035 e p o t * s y m ( M e V ) data (Quad. fit)Quadratic fit (68% CL) FIG. 7. Results for the symmetry energy, e sym ( n ), its potential contribution e potsym , and the effective potential e pot ∗ sym usingEqs. (41), (39), and (40). The meaning of the individual bands and points is the same as in Fig. 6. In the right panel, the light-and dark-shaded red bands and the data (blue points) correspond to calculations where the Landau mass is represented by alinear polynomial. The black squares (without error bars) and the black dashed lines (enclosing a band) represent calculationswhere the Landau mass is represented by a quadratic fit. while the contributions of the potential energies e potsym , and e pot ∗ sym are obtained from e sym according to the followingdefinitions, e potsym ( n ) = e sym ( n ) − t PNM ( n ) + t SNM ( n ) , (39) e pot ∗ sym ( n ) = e sym ( n ) − t ∗ PNM ( n ) + t ∗ SNM ( n ) . (40)We use the fits of the Landau mass, see Sec. III C, todetermine t ∗ , including its uncertainties, as explained inthe following.We present these quantities in Fig. 7 as functions of thedensity n . For e sym , the data points are obtained fromthe PNM and SNM data, while the error bars are definedas the arithmetic average of the PNM and SNM errorbars. For the model, we employ the symmetry energydetermined from the MM, which is defined as e MMsym ( n ) = e MMPNM ( n ) − e MMSNM ( n ) . (41)The results shown in Fig. 7 are obtained from the bestfits to SNM and PNM (scaling 3 ∗ in Table III), wherethe width of the bands is defined as the arithmetic av-erage of the widths in SNM and PNM. The model re-sults, therefore, depend on the choice of prior in SNMand PNM, in particular, on the prior knowledge of thesaturation density and energy considered in SNM, seethe discussion of Fig. 6. For this reason, the MM un-certainty for the symmetry energy is slightly smallerthan the uncertainty of the data in Fig. 7. At nuclearsaturation density, n sat = 0 . − , the data sug-gest e sym = 30 . e sym = 31 . . ± . LO, 30 . E ). Similarly, we predict L sym = 46 . L sym = 58 . L sym = 58 . L sym , however,is very sensitive to the densities at which the value is ex-tracted as well as to the interactions employed, see, e.g.,Ref. [21].The data for e potsym and e pot ∗ sym are obtained from e sym using Eqs. (39) and (40). In the case of e pot ∗ sym , the ef-fective mass is fixed to be the best fit using either thelinear or the quadratic density expansion (depicted bydashed lines in the right panel), and the uncertainty of e pot ∗ sym is defined as the arithmetic average of the uncer-tainties of e sym , t ∗ PNM and t ∗ SNM . Therefore, the un-certainty of e pot ∗ sym also includes the uncertainties in theLandau mass parameters κ sat and κ PNM . We observethat there is a large impact of the Landau mass on e pot ∗ sym ,compared to e potsym with the bare mass. At n sat , we ob-tain e potsym = 18 . e pot ∗ sym = 25 . −
40% as discussed inthe introduction. These numbers are compatible with theexpectations for the complementary contribution fromthe kinetic energy. We find t PNM ( n empsat ) − t SNM ( n empsat ) =13 . t ∗ PNM ( n empsat ) − t ∗ SNM ( n empsat ) = 5 . e pot ∗ sym , we expect a difference when using either aLandau mass that is linear or quadratic in density, seeFig. 5. In Fig. 7 we show two results for e pot ∗ sym . Theblue data points and the dark-shaded (light-shaded) redbands correspond to the results at 68% (95%) confidencelevel when using a Landau mass linear in density. The3black squares and the black-dashed lines, encompassingthe corresponding 68% confidence interval, represent cal-culations with a Landau mass quadratic in density. In-terestingly, the values for e pot ∗ sym obtained from these twofunctions for the Landau mass differ only by about 1 . B. Quadratic contribution to the symmetry energy
The quadratic contribution to the symmetry energy, e sym , , is defined in Eq. (2) as the local curvature inthe isospin-asymmetry parameter δ in SNM. In the fol-lowing, we extract e sym , using this expansion aroundSNM, but also suggest obtaining e sym , from an expan-sion around PNM. We demonstrate that both definitionsprovide comparable results.
1. Expansion around SNM
The quadratic contribution to the symmetry energy isdefined by Eq. (2) relative to the SNM EOS. To deter-mine this contribution directly from the MBPT data, weemploy Eq. (1) up to the fourth order in δ , and fit thecoefficients e sym , and e sym , using optimize.curve fit fromthe Python package SciPy . The fits are performed forisospin asymmetries in the range, δ = 0 . − .
5. We havechecked that the results are insensitive (within variationsof about 0 .
1) to the upper limit of this range—as long asit is chosen to be > . e sym , ( n ) as a function ofthe density using the MM contribution to the quadraticsymmetry energy, e MMsym , ( x ) = 59 t SNM ( x ) + N X j =0 x j j ! (cid:20) v sym2 ,j u j ( x, δ = 0) − v SNM,j (cid:0) u j ( x, δ = 0) − (cid:1) (1 + 3 x ) b sym (cid:21) , (42)where the parameters v sym2 ,i are determined using Lep-age’s Python package lsqfit , as before for other quanti-ties. They are related to the quadratic symmetry energyNEPs as follows: v sym2 , = E sym2 − t satSNM [1 + κ sat + 3 κ sym + κ sat,2 ] , (43) v sym2 , = L sym2 − t satSNM [2 + 5( κ sat + 3 κ sym ) + 8 κ sat,2 ] ,v sym2 , = K sym2 − t satSNM [ − κ sat + 3 κ sym ) + 20 κ sat,2 ] ,v sym2 , = Q sym2 − t satSNM [4 − κ sat + 3 κ sym ) + 40 κ sat,2 ] ,v sym2 , = Z sym2 − t satSNM [ − κ sat + 3 κ sym ) − κ sat,2 ] . These relations generalize the ones in Ref. [16] for aquadratic density-dependent Landau mass. In Eq. (42),the parameter b sym ≡ b PNM − b SNM is fixed by the 3 ∗ fit.The contributions to the symmetry energy due to thepotential energy are determined from the following ex-pressions, e potsym , ( n ) = e sym , ( n ) − t SNM ( n ) , (44) e pot ∗ sym , ( n ) = e sym , ( n ) − t SNM ( n ) (cid:20) κ sat nn sat + κ sat , (cid:18) nn sat (cid:19) + 3 κ sym nn sat (cid:21) . (45)Our results for e sym , , e potsym , , and e pot ∗ sym , are shown inthe first row of Fig. 8. At n empsat , we find e sym , ( n empsat ) =30 . e potsym , ( n empsat ) = 17 . e pot ∗ sym , ( n empsat ) = 26 . .
7) MeV (with the linear density-dependent model for the Landau mass). The large valueof e pot ∗ sym , ( n empsat ), almost 90% of e sym , ( n empsat ) originatesfrom the isospin dependence of the Landau mass, en-coded by κ sym .From the fit model (42), we obtain an estimate for theNEPs that govern the quadratic contribution to the sym-metry energy at the inferred value of n sat . The values aregiven in the first row of Table IV. Our result for E sym , isabout 1 MeV below the total symmetry energy, E sym —the difference is due to non-quadratic contributions.4 e s y m , ( M e V ) Delta expansion 102030 e p o t s y m , ( M e V ) Delta expansion 102030 e p o t * s y m , ( M e V ) Delta expansion fit (68% CL)fit (95% CL)data (68% CL) e P N M s y m , ( M e V ) Eta expansion 102030 e p o t , P N M s y m , ( M e V ) Eta expansion 102030 e p o t * , P N M s y m , ( M e V ) Eta expansion0.00 0.05 0.10 0.15 0.20 n (fm )202 D i ff e r e n c e ( M e V ) n (fm )202 D i ff e r e n c e ( M e V ) n (fm )202 D i ff e r e n c e ( M e V ) FIG. 8. Comparison of the extracted e sym , , e potsym , , and e pot ∗ sym , using Eqs. (42), (44), and (45) (first row) or via an expansionaround PNM, i.e., using Eq. (50) (second row). The third row shows the difference between the δ and η expansions.TABLE IV. Posteriors of empirical parameters obtained from the analysis of the δ and η expansions for e sym , . Valuesinside parentheses represent error bars at the ± σ level. Results for the following three Skyrme-type interactions are given:NRAPR [3], LNS5 [84] and SAMI [85].Expansion n sat E sym , L sym , K sym , Q sym , Z sym , b sym (fm − ) (MeV) (MeV) (MeV) (MeV) (MeV)Prior - 31 . − − − δ (SNM) 0 . † . − − †† η (PNM) 0 . † . − − †† NRAPR [3] 0 .
161 32 .
78 60 −
123 312 − .
160 29 .
15 51 −
119 286 − .
159 28 .
16 44 −
120 372 − † Priors taken from the SNM posteriors in Table III. †† Fixed value.
2. Expansion around PNM
An alternative approach is to determine the contribu-tion e sym , from an expansion around PNM. Since theMBPT approach used here is able to explore asymmetricmatter with arbitrary δ , we can test the accuracy of thisalternative expansion.To this end, we introduce the parameter η = 1 − δ = 2 n p /n , (46) which is twice the proton fraction. Equation (1) can nowbe re-expressed in terms of the this parameter, e ( η ) = e PNM − e sym , + 2 e sym , ) η + ( e sym , + 6 e sym , ) η − e sym , η + e sym , η + O ( η ) . (47)From Eq. (47), it follows then e PNMsym , ( n ) = − ∂e∂η (cid:12)(cid:12)(cid:12)(cid:12) η =0 − ∂ e∂η (cid:12)(cid:12)(cid:12)(cid:12) η =0 . (48)5We determine e PNMsym , and e PNMsym , from fitting the func-tion e ( n, η ) = e PNM ( n ) + a ( n ) η + a ( n ) η + O ( η ) . (49)with e PNMsym , ( n ) = −
14 [3 a ( n ) + 2 a ( n )] , (50) e PNMsym , ( n ) = 18 [ a ( n ) + 2 a ( n )] , (51)at each density to the computed energies per particleat η = 0 . , . , . , and 0 .
3. Again, we also perform aBayesian fit using Eq. (42). The two quantities are shownin the second row of Fig. 8, together with the potentialterms e pot , PNMsym , and e pot ∗ , PNMsym , . The NEPs obtained fromEq. (42) are given in the second row of Table IV. Note,that the differences between the NEPs extracted aroundSNM [using Eq. (2)] and around PNM [using Eq. (48)] aremuch smaller than the uncertainties of these NEPs, whichdemonstrates that the two approaches are consistent withone another. This is further illustrated in the third rowof Fig. 8, where the difference between e SNMsym , and e PNMsym , is shown to be consistent with zero and a small width ofabout 1.5 MeV at n sat .Determining the quadratic contribution to the symme-try energy from an expansion around PNM might be ben-eficial because PNM can usually be computed with muchhigher accuracy since the uncertainties in the 3N inter-actions are reduced. Furthermore, such an extraction isuseful for microscopic approaches in which a small pro-ton impurity can be treated more easily than SNM, e.g.,the auxiliary-field diffusion Monte Carlo approach [21]. C. Non-quadratic contribution to the symmetryenergy e sym ,nq and e sym , We now evaluate the contribution of the non-quadraticterms, defined by Eq. (19), using the expansions aroundSNM and PNM, respectively. For the global symme-try energy, we use our model (41), while for describ-ing the quadratic contribution we use the fit (42). Fig-ure 9 shows our results for the expansion around SNM(blue) and the expansion around PNM (red). Both ex-pansions agree, and the differences are smaller than theuncertainties by a factor of 2 −
3. We also show re-sults for the six Hamiltonians. Their spread is muchsmaller than the uncertainties of the data or the model.This is because the latter are computed as arithmeticaverages of the error bars of the global symmetry en-ergy and the quadratic contribution. At n empsat , we ob-tain from our model e sym ,nq = 1 . e pot sym ,nq =0 . e pot ∗ sym ,nq = − . e sym ,nq = 0 . e pot sym ,nq = 0 . e pot ∗ sym ,nq = − . −
5% to the symmetry energy. They originate mainlyfrom the kinetic energy, since e pot sym ,nq remains close to zero across all densities. The associated uncertaintiesrepresent the model dependence explored in the presentapproach. Our estimates for the non-quadratic contribu-tions to the symmetry energy and the NEP are summa-rized in Table V. We also compare the present NEPs tothe three selected Skyrme interactions, showing a goodagreement between the microscopic results and the EDFapproaches.We calculate the quartic contribution to the symme-try energy e sym , using Eq. (51), and show the result-ing NEPs in Table V. We find that the quartic term tothe symmetry energy accounts for about 60 −
70% ofthe total non-quadratic contribution, while the remain-ing 30 −
40% originate from higher-order contributions.The convergence of these additional contributions is dis-cussed in Ref. [36]. We stress that this does not includeany logarithmic contribution because such a contributionwould vanish in PNM.A recent analysis based on a general EDF approach—which was optimized to the properties of finite nuclei—concluded that quartic terms ∝ δ have little impacton nuclei [92]. The result was interpreted as a conse-quence of the fact that the asymmetry δ in finite nu-clei is small: for Z > − .
33 and+0 .
38 in the latest Atomic Mass Data Center mass tableAME2016 [93]. A quartic term was, however, found tobe important to correctly reproduce the PNM energy perparticle. In order to reproduce the PNM energy per par-ticle predicted in Ref. [30], Ref. [92] found a quartic termof e sym , = 2 .
635 MeVat n = 0 . − . This term was theonly non-quadratic contribution considered in Ref. [92],and is consistent within our upper 68% confidence inter-val for the non-quadratic contribution to the symmetryenergy. The higher value for e sym , obtained in Ref. [92]might be related to the larger value in the PNM energyper particle obtained in Ref. [30], as shown in Fig. 1.This affects the symmetry energy because the contribu-tion e sym , is needed to correctly reproduce the PNMEOS. Both, the PNM energy per particle in Ref. [30] and e sym , obtained in Ref. [92] are ∼ D. Logarithmic contribution to the symmetryenergy e sym,log The leading-order logarithmic contribution to the sym-metry energy, see Eq. (17), was suggested to be of theform δ log | δ | [35, 36]. It, therefore, vanishes in bothSNM and PNM, and data at finite isospin asymmetryare required to determine its magnitude. Such a loga-rithmic term would appear by a characteristic arch-likestructure in the δ -dependent residuals between the dataand a model without the logarithmic term. Figure 10shows these residuals as a function of δ at three different6 n (fm )10123 e s y m , n q ( M e V ) H1H2H3H4H5H7 n (fm )10123 e p o t s y m , n q ( M e V ) Data (delta) (68% CL)Data (eta) (68% CL) n (fm )10123 e p o t * s y m , n q ( M e V ) model (delta) (68% CL)model (eta) (68% CL) FIG. 9. Non-quadratic terms e sym ,nq calculated via an expansion around SNM (blue) and PNM (red). For the right panel, theLandau mass is described by a linear fit. The coloured lines depict results for the six individual Hamiltonians.TABLE V. Posteriors of non-quadratic and quartic NEPs obtained from the analysis of e sym ,nq using the δ and η expansions,and e sym , obtained from the η expansion only. For the extraction of e sym ,nq via the η expansion, the values inside the squarebrackets are obtained from a fit to the data in order to provide a direct comparison to the corresponding analysis of e sym , .Values in parenthesis represent the ± σ uncertainties. The NEPs for the following three Skyrme-type interactions are given:NRAPR [3], LNS5 [84] and SAMI [85].Non-quadratic E sym ,nq L sym ,nq K sym ,nq Q sym ,nq Z sym ,nq contribution (MeV) (MeV) (MeV) (MeV) (MeV)SNM 1 . − . − (cid:2) . (cid:3) (cid:2) . (cid:3) (cid:2) − (cid:3) (cid:2) (cid:3) (cid:2) (cid:3) NRAPR [3] 1 .
40 5 6 − − .
70 6 9 − .
08 3 2 2 − E sym , L sym , K sym , Q sym , Z sym , contribution (MeV) (MeV) (MeV) (MeV) (MeV)PNM 1 . . − .
95 3 4 − − .
17 5 6 − .
70 2 2 1 − densities. For asymmetric matter, we use y model ( n, δ ) = y SNM ( n ) + y sym , ( n ) δ + y sym,nq ( n ) δ , (52)where e SNM , e sym , and e sym,nq are given by Eqs. (30),(42), and (19). Note that in the model (52), the fourth- order δ term also includes possible higher-order contri-butions (like, for instance, a δ term) contained in theterm y sym,nq . The different panels in Fig. 10 show theresults at three densities, n = 0 .
06, 0 .
12, and 0 .
16 fm − (from the top to the bottom panel), and for the threechoices for the variable y : e , e pot , and e pot* (from left7 R = M o d e l - D a t a ( M e V ) n = 0.06 fm y = e n = 0.06 fm y = e pot n = 0.06 fm y = e pot* R = M o d e l - D a t a ( M e V ) n = 0.12 fm y = e n = 0.12 fm y = e pot Mean R (Delta)Mean R (Eta) n = 0.12 fm y = e pot* ± (R) (Delta)± (R) (Eta) R = M o d e l - D a t a ( M e V ) n = 0.16 fm y = e H1H2H3 n = 0.16 fm y = e pot H4H5H7 n = 0.16 fm y = e pot* FIG. 10. Residuals of the model, see Eq. (52), with respect to the data as a function of the asymmetry parameter δ for differentvalues of the density: n = 0 .
06 fm − (first row), n = 0 .
12 fm − (second row), and n = 0 .
16 fm − (third row). The resultsare shown for the two different calculations of y sym,2 and y sym,nq —the expansions around PNM (red) and SNM (blue). Thedifferent columns correspond to e , e pot , and e pot , ∗ . The coloured lines depict the residuals of the fit for each Hamiltonian. Inthe last column, the black-dashed lines represent the upper and lower limits of the uncertainty in the residuals, respectively,by disregarding the uncertainties in the effective masses. to right) as a function of the isospin asymmetry δ . Thesquares (shaded bands) represent the mean (68% confi-dence level) of the residuals. The presence of logarithmicterms would appear as a systematic deviation of the theseresiduals from zero in asymmetric matter. However, thisis not what we observe at the three considered densi-ties and for all energy observables. Instead, we find theresiduals to be compatible with zero and almost flat as afunction of the isospin asymmetry. This is also the casefor the results for each Hamiltonian, which we show ascoloured lines. The results for the individual Hamilto-nians vanish on average, but the uncertainty bands re-main quite sizable, about ± − ∼ . Contribu- This is the n3lo450 interaction constructed in Refs. [37, 57, 80],which is not considered in this work. tions of the order of ∼ . VII. IMPACT ON THE NEUTRON-STARCRUST-CORE TRANSITION
In this section, we study the impact of the non-quadratic contribution to the symmetry energy on thecrust-core transition in neutron stars, for which the sym-metry energy plays an important role [94–96]. This tran-sition occurs at the core-crust transition density n cc withan isospin asymmetry δ cc that is determined by the beta-equilibrium. The parameters n cc and δ cc can be obtainedfrom uniform matter by determining the density at whichmatter becomes unstable with respect to density fluctu-ations (spinodal instability) [94].In multi-component matter, e.g., matter that consistsof neutrons and protons, this spinodal density is deter-mined from the curvature (Hessian) matrix C , defined asthe second derivative of the energy density with respectto the component densities [94]. From the eigenvaluesof C , one can determine the stability of matter: if alleigenvalues are positive, i.e. if C is positive semi-definite,the matter exhibits a local stability against density fluc-8 n (fm ) only + FIG. 11. Predictions for beta equilibrium in low-density uni-form matter obtained by solving Eq. (54), and for the spin-odal density by solving Eq. (53). The intersection denotes thecrust-core transition, as indicated by a dot in the inset. Thequadratic approximation (red band) is compared to the casewhere quartic contributions are included (blue band). tuations of all components in any combination, while achange of sign for any eigenvalue triggers an instabilitywith respect to density fluctuations indicated by its asso-ciated eigenvector. The change of sign of the eigenvaluescan be extracted from the determinant of C , which readsin nuclear matter,det C ( n, δ ) = ∂µ n ∂n n ∂µ p ∂n p − ∂µ n ∂n p ∂µ p ∂n n , (53)where µ n and µ p are the neutron and proton chemicalpotentials. For simplicity, we have neglected the finite-size contribution from the Coulomb interaction as wellas the gradient density terms induced by the finite rangeof nuclear interactions. It is expected that these termsreduce the spinodal density by only about 0 .
01 fm − [94,96].In sub-saturation asymmetric matter, the equilibriumstate is the state that satisfies the chemical potentialequilibrium µ n = µ p + µ e , at fixed baryon number n = n n + n p and charge neutrality n e = n p . At zerotemperature, and considering relativistic electrons, thissystem of equations reduces to a single non-linear equa-tion, s m e + (cid:18) π − δ β ) n (cid:19) / = 2 ∂e ( n, δ β ) ∂δ , (54)whose solution, δ β ( n ), is obtained by using a combinationof the bisection and the secant methods implemented inthe Python package of Ref. [97]. Then, we define thecrust-core transition as the solution ( n cc , δ cc ) to both,the instability onset criterion, det C = 0, and the betaequilibrium condition, e.g., Eq. (54). Equivalently, n cc is defined as the spinodal density in beta equilibrium, TABLE VI. Crust-core transition density and isospin asym-metry, n cc and δ cc , respectively, for the purely quadratic case( δ only) and for the case including the quartic contribution( δ + δ ), see Eq. (52).Model n cc (fm − ) δ cc δ only 0 . . δ + δ . . where δ cc = δ β ( n cc ). Figure 11 shows the intersectionbetween these two determinations. We investigate thepurely quadratic approximation for the symmetry energy,with the NEPs given in Table IV, and with the quarticterms from Table V included. For all cases, the refer-ence MM in SNM is determined by the best fit given inTable III for the scaling 3 ∗ .When we include quartic contributions, the spinodaldensity in neutron-rich matter is increased compared tothe case where only the quadratic term is considered.This is because the quartic term increases the symme-try energy. For the same reason, the isospin asymme-try is decreased when non-quadraticities are included.Our results are summarized in Table VI, and depictedin Fig. 11 by the blue and red points. They are in agree-ment with the predictions of, e.g., Refs. [95, 96] with L sym , ≈
45 MeV. From the comparison of our resultswith and without the quartic term, we find that the tran-sition density changes by ∼
5% while δ cc changes by only ∼ VIII. SUMMARY AND CONCLUSIONS
We have analyzed the properties of asymmetric nuclearmatter based on MBPT calculations [20] for six com-monly used chiral EFT Hamiltonians with NN and 3Ninteractions. The global symmetry energy, i.e., the dif-ference between EOS in the limits of PNM and SNM, aswell as its quadratic and quartic contributions have beendetermined with theoretical uncertainty estimates. Wehave calculated the quadratic contribution to the sym-metry energy from the usual expansion around SNM, andhave also employed a non-standard approach using an ex-pansion for small proton fractions around PNM. The twoapproaches are in excellent agreement. Furthermore, wehave investigated the strength of the non-quadraticitiesas well as their model dependence. The quartic contri-bution to the symmetry energy was found to be about1 . . δ , indicating no systematic deviation fromzero as expected for a logarithmic contribution. However,we also saw that present uncertainties, indicated by thedispersion of the six Hamiltonians of about 1 − ∼ δ decreased by ∼
1% when non-quadraticities are included.Hence, these contributions are only small corrections, butneed to be included for a precise calculation of the core-crust transition properties.To gauge the full theoretical uncertainties of the non-quadratic contributions to the symmetry energy, futureanalyses need to explore a wider range of nuclear inter-actions and additional asymmetric-matter calculationsusing different many-body approaches and regulariza-tion schemes. In particular, this requires the devel-opment of improved chiral NN and 3N interactions upto N LO [56, 98, 99], which will enable order-by-orderanalyses of the neutron-rich matter EOS with statis-tically meaningful uncertainty estimates derived fromchiral EFT [22, 48]. Our work provides a framework,e.g., Python codes [28] and Supplemental Material re-lated to our data, for future investigations of the isospin-dependence of nuclear matter.
ACKNOWLEDGMENTS
We thank K. Hebeler, J. Lattimer, and A. Schwenkfor fruitful discussions. R.S. is supported by the PHASTdoctoral school (ED52) of
Universit´e de Lyon . R.S. andJ.M. are both associated to the CNRS/IN2P3 NewMACproject, and are also grateful to PHAROS COST Ac-tion MP16214 and to the LABEX Lyon Institute of Ori-gins (ANR-10-LABX-0066) of the
Universit´e de Lyon forits financial support within the program
Investissementsd’Avenir (ANR-11-IDEX-0007) of the French govern-ment operated by the National Research Agency (ANR).C.D. acknowledges support by the Alexander von Hum-boldt Foundation through a Feodor-Lynen Fellowshipand the US Department of Energy, the Office of Science,the Office of Nuclear Physics, and SciDAC under awardsDE-SC00046548 and DE-AC02-05CH11231. The workof I.T. was supported by the U.S. Department of En-ergy, Office of Science, Office of Nuclear Physics, undercontract No. DE-AC52-06NA25396, by the NUCLEI Sci-DAC program, and by the LDRD program at LANL.
Appendix A: Details of the Fitting Procedure
To obtain the parameter values, we perform in this pa-per a Bayesian analysis using the nonlinear least-square fitting package lsqfit developed in Python [100]. lsq-fit uses scipy ’s least-squares minimization routine tooptimize our model with respect to data y i (here, thepredictions from each of the six Hamiltonians consideredin this work). The index i runs over all data samplingpoints (e.g., the density grid). The propagation of theGaussian uncertainties from the model parameters, p α ,to the functions of those parameters, f ( { p α } ), requirescalculating derivatives of those functions with respect tothe parameters, which is achieved by automatic differen-tiation. In this work, the spread of the six Hamiltoniansis interpreted as a fair representation of the width σ ofa normal distribution that reflects theoretical uncertain-ties.The fitting procedure requires the minimization ofthe objective function, hereafter noted the χ function,which, in general, receives contributions from both the in-put data ( χ ) and the prior information on the modelparameters ( χ ). We can write these two contribu-tions as χ = X ij ∆ y ( p ) i cov − ij ∆ y ( p ) j , (A1) χ = X α ( p α − µ α ) σ α , (A2)with the total χ being the sum of the two terms, χ = χ + χ . (A3)In Eq. (A1), ∆ y ( p ) i = f ( { p α } ) i − E ( y i ), where the ex-pectation value E ( y i ) is defined as E ( y i ) = 16 X µ =1 ( y i ) µ , (A4)the summation index µ runs over the six nuclear Hamil-tonians, and cov ij is the co-variance matrix between thedata points y i and y j , which is given bycov ij ≡ cov( y i , y j ) = E ( y i y j ) − E ( y i ) E ( y j ) . (A5)The correlation matrix corr ij is then defined by its matrixelements, corr ij = cov ij √ cov ii cov jj . (A6)If the data are independent from each other, the co-variance and correlation matrices are diagonal and χ is associated with a normal distribution, as it is for χ . Here, the data are not independent because ofcorrelations in density: the knowledge of the predic-tions of the Hamiltonians at a few density points canbe used to determine other points by interpolation or—to some extent—by extrapolation. We, thus, expect anon-diagonal co-variance matrix. While the Hamiltoni-ans are by construction strongly correlated, we do notaccount for their correlations in the co-variance matrix.0It only includes the correlations between the differentdata points, see Eq. (A5).In Eq. (A2), µ α is the prior mean value of the param-eter p α , and σ α is its standard deviation. Note that, inthis work, we only consider uncorrelated Gaussian priordistributions for our fit parameters.The best-fit values of the parameters ¯ p α are those thatminimize the total χ . The inverse co-variance matrixcorresponding to the posterior distribution of the best-fitparameters is defined as,(cov − p ) αβ = ∂χ ∂p α ∂p β (¯ p ) . (A7)This expression is used by lsqfit to define the uncertain-ties in the fit parameters p α . Appendix B: Correlations in the data sample
The data analysis performed in this work involves para-metric fits to data that are highly correlated across densi-ties (or equivalently, Fermi momenta). It was found thatthis correlation had to be taken into account in order toachieve compatibility between the obtained posterior dis-tributions and the data, as shown in Figs. 5, 6 and 8. Inthis appendix, we provide an estimate of this correlation.As mentioned in Appendix A, the correlation in the in-put data is captured by the covaraince matrix (A5). Thestrong positive correlation in the data points across den-sities results in large off-diagonal elements with respectto the diagonal ones. We analyse the strength of theseoff-diagonal elements as follow. We first diagonalize thecorrelation matrix (A6) and obtain the eigenvalues andassociated eigenvectors. In the case of maximum corre-lation all but one eigenvalue are 0. In our case, we foundthat the largest eigenvalue contributes to about 96% ofthe trace of the correlation matrix. Then, the eigenvectorcorresponding to this dominant eigenvalue captures thecorrelation (or mixing) between the data, which needs tobe quantified. Let λ be this eigenvector with components( λ , . . . , λ N D ), where N D is the number of data points.In the case of maximum correlation, the quantities | λ | , . . . , | λ N D | have zero dispersion about their average value,and in the case of minimum correlation, this dispersiontakes a maximum value. We thus define the dispersionabout the mean as d = 1 N D N D X i =1 ( | λ i | − ¯ λ ) , (B1)where ¯ λ is the average value of the | λ i | . In the case ofa maximum correlation d = 0, whereas in the case ofminimum correlation d = N D [1 − N D ].We then define the level of correlation as l corr ≡ − N D d , (B2)where l corr = 1 in the case of maximum correlation and,for zero correlation l corr = 1 /N D , which approaches 0 in the limit of infinitely many data points. This parameteris similar to the correlation coefficient, ρ which one in-troduces in the case of identical off-diagonal elements inthe covariance matrix.The values of l corr for the various fits presented in thepaper are as follows. In Sec. III C, fits to the effectivemasses were presented in Fig. 5. For the linear fit, l corr is0 .
96 in SNM and 0 .
92 in PNM; for the quadratic fit, 0 . .
54 in PNM. For the fits shown in Fig. 6, l corr is 0 . , . , .
46 in the scalings 1, 2 and 3, respec-tively, in SNM. Similarly in PNM, it is 0 . , . , . l corr = 0 .
50 for theanalysis with the expansion around SNM and l corr = 0 . Appendix C: Goodness of the fit and Q–Q plots
In this section we quantify the goodness of the fitsusing the Q–Q plot method. We consider the two il-lustrative fits in Fig. 12: the fit of the Landau mass inSNM using a quadratic density functional (top panel inFig. 5) and the energy per particle in PNM as a functionof density (lower right panel in Fig. 6).For the analysis of the Landau mass, we found χ / dof ≈ .
1, and for the energy per particle, we findthat χ / dof ≈ .
4. These numbers are already indica-tions for a successful fit.We further investigate the goodness of the fit by exam-ining the normalized residuals of the fit, which are definedas the differences between the data and the best fit model.These differences are expressed in the basis formed by theeigenvectors of the correlation matrix (A6). Then, eachcomponent of the eigenvectors are divided by the squareroot of the corresponding eigenvalue. The assumption ofthis analysis for a good fit is that the normalized residualsshould be uncorrelated and randomly distributed aroundthe mean value following a normal distribution.Q–Q plots are a way for testing that the residuals fol-low a normal distribution. They are obtained as follow.First, the residuals are ordered from the smallest to thelargest on the y -axis, then they are plotted against anordered list of samples drawn from a normal distribution(centered at zero, with a width of one). If the residualsare perfectly normal distributed, then the result alignson a straight line with slope 1. This alignment can becaptured by the coefficient of determination, R definedas R = 1 − P ni =1 ( r i − ˆ r i ) P ni =1 ( r i − ¯ r ) (C1)where n is the number of residuals, r i are the residuals,ˆ r i are the values expected for those residuals and ¯ r =1 Theoretical quantiles O r d e r e d f i t r e s i d u a l s R =0.89
Theoretical quantiles O r d e r e d f i t r e s i d u a l s R =0.83
FIG. 12. Q–Q plot for the analysis of the Landau mass inSNM (top row) and energy per particle in PNM (bottom row).The reported R score, see Eq. (C1), is a measure of how wellthe blue points lie on the dashed red line. n P ni =1 r i . The best possible R score is 1 and it can benegative for arbitrarily worse fits.The top panel of Fig. 12 shows the Q–Q plot for ouranalysis of the Landau mass in SNM. The ideal straightline with unit slope is shown as the dashed-red line, whilethe ordered residuals are shown as blue circles. The bluesolid line is a fit to the blue points. The R score, seeEq. (C1), is 0 .
89, which leads us to the conclusion thatthe ordered fit residuals are consistent with a normal dis-tribution with a mean one, so close to χ / dof = 1. Wealso give similar results for the analysis of the energy perparticle in PNM in the bottom panel of Fig. 12. Here,the coefficient of determination is 0 .
83, indicating a goodfit. [1] M. Oertel, M. Hempel, T. Kl¨ahn, and S. Typel, Rev.Mod. Phys. , 015007 (2017).[2] G. Fiorella Burgio and A. F. Fantina, “Nuclear equa-tion of state for compact stars and supernovae,” in ThePhysics and Astrophysics of Neutron Stars , edited byL. Rezzolla, P. Pizzochero, D. I. Jones, N. Rea, andI. Vida˜na (Springer International Publishing, Cham,2018) pp. 255–335.[3] A. W. Steiner, M. Prakash, J. M. Lattimer, and P. J. El-lis, Phys. Rept. , 325 (2005), arXiv:nucl-th/0410066.[4] M. Baldo and G. Burgio, Prog. Part. Nucl. Phys. ,203 (2016), arXiv:1606.08838 [nucl-th].[5] K. Hebeler, J. Lattimer, C. Pethick, and A. Schwenk,Phys. Rev. Lett. , 161102 (2010), arXiv:1007.1746[nucl-th].[6] S. Gandolfi, J. Carlson, and S. Reddy, Phys. Rev. C , 032801 (2012), arXiv:1101.1921 [nucl-th].[7] A. Steiner and S. Gandolfi, Phys. Rev. Lett. , 081102(2012), arXiv:1110.4142 [nucl-th]. [8] C. Horowitz, E. Brown, Y. Kim, W. Lynch, R. Michaels,A. Ono, J. Piekarewicz, M. Tsang, and H. Wolter, J.Phys. G , 093001 (2014), arXiv:1401.5839 [nucl-th].[9] N. S. A. C. (NSAC), , (2015).[10] S. Abrahamyan et al. , Phys. Rev. Lett. , 112502(2012), arXiv:1201.2568 [nucl-ex].[11] C. J. Horowitz, K. S. Kumar, and R. Michaels, Eur.Phys. J. A , 48 (2014), 1307.3572.[12] X. Roca-Maza, M. Brenna, B. Agrawal, P. Bortignon,G. Col, L.-G. Cao, N. Paar, and D. Vretenar, Phys.Rev. C , 034301 (2013), arXiv:1212.4377 [nucl-th].[13] M. Tsang et al. , Phys. Rev. C , 015803 (2012),arXiv:1204.0466 [nucl-ex].[14] J. M. Lattimer and A. W. Steiner, Eur. Phys. J. A ,40 (2014), arXiv:1403.1186 [nucl-th].[15] J. M. Lattimer and Y. Lim, Astrophys. J. , 51(2013). [16] J. Margueron, R. Hoffmann Casali, and F. Gulminelli,Phys. Rev. C , 025805 (2018), arXiv:1708.06894 [nucl-th].[17] M. Dutra, O. Lourenco, J. Sa Martins, A. Delfino,J. Stone, and P. Stevenson, Phys. Rev. C , 035201(2012), arXiv:1202.3902 [nucl-th].[18] M. Dutra, O. Loureno, S. Avancini, B. Carlson,A. Delfino, D. Menezes, C. Providncia, S. Typel,and J. Stone, Phys. Rev. C , 055203 (2014),arXiv:1405.3633 [nucl-th].[19] I. Tews, T. Krger, K. Hebeler, and A. Schwenk, Phys.Rev. Lett. , 032504 (2013), arXiv:1206.0025 [nucl-th].[20] C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev.C , 054314 (2016), arXiv:1510.06728 [nucl-th].[21] D. Lonardoni, I. Tews, S. Gandolfi, and J. Carlson,Phys. Rev. Res. , 022033 (2020), arXiv:1912.09411[nucl-th].[22] C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R.Phillips, (2020), arXiv:2004.07232 [nucl-th].[23] A. W. Steiner, Phys. Rev. C , 045808 (2006),arXiv:nucl-th/0607040.[24] L.-W. Chen, B.-J. Cai, C. M. Ko, B.-A. Li, C. Shen, andJ. Xu, Phys. Rev. C , 014322 (2009), arXiv:0905.4323[nucl-th].[25] B.-J. Cai and L.-W. Chen, Phys. Rev. C , 024302(2012), arXiv:1111.4124 [nucl-th].[26] B.-J. Cai and B.-A. Li, Phys. Rev. C , 011601 (2015),arXiv:1503.01167 [nucl-th].[27] W. Seif and D. Basu, Phys. Rev. C , 028801 (2014),arXiv:1401.0090 [nucl-th].[28] R. Somasundaram, “Somasundaram-rahul/nuclear-symmetry-energy,” (2020).[29] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev.Mod. Phys. , 121 (2003).[30] G. Wlaz lowski, J. Holt, S. Moroz, A. Bulgac, andK. Roche, Phys. Rev. Lett. , 182503 (2014),arXiv:1403.3753 [nucl-th].[31] I. Bombaci and U. Lombardo, Phys. Rev. C , 1892(1991).[32] C.-H. Lee, T. Kuo, G. Li, and G. Brown, Phys. Rev. C , 3488 (1998).[33] T. Frick, H. Muther, A. Rios, A. Polls, and A. Ramos,Phys. Rev. C , 014313 (2005), arXiv:nucl-th/0409067.[34] I. Vidana, C. Providencia, A. Polls, and A. Rios, Phys.Rev. C , 045806 (2009), arXiv:0907.1165 [nucl-th].[35] N. Kaiser, Phys. Rev. C , 065201 (2015),arXiv:1504.00604 [nucl-th].[36] C. Wellenhofer, J. W. Holt, and N. Kaiser, Phys. Rev.C , 055802 (2016), arXiv:1603.02935 [nucl-th].[37] D. R. Entem and R. Machleidt, Phys. Rev. C , 041001(2003).[38] K. Hebeler, S. K. Bogner, R. J. Furnstahl, A. Nogga,and A. Schwenk, Phys. Rev. C , 031301 (2011).[39] M. C. M. Rentmeester, R. G. E. Timmermans, and J. J.de Swart, Phys. Rev. C , 044001 (2003).[40] H.-W. Hammer, A. Nogga, and A. Schwenk, Rev. Mod.Phys. , 197 (2013).[41] K. Hebeler, J. Holt, J. Menndez, and A. Schwenk,Annual Review of Nuclear and Particle Science , 457 (2015), https://doi.org/10.1146/annurev-nucl-102313-025446.[42] I. Tews, Z. Davoudi, A. Ekstrm, J. D. Holt, and J. E.Lynn, J. Phys. G , 103001 (2020), arXiv:2001.03334 [nucl-th].[43] E. Epelbaum, H.-W. Hammer, and U.-G. Meißner, Rev.Mod. Phys. , 1773 (2009), arXiv:0811.1338 [nucl-th].[44] R. Machleidt and D. R. Entem, Phys. Rept. , 1(2011), arXiv:1105.2919 [nucl-th].[45] S. Bogner, R. Furnstahl, and A. Schwenk, Progress inParticle and Nuclear Physics , 94 (2010).[46] R. J. Furnstahl and K. Hebeler, Reports on Progress inPhysics , 126301 (2013).[47] E. Epelbaum, H. Krebs, and U. Meiner, Eur. Phys. J.A , 53 (2015), arXiv:1412.0142 [nucl-th].[48] C. Drischler, J. Melendez, R. Furnstahl, andD. Phillips, (2020), arXiv:2004.07805 [nucl-th].[49] A. Carbone, A. Rios, and A. Polls, Phys. Rev. C ,044302 (2013), arXiv:1307.1889 [nucl-th].[50] A. Carbone, A. Cipollone, C. Barbieri, A. Rios,and A. Polls, Phys. Rev. C , 054326 (2013),arXiv:1310.3688 [nucl-th].[51] G. Hagen, T. Papenbrock, A. Ekstr¨om, K. Wendt,G. Baardsen, S. Gandolfi, M. Hjorth-Jensen, andC. Horowitz, Phys. Rev. C , 014319 (2014),arXiv:1311.2925 [nucl-th].[52] A. Ekstr¨om, G. R. Jansen, K. A. Wendt, G. Hagen,T. Papenbrock, B. D. Carlsson, C. Forss´en, M. Hjorth-Jensen, P. Navr´atil, and W. Nazarewicz, Phys. Rev. C , 051301 (2015).[53] K. Hebeler and A. Schwenk, Phys. Rev. C , 014314(2010).[54] C. Wellenhofer, J. W. Holt, N. Kaiser, and W. Weise,Phys. Rev. C , 064009 (2014), arXiv:1404.2136 [nucl-th].[55] J. Holt and N. Kaiser, Phys. Rev. C , 034326 (2017),arXiv:1612.04309 [nucl-th].[56] C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev.Lett. , 042501 (2019), arXiv:1710.08220 [nucl-th].[57] L. Coraggio, J. W. Holt, N. Itaco, R. Machleidt, L. E.Marcucci, and F. Sammarruca, Phys. Rev. C ,044321 (2014).[58] J. W. Holt, N. Kaiser, and W. Weise, Prog. Part. Nucl.Phys. , 35 (2013), arXiv:1304.6350 [nucl-th].[59] A. Gezerlis, I. Tews, E. Epelbaum, S. Gandolfi,K. Hebeler, A. Nogga, and A. Schwenk, Phys. Rev.Lett. , 032501 (2013), arXiv:1303.6243 [nucl-th].[60] A. Gezerlis, I. Tews, E. Epelbaum, M. Freunek, S. Gan-dolfi, K. Hebeler, A. Nogga, and A. Schwenk, Phys.Rev. C , 054323 (2014), arXiv:1406.0454 [nucl-th].[61] I. Tews, S. Gandolfi, A. Gezerlis, and A. Schwenk, Phys.Rev. C , 024305 (2016), arXiv:1507.05561 [nucl-th].[62] J. Lynn, I. Tews, J. Carlson, S. Gandolfi, A. Gezerlis,K. Schmidt, and A. Schwenk, Phys. Rev. Lett. ,062501 (2016), arXiv:1509.03470 [nucl-th].[63] C. Drischler, V. Som`a, and A. Schwenk, Phys. Rev. C , 025806 (2014).[64] C. Wellenhofer, J. W. Holt, and N. Kaiser, Phys. Rev.C , 015801 (2015).[65] J. Simonis, K. Hebeler, J. Holt, J. Menendez,and A. Schwenk, Phys. Rev. C , 011302 (2016),arXiv:1508.05040 [nucl-th].[66] G. Hagen, A. Ekstr¨om, C. Forss´en, G. R. Jansen,W. Nazarewicz, T. Papenbrock, K. A. Wendt,S. Bacca, N. Barnea, B. Carlsson, C. Drischler,K. Hebeler, M. Hjorth-Jensen, M. Miorelli, G. Orlan-dini, A. Schwenk, and J. Simonis, Nat. Phys. , 186(2016), 1509.07169. [67] R. F. Garcia Ruiz, M. L. Bissell, K. Blaum,A. Ekstr¨om, N. Fr¨ommgen, G. Hagen, M. Ham-men, K. Hebeler, J. D. Holt, G. R. Jansen,M. Kowalska, K. Kreim, W. Nazarewicz, R. Neu-gart, G. Neyens, W. N¨ortersh¨auser, T. Papenbrock,J. Papuga, A. Schwenk, J. Simonis, K. A. Wendt, andD. T. Yordanov, Nat. Phys. , 594 (2016), 1602.07906.[68] G. Hagen, G. R. Jansen, and T. Papenbrock, Phys.Rev. Lett. , 172501 (2016), 1605.01477.[69] J. Simonis, S. R. Stroberg, K. Hebeler, J. D. Holt,and A. Schwenk, Phys. Rev. C , 014303 (2017),1704.02915.[70] J. Birkhan, M. Miorelli, S. Bacca, S. Bassauer, C. A.Bertulani, G. Hagen, H. Matsubara, P. von Neumann-Cosel, T. Papenbrock, N. Pietralla, V. Y. Ponomarev,A. Richter, A. Schwenk, and A. Tamii, Phys. Rev. Lett. , 252501 (2017), 1611.07072.[71] T. D. Morris, J. Simonis, S. R. Stroberg, C. Stumpf,G. Hagen, J. D. Holt, G. R. Jansen, T. Papenbrock,R. Roth, and A. Schwenk, Phys. Rev. Lett. , 152503(2018).[72] J. Holt, S. Stroberg, A. Schwenk, and J. Simonis,(2019), arXiv:1905.10475 [nucl-th].[73] M. Mougeot et al. , Phys. Rev. C , 014301 (2020),arXiv:2006.02712 [nucl-ex].[74] S. Kaufmann et al. , Phys. Rev. Lett. , 132502(2020), arXiv:2003.06353 [nucl-ex].[75] J. W. Holt, N. Kaiser, and W. Weise, Phys. Rev. C ,024002 (2010), arXiv:0910.1249.[76] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall,Phys. Rev. C , 1804 (1998).[77] R. B. Wiringa, V. Stoks, and R. Schiavilla, Phys. Rev.C , 38 (1995), arXiv:nucl-th/9408016.[78] B. Pudliner, V. Pandharipande, J. Carlson, and R. B.Wiringa, Phys. Rev. Lett. , 4396 (1995), arXiv:nucl-th/9502031.[79] L. Coraggio, A. Covello, A. Gargano, N. Itaco, D. En-tem, T. Kuo, and R. Machleidt, Phys. Rev. C ,024311 (2007), arXiv:nucl-th/0701065.[80] L. Coraggio, J. Holt, N. Itaco, R. Machleidt, andF. Sammarruca, Phys. Rev. C , 014322 (2013),arXiv:1209.5537 [nucl-th].[81] J. Jeukenne, A. Lejeune, and C. Mahaux, Phys. Rept. , 83 (1976). [82] K. Hassaneen and H. Muther, Phys. Rev. C , 054308(2004), arXiv:nucl-th/0408035.[83] E. van Dalen, C. Fuchs, and A. Faessler, Phys. Rev.Lett. , 022302 (2005), arXiv:nucl-th/0502064.[84] L. Cao, U. Lombardo, C. Shen, and N. Van Giai, Phys.Rev. C , 014313 (2006), arXiv:nucl-th/0507071.[85] X. Roca-Maza, G. Colo, and H. Sagawa, Phys. Rev. C , 031306 (2012), arXiv:1205.3958 [nucl-th].[86] W. Zuo, I. Bombaci, and U. Lombardo, Phys. Rev. C , 024605 (1999), arXiv:nucl-th/0102035.[87] W. Zuo, A. Lejeune, U. Lombardo, and J. Mathiot,Nucl. Phys. A , 418 (2002), arXiv:nucl-th/0202076.[88] B. Fornberg, Math. Comp. , 699 (1988).[89] J. Margueron and F. Gulminelli, Phys. Rev. C ,025806 (2019), arXiv:1807.01729 [nucl-th].[90] J. Margueron, R. Hoffmann Casali, and F. Gulminelli,Phys. Rev. C , 025806 (2018), arXiv:1708.06895 [nucl-th].[91] B.-A. Li, P. G. Krastev, D.-H. Wen, and N.-B. Zhang,Eur. Phys. J. A , 117 (2019), arXiv:1905.13175 [nucl-th].[92] A. Bulgac, M. M. Forbes, S. Jin, R. Navarro Perez,and N. Schunck, Phys. Rev. C , 044313 (2018),arXiv:1708.08771 [nucl-th].[93] G. Audi, F. G. Kondev, M. Wang, W. Huang, andS. Naimi, Chinese Physics C , 030001 (2017).[94] C. Pethick, D. Ravenhall, and C. Lorenz, NuclearPhysics A , 675 (1995).[95] C. Ducoin, J. Margueron, and C. Providencia, EPL ,32001 (2010), arXiv:1004.5197 [nucl-th].[96] C. Ducoin, J. Margueron, C. Providencia, and I. Vi-dana, Phys. Rev. C , 045810 (2011), arXiv:1102.1283[nucl-th].[97] P. Lepage, C. Gohlke, and D. Hackett, “gplepage/gvar:gvar version 11.7,” (2020).[98] J. Hoppe, C. Drischler, K. Hebeler, A. Schwenk,and J. Simonis, Phys. Rev. C , 024318 (2019),arXiv:1904.12611 [nucl-th].[99] T. Hther, K. Vobig, K. Hebeler, R. Machleidt,and R. Roth, Phys. Lett. B , 135651 (2020),arXiv:1911.04955 [nucl-th].[100] P. Lepage and C. Gohlke, “gplepage/lsqfit: lsqfit version11.5.1,” (2020). onstraints on the nuclear symmetry energy from asymmetric-matter calculationswith chiral NN and 3N interactions (Supplementary Material) R. Somasundaram, ∗ C. Drischler,
2, 3, † I. Tews, ‡ and J. Margueron § Univ Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3,IP2I Lyon, UMR 5822, F-69622, Villeurbanne, France Department of Physics, University of California, Berkeley, CA 94720, USA Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA (Dated: September 11, 2020)We briefly describe the GitHub repository where the data and Python codes used in this paperare provided under the MIT license.
I. DESCRIPTION OF THE GITHUBREPOSITORY
The repository [1] contains the Python codes used toperform the analysis as well as generate all the figurespresented in the publication. This repository is publiclyavailable on GitHub [1] and distributed under the MITlicense.The model parameters are obtained from a Bayesiananalysis using the nonlinear least-square fitting lsqfit package developed in Python [2], which usesthe scipy least-squares minimization routine to optimizeour model with respect to data.
II. LIST OF FOLDERS
There are two main directories at the root of the repos-itory: data and results. • data: contains the predictions (that we treat asdata) from each of the six Hamiltonians consideredin this work. They are given as ASCII files. – data/Effective mass: consists of the singleparticle energies for the six Hamiltonians insymmetric and neutron matter. – data/EOS Drischler: contains the data filesfor the energy per particle for different valuesof iso-spin asymmetry. • results: contains the figures as pdf files. The follow-ing sub-folders contain the various figures presentedin the paper as follows, – results/3 scales: Figs. 1 and 6. – results/crust core: Fig. 11. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] – results/effective mass: Figs. 2, 3, 4 and 5. – results/esym esym2: Figs. 7 and 8. – results/non quadraticities: Figs. 9 and 10. – results/goodness: plots presented in the ap-pendix of the paper. III. LIST OF PYTHON CODES
The Python codes are located in the main folder. Theydepend on the following Python libraries: scipy, numpy,matplotlib, gvar, lsqfit [2], sklearn.metrics. • EM.py: performs the analysis of the single particleenergies and produces the plots stored in the folderresults/effective mass. • SM NM.py: performs the analysis of the energy perparticle in symmetric and neutron matter and pro-duces the plots stored in the folder results/3 scales. • symmetry energy.py: calculates the global sym-metry energy and creates the corresponding plotstored in the folder results/esym esym2. • quadratic symmetry energy.py: calculates thequadratic symmetry energy and creates thecorresponding plot stored in the folder re-sults/esym esym2. • non quadraticities.py calculates the non-quadraticcontributions to the symmetry energy as well as thefinal fit residuals and creates the plots stored in thefolder results/non quadraticities. • crust core.py calculates the crust-core transi-tion and produces the plot in the folder re-sults/crust core. IV. LAUNCH THE PYTHON SCRIPT
The scripts are written in Python 3, and can belaunched from a terminal as: a r X i v : . [ nu c l - t h ] S e p python3 ScriptName.py where “ScriptName.py” is one of the scripts listed inSec. III.where “ScriptName.py” is one of the scripts listed inSec. III.