Setting $\textit{ab initio}$ uncertainties on finite-temperature properties of neutron matter

We present a comprehensive error band on neutron matter properties at finite temperature (finite-T) including uncertainties on the nuclear interaction, the many-body method convergence, and the thermodynamical consistency of the approach. This study provides $\textit{ab initio}$ predictions for finite-T neutron matter employing chiral interactions which are selected on the basis of their performance in both finite nuclei and infinite matter at zero temperature. Since proper theoretical uncertainties at finite-T are still generally lacking, the band provided here represents a significant step forward in setting first-principles constraints on thermal aspects of the nuclear matter equation of state.

In low-energy nuclear physics one is confronted with systems which size ranges from light isotopes (few fm) to compact stars (tens of km). The unifying relation between these apparently different realms, finite nuclei and neutron stars, is the nuclear matter equation of state (EoS) [1]. In this regime, the confinement properties of quantum chromodynamics (QCD) are such that protons and neutrons are the common relevant degrees of freedom at play and, therefore, the knowledge of nuclear forces is key to understand the above phenomena. Deriving nuclear interactions directly from QCD is highly non trivial and nuclear physicists are forced to work the problem around by combining experimental information on finite nuclei, astrophysical observations and theoretical models. Consequently, constraints from different sources try to establish this structural equation for nuclear matter, among which are theoretical calculations of nuclear structure and reactions, systematics of nuclear masses and properties of isolated or merging neutron stars (see Refs. [2,3] for recent reviews).
Ab initio nuclear theory plays a major role in this endeavor. For instance, both predictions of the neutron distribution in 48 Ca [4], and calculations of neutron matter around saturation density [5], have put a window on the size of a neutron star, hence on its EoS, solely based on nuclear physics. Such first-principle predictions are further validated by the fact that the same nuclear forces are employed in finite nuclei and nuclear matter, where they have proved to work equally well [6][7][8][9][10][11]. Thanks to the first detected neutron-star merger event GW170817 [12], outcomes of the analysis combined with neutron matter calculations [13], and also perturbative-QCD [14], are putting challenging constraints on the radius of a neutron star. On top of that, the forthcoming measurements on the mass/radius of neutron stars via the X-ray telescope NICER, will help locate both the high-and low-density regimes of the nuclear EoS [15]. These new developments provide a clear path to nail down the correct relation between pressure and energy density of nuclear matter in the zero temperature (zero-T) limit.
The situation is more uncertain when finite temperature (finite-T) is included in the EoS, condition which has to be considered when studying binary neutron star mergers or matter formed in heavy-ion collisions. Despite several first-principle calculations of nuclear matter at finite-T over the past forty years [16][17][18][19][20][21][22][23][24][25][26][27][28][29], there are only few attempts to include such results in EoSs for astrophysical simulations. Most applications have concentrated on the study of proto-neutron stars [30][31][32] or corecollapse supernovae [33]. For binary neutron star mergers, the most widely used finite-T EoSs are usually based on Skyrme-like models, or on relativistic mean field theory (for a full list of general purpose EoS see Refs. [2,3]). Another approach has been that of describing the thermal part of the EoS as an ideal fluid [34]. Nonetheless, this latter approximation can prove to be quite crude for mimicking thermal properties of the nuclear EoS [35][36][37]. Recently, finite-T EoSs for neutron stars have appeared based on quite different approaches [38][39][40]. However, no detailed uncertainty analysis of nuclear matter properties that can constrain the finite-T nuclear EoS is available to date. This work addresses this latter issue by performing a first-principles study of neutron matter at finite-T, that combines for the first time uncertainties on the nuclear interaction, the many-body approximation and the thermodynamical consistency of the approach. Since ab initio results have shown to be quite constraining at nuclear densities to select acceptable zero-T EoSs from other approaches [5], this study wants to be a first step towards a similar selection based on nuclear physics including finite-T. Furthermore, we present a first complete analysis of the neutron effective mass at finite-T, quantity which has been shown to be crucial in determining the thermal behavior of the nuclear EoS [36,37,41,42].
Within a low-energy non-relativistic framework, we make use of the many-body self-consistent Green's function (SCGF) method [43] to investigate the properties of infinite nuclear matter employing chiral nuclear interactions [44][45][46]. This approach takes into account beyond mean-field effects constructing a fully-dressed nucleon propagator [47]. Via solution of the Dyson equa-tion, G(p, ω) = G 0 (p, ω) + G 0 (p, ω)Σ (p, ω)G(p, ω) , (1) a self-consistent description for the dressed propagator G in terms of single-particle momentum p and energy ω is found, build upon its free version G 0 . The nonperturbative self-energy Σ is constructed at each iterative step in the solution of the Dyson equation through a resummation of intermediate scattering diagrams, within the so called ladder approximation [21,23,24]. Eq. (1) can be recast as a solution for the nucleon spectral function A given by the formula [24]: (2) A self-consistent solution is found when the spectral function which enters the calculation of the self-energy Σ equals the one obtained solving Eq. (2). With this spectral function it is possible to access the total energy per nucleon of the system via the Galitskii-Migdal-Koltun sumrule [48]: (3) ν is the degeneracy of the system, 2 for pure neutron matter (PNM) and 4 for symmetric nuclear matter (SNM); n is the total number density; Ŵ is the expectation value of the three-body operator. The SCGF method is implemented directly at finite-T and it provides in principle a thermodynamically consistent description of the many-body system, i.e. the microscopic and macroscopic (thermodynamical) estimates of physical properties should equal one another [49,50].
We start by analyzing in Fig. 1 the zero-T energy per nucleon of SNM obtained employing Eq. (3). We use two plus three body chiral Hamiltonians which have proven successful in predicting finite nuclei properties: the N2LO sat reproduces nuclear radii and binding energies up to 40 Ca [6,7]; the (EM) potentials [52], with labels λ 2N /Λ 3N being the low-resolution scale on the twobody force and the cutoff on the three-body force, reproduce reasonably well the ground state energies of closedshell nuclei, two-neutron separation energies and 2 + excited states up to 78 Ni [53]; we have added two more interactions, dubbed N3LO/N2LO for the two-/three-body chiral force order, which, a part from being fit to properties of light nuclei as the previous potentials, are also fit to the triton beta decay, with label (Λ) being the nonlocal cutoff on the current [54]. Details for each Hamiltonian can be found in the corresponding Refs. [6], [52], [54]. We obtain a spread of ∼7 MeV in energy range at saturation density n sat = 0.16fm −3 . However predicted saturation points build a Coester-like line, as also seen in Ref. [8]. We use the empirical saturation box in Fig. 1 [51]. Light/grey band is the accepted zero-T uncertainty band (see text for details).
select five interactions out of the original seven, based on being the predicted saturation point consistent with either the density or energy ranges of the box (a light/grey band highlights the chosen interactions). The box in Fig. 1 is given by 12 Skyrme functional calculations, constrained by properties of doubly magic nuclei and ab initio low-density neutron matter [51]. The predicted saturation points lie then in a range of densities n=[0.16-0.18]fm −3 and energies E/A=-[13.5-16.2]MeV.
The selected interactions are used in Fig. 2 to calculate the symmetry energy SYM (highlighted with a band), which is obtained as the difference between the energies of PNM and SNM, also shown in the figure. All SYM curves stand together except for the N2LO sat calculation. This is caused by an extremely soft PNM energy per nucleon. We delimit the predicted symmetry energies from below exploiting the unitary-gas limit (dotted line) corresponding to the conservative choice of Ref. [55]. While all calculations respect this lower bound, the N2LO sat violates it in a region of densities from n sat /2 to 1.25n sat . Furthermore, it is the only case which does not match the comprehensive uncertainty interval given by Oertel et al. [2], which constrains the symmetry energy at n sat employing a large number of theoretical and experimental calculations from different sources. The uncertainty interval at 2n sat provided by Zhang et al. [56], matched by all calculations except for N2LO sat , is extracted constraining an explicitly isospin-dependent parametric EoS with the data analysis of GW170817 [58]. It is interesting to note that, except for a region with densities below ∼0.12fm −3 , it is not possible to reconcile our theoretical predictions with the symmetry energy extracted from the measured 197 Au+ 197 Au reaction data from the FOPI-LAND (lighter/yellow band in Fig. 2) and ASY-EOS (darker/red band) experiments [57]. The experi-  Fig. 1. Dotted line is the conservative unitarygas limit given by Ref. [55]. Intervals at nsat and 2nsat come from Ref. [2] and [56]. Lighter/yellow and darker/red bands define respectively the FOPI/LAND and ASY-EOS results extracted from reaction experiments [57].
mental extracted value of 34MeV at n sat stands higher than the entire ab initio uncertainty band we provide, with growing departure as density increases [57]. The explicit inclusion of the ∆ degree of freedom could help improve this issue [59].
In spite of the apparently poor performance of N2LO sat , we choose to use the entire band of the five interactions to calculate the zero-T pressure in PNM. This choice is also based on the outstanding performance of N2LO sat in finite nuclei [4,7,9]. Figure 3 compares our results with the constrained bands at a 90% and 50% confidence level as given by the analysis of GW170817 [58]. All calculations are very close to the LIGO/Virgo bands at low densities (results from Ref. [58] start at ∼ n=0.06fm −3 ), except for the N2LO sat which remains very soft. However, as it appears clearly from Fig. 3, the N2LO sat pressure grows fast approaching the other calculations around n sat . The ab initio calculations then stay within the 50% confidence level band all the way up to ∼ 2n sat . This comparison underlines the fact that, even though the PNM energy predicted could be low, as in N2LO sat , its derivative to obtain the pressure can be in any case quite steep.
Based on all of the above considerations, we select two final chiral interactions to generate an uncertainty band to study finite-T properties of neutron matter: the upper limit given by the N2LO sat and lower limit by the 2.0/2.0(EM). We restrict the upper limit of the grey band  Fig. 2, legend follows that in Fig. 1. Small and big dotted lines delimit the LIGO/Virgo 90% and 50% confidence level bands on the PNM pressure as obtained from the analysis of GW170817 [58].
in Fig.1 to N2LO sat on the basis of the closest predictions to the SNM empirical saturation point, and given that the error at 2n sat encompasses practically all the five previously selected ones, so we are sure to keep the uncertainty at high densities as conservative as possible (this applies also for symmetry energy and pressure in Figs. 2-3). Furthermore, our choice is corroborated by a reliable prediction of the SNM liquid-gas phase transition obtained using these interactions [28]. We present in Fig. 4 specific PNM properties at T=30 MeV, chosen as a representative temperature, employing the above mentioned chiral interactions. In panel (a) we show the pressure as a function of number density. To understand the uncertainty we have on the many-body calculation and test the method convergence, we present three different approximations: self-consistent Hatree-Fock (SCHF), self-consistent second order (SC2O) and SCGF, where the self-energy is truncated respectively at first-, second-and all-orders (ladder ). Furthermore, we test the thermodynamical consistency of our calculations by comparing the SCGF pressure obtained via the microscopic chemical potentialμ, i.e.P = n(μ − F/A), with respect to the pressure obtained thermodynamically via derivative of the free-energy per nucleon F/A, i.e. P = n 2 ∂F/A ∂n , dubbed SCGF (therm) (see Ref. [28] for further details). On one hand we see how, for both chiral interactions, the SC2O calculations stand practically on top of the full SCGF ones, meaning that self-energy terms beyond second-order are quite small for these interactions. On the other hand, the differences arising between the microscopic and macroscopic derivation of the pressure, SCGF vs SCGF(therm), are more visible. Thermodynamical consistency holds when only two-body forces are considered [24], so this issue is related to the inclusion of three-body forces and specifically to the approximation with which we calculate Ŵ in Eq. (3) [28,60]. It is instructive to see how even the SCHF first order calculations stand very close to the full SCGF ones, providing at 2n sat a band of ∼1.5MeV for the full many-body method uncertainty, compared to ∼6 MeV coming from the chiral interaction error band.
Panel (b) of Fig. 4 displays the neutron effective mass at T=30 MeV calculated at the Fermi momentum p = p F and Fermi energy ω = ε(p F ) for each density. Contrary to the pressure, a big discrepancy is visible here between the SCHF and the SC2O/SCGF calculations. A part from the differences in the self-energy truncations, this behavior is to be ascribed also to the fact that at the SCHF level we only have one term for the effective mass, m * /m = m k = 1 + m p ∂ReΣ(p,ε(p)) ∂p −1 , while for both SC2O and SCGF we also have a contribution from the m ω = 1 − ∂ReΣ(p,ω) ∂ω , which leads to a total effective mass of m * /m = m k m ω . These two quantities, m k and m ω , incorporate the properties of the varying singlenucleon self-energy in terms of momentum and energy respectively, measuring its non locality either in space or time [61]. Similar to the pressure in panel (a), the SC2O m * is already a good approximation of the full SCGF results. The rising of the effective mass with density after reaching a certain minima is caused by the inclusion of three-body forces [62][63][64].
To understand more in depth the behavior of the effective mass, we present in panel (c) and (d) of Fig. 4 the single-particle potential as a function of momentum and the real self-energy as a function of energy, which contribute to the calculation of the m k and m ω respectively. The single-particle potential in panel (c) corresponds to an "on-shell" real self-energy, ReΣ(p, ε(p)), obtained solving a self-consistent equation for the single-particle energy [24,65]. A more repulsive single-particle potential and a steeper behavior for the SCHF case is what causes the m k to become smaller, and consequently obtain a smaller effective mass with respect to the SC2O/SCGF cases (see panel (b)). In panel (d) of Fig. 4 we plot the self-energy as a function of energy calculated at p F for n sat . The blue dots in the inset show the points, ω = ε(p F ), where the derivative of the self-energy is performed for the calculation of the m ω , highlighting quite a different steepness for the two interactions. The inclusion of the m ω is found to be fundamental to reproduce the behavior of thermal effects in finite-T EoSs [37].
To conclude, we employed state-of-the-art chiral interactions, proven successful in describing finite nuclei and zero-T infinite matter, to set a global uncertainty band on finite-T properties of neutron matter. We provided a comprehensive error analysis that accounts for uncertainties in the interaction, the many-body method truncation and the thermodynamical consistency of the approach. The major uncertainty on the PNM finite-T pressure is related to the chiral interaction, being this at 2n sat four times larger than the error associated to the many-body method. On the contrary, this behavior is reversed for microscopic properties, such as the effective mass or the single-particle potential, where the full many-body error is more than twice the chiral interaction uncertainty band at 2n sat . This underlines the fact that beyond first-order calculations of the nucleon effective mass are mandatory, especially in view of recent studies which show how parametrized functions of this quantity can mimick the thermal part of the EoS [37]. Furthermore, the value of the effective mass has been shown to be crucial for the onset of explosion in core-collapse supernovae [41,42]. This study of ab initio uncertainties on finite-T properties of neutron matter represents a first fundamental step to reliably constrain the nuclear EoS employed in astrophysical simulations of merging neutron stars [35,[66][67][68][69].
The Author is grateful to Carlo Barbieri for a thorough reading of the manuscript and for much needed comments. Fundamental discussions with Arnau Rios on the numerical calculations are acknowledged. Calculations for this research were conducted on the Lichtenberg high-performance computer of the TU Darmstadt.