The Equation of State of Neutron-Rich Matter at Fourth Order of Chiral Effective Field Theory and the Radius of a Medium-Mass Neutron Star

We report neutron star predictions based on our most recent equations of state. These are derived from chiral effective field theory, which allows for a systematic development of nuclear forces, order by order. We utilize high-quality two-nucleon interactions and include all three-nucleon forces up to fourth order in the chiral expansion. Our ab initio predictions are restricted to the domain of applicability of chiral effective field theory. However, stellar matter in the interior of neutron stars can be up to several times denser than normal nuclear matter at saturation, and its composition is essentially unknown. Following established practices, we extend our microscopic predictions to higher densities matching piecewise polytropes. The radius of the average-size neutron star, about 1.4 solar masses, is sensitive to the pressure at normal densities, and thus it is suitable to constrain ab initio theories of the equation of state. For this reason, we focus on the radius of medium-mass stars. We compare our results with other theoretical predictions and recent constraints.


I. INTRODUCTION
The equation of state (EoS) of neutron-rich matter is at the forefront of nuclear astrophysics because of its role in shaping the properties of neutron stars. Recently, interest in compact stars has increased considerably as we have entered the "multi-messenger era" of astrophysical observations. The recent GW170817 neutron star merger event has yielded new and independent constraints on the radius of the canonical mass neutron star [1,2]. Astronomy with gravitational waves provides additional opportunities to explore these exotic systems and other yet unknown regimes in the Cosmos.
The structure of a neutron star probes a very large range of densities, from the density of iron up to several times the nuclear matter saturation density, and thus no theory of hadrons can be considered reliable if extended to those regions. On the other hand, contemporary ab initio theories of nuclear and neutron matter at normal densities can be taken as the baseline for any extension or extrapolation method, which will unavoidably involve phenomenology. We recall that the radius of a 1.4 M is sensitive to the pressure at normal densities, see Ref. [3], and thus it is a suitable constraint for microscopic theories of the EoS at those densities where they are reliable. Chiral effective field theory (EFT) [4,5] provides a path to a consistent development of nuclear forces. Symmetries relevant to low-energy QCD are incorporated in the theory, in particular chiral symmetry. Thus, although the degrees of freedom are pions and nucleons instead of quarks and gluons, there exists a solid connection with the fundamental theory of strong interactions through its symmetries.
Chiral EFT employs a power counting scheme in which the progression of two-and many-nucleon forces is constructed following a clear and systematic hierarchy. This allows for the inclusion of all three-nucleon forces (3NFs) which appear at a given order, thus eliminating the inconsistencies which are unavoidable when adopting mesontheoretic or phenomenological forces. Finally, it provides a clear method for controlling the truncation error on an order-by-order basis. In the simplest approach, the latter can be expressed as the difference between the quantity computed at a given order of the chiral expansion and the one obtained at the next order. We will revisit this point in Section IV A.
We will report neutron star predictions based on our most recent EoS, which includes all subleading 3NFs [6,7]. For the reasons stated above regarding the limitations of chiral EFT, we focus on the average-mass neutron star, rather than the maximum mass of a sequence, as the latter is much more sensitive to the polytropic extensions we perform in order to complete the EoS. We also discuss the proton fraction present in our β-stable EoS, because of the radius of the canonical-mass neutron star predicted from the set of EoS applied in Ref. [49] is predicted to be in the range 10.45-12.66 km. From a variety of techniques, based on experimentally determined quantities correlated to symmetry energy parameters, the radius is determined to be between 10.7 and 13.1 km [49][50][51][52], while using a range of theoretical models a limit of 9.7 to 13.9 km is obtained [49,53,54]. Recent surveys of neutron star physics and theoretical approaches include Refs. [55][56][57]. Exotic matter in stars is addressed, for instance, in Ref. [58].
On theoretical grounds, the largest mass was predicted to be 3.2 M [59], based on only three assumptions: (1) General Relativity is the appropriate theory to describe these massive stars; (2) the EoS is constrained by Le Chatelier's principle (∂P /∂ ≥ 0); (3) the causality condition, which constrains the speed of sound in dense matter to remain below the speed of light. While such massive neutron star may be theoretically possible, none has been observed in this mass range.
It is interesting to note the small range of values for the radius across the mass range of neutron stars. Heavier neutron stars have larger central densities and thus the star becomes comparatively more compact, resulting in a very narrow mass-radius range, in contrast to main-sequence stars whose masses and radii span several orders of magnitude.

III. DESCRIPTION OF THE CALCULATION
The EoS for neutron and symmetric matter are obtained at the leading-order in the hole-line expansion-namely, via a non-perturbative calculation of the particle-particle ladder. As pointed out in Refs. [6,7], the third-order hole-hole diagram was considered and found to be very small at normal density [60]. In comparison, the particle-hole diagram is larger, but its contribution is still relatively minor-we estimate an uncertainty in the order of ±1 MeV on the potential energy per particle at saturation. The single-particle potentials are computed self-consistently with the G-matrix, employing a continuous spectrum.
Next, we proceed to describe the input two-nucleon forces (2NFs) and 3NFs.

A. The Two-Nucleon Force
The 2NF we apply are from Ref. [61], a family of high-quality potentials from leading order (LO) to fifth order (N 4 LO) of the chiral EFT. At fifth order, the nucleon-nucleon (N N ) data below pion production threshold are reproduced with the excellent precision of χ 2 /datum = 1.15.
The interactions in this set are more internally consistent than those from the previous generation [62]. Furthermore, the long-range part of these potentials is tightly constrained by the πN low-energy constants (LECs) from the Roy-Steiner analysis of Ref. [63]. This analysis is sufficiently accurate to render errors in the πN LECs essentially negligible for the purpose of quantifying the uncertainty.
To suppress high-momentum components in the potential when solving the non-perturbative Lippmann-Schwinger equation, one must apply a regulator, for which we choose the non-local form: expressed in terms of the final (p ≡ | p |) and initial (p ≡ | p |) nucleon momenta in their center-of-mass system. With the choice Λ = 450 MeV-which we maintain throughout this paper-the potentials are soft according to the Weinberg eigenvalue analysis of Ref. [64] and the perturbative calculations of infinite matter from Ref. [65].

B. The Three-Nucleon Force
In the framework of the ∆-less chiral EFT (which we apply), the first occurrence of 3NF is seen at the third order. The leading 3NF consists of three components [66]: the long-range two-pion-exchange (2PE) graph, which depends on LECs c 1 , c 3 , and c 4 , the medium-range one-pion-exchange (1PE) diagram, carrying the LEC c D , and a short-range contact term, containing the LEC c E . We recall that the terms depending on c 4 , c D , and c E do not contribute in neutron matter [67].
In infinite matter, it is possible to construct approximate expressions for the 3NF as density-dependent effective twonucleon interactions as derived in Refs. [68,69]. These can be written in terms of the well-known non-relativistic twobody nuclear force operators and, thus, can be easily implemented in the N N partial wave formalism for the G-matrix, which leads to the EoS. One must be careful to avoid overcounting in the 3NF when the latter is represented as densitydependent potentials, which is accomplished by including the appropriate combinatoric factor in the calculation of the energy per particle, as done in Ref. [6,7].  The topologies of the effective density-dependent two-nucleon interactions originate from the corresponding chiral 3NF. At N 2 LO, there are six one-loop topologies. Of those, three come from the 2PE graph of the chiral 3NF and contain the LECs c 1,3,4 , which also appear in the 2PE part of the N N interaction. Of the remaining three one-loop topologies, two originate from the 1PE diagram of the 3NF, and depend on the LEC c D . The last one-loop diagram is the 3NF contact term, with LEC c E .
We also include the subleading (N 3 LO) 3NF, derived in Ref. [70,71]. References [65,[72][73][74] report applications of the subleading 3NF in some many-body systems. The long-range part of the 3NF at N 3 LO includes the 2PE topology, which is the longest-range contribution, the two-pion-one-pion exchange (2P1PE) topology, and the ring topology, which represents a pion being absorbed and reemitted from each of the three nucleons. Again, in-medium N N potentials can be obtained from these topologies. For the long-range subleading 3NFs, the expressions are given in Ref. [75] for symmetric nuclear matter (SNM) and in Ref. [76] for neutron matter (NM). The short-range subleading 3NF includes the following topologies: the one-pion-exchange-contact (1P-contact), which ends up giving a vanishing net contribution, the two-pion-exchange-contact (2P-contact), and relativistic corrections, which depend on the C S and the C T LECs of the 2NF and are proportional to 1/M , where M is the nucleon mass. We include relativistic corrections as well and find them to be very small (less than one MeV). The expressions for the in-medium N N potentials generated by the short-range subleading 3NFs are taken from Ref. [77] for SNM and Ref. [78] for NM. Table I displays the LECs we use, which are taken from Ref. [65]. An important point to highlight: the largest part of the subleading two-pion-exchange 3NF has a very similar formal structure to the leading 2PE [79] and thus a large part of the subleading two-pion-exchange 3NF can be effectively accounted for with a shift of the LECs. For this reason, when the subleading 3NFs are included, we replace the c 1 , c 3 , and c 4 LECs shown in Table I  In Figure 1, we show the EoS in NM over four orders, from LO to N 3 LO [6]. Large variations at low orders are of course not surprising, nor is the remarkable impact of the leading 3NF at N 2 LO. The transition to fourth order brings in a slight increase in attraction, as was found from other EFT-based predictions [73]. The overall convergence pattern is encouraging.
In Figure 2, an analogous presentation is provided for SNM. Similar considerations apply with regards to the order-by-order convergence pattern and the 3NF "signature" as in NM.

C. Equation of State for Stellar Matter
In this section we outline the main steps to construct the EoS for stellar matter in β-equilibrium. We define the total energy per baryon as: where Y n,p stands for the neutron or proton fraction. The last term accounts for the baryon rest energy, while e e/µ are the energies (per baryon) of the electrons and muons, respectively. All terms are functions of density. The energy density ( i ), pressure (p i ), and number density (ρ i ) for each particle species, i, at a given Fermi momentum, (k F ) i , can be expressed as: where γ is the spin/isospin degeneracy factor. The energy per volume is obviously related to the energy per particle: = ρ e(ρ) .
The Fermi energy for each species, i, or chemical potential, is given by The fractions of each particle species, is related to the chemical potential through: From the condition of energy minimization, subjected to the constraints of conserved nucleon density and global charge neutrality, one can obtain all particle fractions.

D. Polytropic extrapolation
Chiral predictions have a limited domain of validity, which, in the previous section, we estimated to be about twice the saturation density. The densities within neutron stars can reach five to six times saturation density, and therefore an appropriate method for extrapolating the EoS to these densities must be employed. To accomplish this, we express the high density pressure through polytropes [80]: where α is chosen such as to ensure continuity at the matching density. A comment is in place: while continuity of the pressure is of course preserved, additional considerations are necessary to ensure continuity of the derivative. The latter would be essential to implement thermodynamic consistency of the piecewise EoS, which is beyond our present scope. Note, further, that the presence of discontinuities in the polytropic index is not unusual for the purpose of describing the global features of the star [80]. Following Ref. [3], we match piecewise polytropes to the ab initio predictions as explained next. The microscopic predictions reach a Fermi momentum of 1.6 fm −1 , which corresponds to 2.016 fm −1 in pure neutron matter at the same density, ρ = 0.277 fm −3 . The standard practice is to ensure that the characteristic momentum of the system, p, divided by the cutoff, Λ, is a reasonable expansion parameter. Taking p to be the average momentum in a free Fermi gas of neutrons at the highest density we consider, we obtain a value of 68% for p/Λ, which is a pessimistic estimate, since the average momentum in β-stable matter is smaller than in pure NM. Having chosen the matching density, ρ 1 , we join the pressure predictions with polytropes of different adiabatic index, ranging from 1.5 to 4.5. This range is chosen following guidelines from the literature, in particular Ref. [80], where constraints on phenomenologically parameterized neutron-star equations of state are investigated. To simulate a (likely) scenario where the pressure displays different slopes in different density regimes, we define a second matching density, ρ 2 , approximately equal to 2ρ 1 , at which point a set of polytropes covering the same range of Γ is attached to each of the previous polytropes, yielding a total of 49 possible combinations. This is illustrated in Figure 3. It is important to emphasize that high-density EoS extrapolations are not meant to be a replacement for microscopic theoretical predictions [3] which, at this time, are not feasible at super-high densities. Instead, the spreading of the high-density pressure values from the piecewise variation of the polytrope index allows to probe the sensitivity of lower-density predictions to the much larger uncertainty at high density.
To construct a physical EoS for high densities, we must apply additional constraints. One is the causality limit, which imposes the speed of sound in matter to be less than the speed of light. In terms of P ( ), the causality condition reads: Note that the constraint on the speed of sound is strictly valid only in the absence of dispersion or absorption in stellar matter [81]. Nevertheless, imposing the causality constraint is standard practice when constructing neutron star EoS and we will apply it in this work [32,59,80,81].
Additionally we will only consider polytropes, which can support a maximum mass of at least 2.01 M , to be consistent with the lower limit of the (2.08 ± 0.07) M observation reported in Ref. [82] for the J0740 + 6620 pulsar along with a radius estimate of (12.35 ± 0.75) km. To complete the EoS on the low-density side, we attach a crustal EoS [83,84].

E. Mass-radius relation
With the EoS available over a full range of densities, we move to the mass-radius relation in a neutron star. In this section we will briefly review the relativistic equations for hydrostatic equilibrium, the TOV equations [13,14] and how the mass-radius relationship emerges from them for a given input EoS.
The TOV equation describes a spherically symmetric inertial massive object composed of a perfect fluid in hydrostatic equilibrium. The equation relates the pressure within the star to the mass-energy density as functions of the radial distance from the star's center: A spherical shell of material is related to the energy density at a distance r from the star's center by: The star's gravitational mass (M ) is determined from the radius (R) and the mass-energy density ( (r)): Since the pressure and energy-density are functions of density, for a fixed central density the mass-radius of the star can be determined by Equations (15) and (14). Equation (14) can be integrated numerically by summing over shells of fixed width at incremented distance from the star's center so as to evaluate the total pressure as a function of radial distance. Equation (15) can be integrated in the same fashion, simultaneously, to determine the mass contained within each spherical shell. To accomplish this, we employed the fourth-order Runge-Kutta method. The radial distance at which the pressure effectively vanishes corresponds to the star's radius. Then, Equation (16) provides the total mass enclosed within such radius.

A. Chiral Uncertainty
First, some comments on the estimation of the chiral error. As pointed out in Section III A, errors in the πN LECs are small enough to be neglected when quantifying the combined uncertainty. The truncation error is of course central to the philosophy and the application of chiral EFT, being the indicator of the quality of the convergence pattern.
In the most basic approach, one may argue that, if observable X has been determined at order n and at order n + 1, a reasonable estimate of the truncation error is simply the difference between the value of X at order n and the one at the next order: To estimate the uncertainty at the highest available order, one needs additional considerations. We follow the prescription of Ref. [85]. First, one needs to identify the typical momentum or energy scale for the system being considered, say, p. Defining Q as the largest between p Λ b and mπ Λ b , where Λ b is the breakdown scale of the chiral EFT, about 600 MeV [85]. The truncation error for the observable X at N 3 LO is then defined as: We take p to be the average momentum of a Fermi gas of neutrons, 3 5 k n F .

B. Results
In Figure 4, we show the proton fractions we obtain from leading to fourth order of chiral perturbation theory. Electron fractions are not shown to avoid excessive crowding-prior to the onset of muons, electron and proton fractions are of course equal, and remain close to each other, as it can be inferred from the small values of the muon fractions. Our proton fractions are very small, consistent with the relatively soft symmetry energy, see Figure 5, which brings up the issue of direct Urca processes-the most effective neutrino emission mechanism. Our predictions are far from the direct Urca threshold of approximately 11% around normal density. Instead, we find for the N 3 LO proton fraction a normal density a value of 0.046 ± 0.0035.
To offer the reader a broader overview, we show in Table II the value of the symmetry energy and the slope parameter L with their uncertainty, where L is defined as These values were shown in Ref. [6] and compared to recent constraints. We proceed with pressure predictions. Figure 3 displays the pressure in stellar matter as a function of the number density. The various curves span the range of acceptable combinations of polytropes, as explained in Section III D. In Table III we present predictions for the radius, central density, and speed of sound for M = 1.4 M for those combinations. Several comments are in place. First, we note that, overall, the radius is not very sensitive to the Γ 1 , Γ 2 variations. For values of Γ 1 on the low end of the range, acceptable combinations require values of Γ 2 on the higher end-understandable in terms of maximum mass constraints. Further, for a fixed value of Γ 1 , the sensitivity of the radius to variations in Γ 2 is essentially negligible, and of course it vanishes when the central density is below the second matching density. One more comment regarding causality: While all EoS are causal at the central densities displayed in Table III, some may not be so at the central density of the maximum mass of the sequence. In such cases, the M(R) correlation is typically truncated at the causality limit, retaining the EoS at the lower densities.
The mean value and standard deviation are (R 1.4 = 11.96 ± 0.58) km. The same procedure is applied at the lower orders-LO to N 2 LO-that is, the EoS at each order is extended with polytropes and the mean value of the radius is calculated. With radius predictions available from leading to fourth order, we determine the truncation error via the prescription in Equation (18), where we take p to be the neutron Fermi momentum at saturation density. This choice is reasonable because the average density of a neutron star is comparable to saturation density and because of the strong sensitivity of the radius to the normal density region. From Equation (18) we obtain an uncertainty of ± 0.54 km. Combining the truncation and extrapolation uncertainties in quadrature, we state our estimate of the radius as: in excellent agreement with the LIGO/Virgo range of 11.1 to 13.4 km [47]. Our result is within the range generally found with EoS based on chiral EFT, which is 10 km to 14 km [86,87], accounting for additional theoretical uncertainties, such as those originating from the choice of the many-body method and the implementation of the 3NF [88][89][90][91][92]. Some sensitivity of R 1.4 to the matching density was found [93,94]. Moving the matching density from ρ 0 to 2ρ 0 changed the range to (9.4-12.3) km [94] and to (10.3-12.9) km [93]. In the present analysis, the first matching point is determined by the highest density we reach out with the EFT calculations-a natural matching point. As for the second matching density, Table III shows that the details of the extension at the higher densities has only a minor impact on R 1.4 .
Closely related to the star radius is the tidal deformability, which is a measure of how easily the star deforms under an external tidal field. Therefore, a smaller compact star will have smaller tidal deformability as compared to a larger, less dense star. The tidal deformability is defined as the ratio of the induced quadrupole Q ij to the perturbing tidal field E ij , which explains the sensitivity to the radius [95]. The dimensionless tidal deformability in the mass range that is relevant for GW170817 (and for the present analysis)-1.1M ≤ M ≤ 1.6M -is found to be proportional to is the (dimensionless) star compactness. The dimensionless tidal deformability can then be expressed as: where a = 0.0093 ± 0.0007 [96]. Using our estimate for the radius as in Equation (20), we obtain from Equation (22) a range of Λ between 247 and 474, in excellent agreement with the Bayesian analysis of Ref. [93], where the limits of 68% credibility are found to be 249 and 465, with Λ = 379 the most probable value. Recent analyses of the neutron star radius include Refs. [97,98]. Our calculations differ from those in several ways, most importantly, the derivation of the microscopic part of the EoS is based on local chiral interactions at N 2 LO-by contemporary standards, high-precision description of N N data requires the construction of potentials at the fourth order. Furthermore, all subleading 3NF at the same order have been available for some time. For these reasons, modern chiral EFT-based nuclear forces should be consistently at N 3 LO. Other differences include the many-body theory-local interactions are constructed to be used in QMC calculations-and the derivation of the β-stable EoS,

V. SUMMARY AND CONCLUSIONS
A fully microscopic EoS up to central densities of the most massive stars-potentially involving non-nucleonic degrees of freedom and phase transitions-is not within reach. Nevertheless, neutron stars are powerful natural laboratories for constraining theories of the EoS. One must be mindful about the theory's limitations and the best ways to extract useful information from the observational constraints. Concerning the latter, there is no doubt that The golden age of neutron-star physics has arrived [101].
Recently, we developed EoS for NM and SNM based on high-quality 2NF at N 3 LO and including all subleading 3NF. These were used to construct the symmetry energy and the EoS in β-equilibrated matter, from which proton and lepton fractions are easily extracted. The stellar matter EoS was then extended to densities inaccessible to chiral EFT by matching it with piecewise polytropes of different adiabatic index. From those combinations, we dropped the EoS which did not satisfy the maximum mass constraint. The causality condition is also applied.
Constraints on the radius of a medium-mass neutron star, R 1.4 , are becoming more stringent, with the current uncertainty reported at about 2 km. Furthermore, R 1.4 is known to be sensitive to the pressure in neutron-rich matter near normal densities, accessible to modern effective field theories of nuclear forces. For these reasons, we focused on predicting R 1.4 with proper uncertainty quantification. From reports in the literature, our predicted range would increase by about 1 km on either side when additional theoretical uncertainties are included.
Based on our analysis in Section IV B, we are confident that the estimate given in Equation (20), approximately (12 ± 1) km, is characteristic of EFT predictions based on high-quality 2NF and properly calibrated (leading and subleading) 3NF [102]. The range currently cited for chiral EFT-based predictions of R 1.4 is between 10 km and 14 km, accounting for additional theoretical uncertainties. In fact, it is interesting to notice that the extensive analysis from Ref. [93], where 300,000 possible EoS were generated, provides a range for R 1.4 between 10.0 and 12.7 km, with 12.0 km being the most probable value. We recall that the outcome of PREX II gives 13.33 km as the lower limit for the radius, which is problematic to reconcile with a multitude of microscopic predictions [8].
In conclusion, we reiterate that gravitational wave astronomy offers new exciting opportunities for nuclear astrophysics. Even though chiral EFT cannot reach out to the extreme-density and yet unknown regimes at the core of these remarkable stars, continuously improved ab initio calculations of the nuclear EoS are an essential foundation for interpreting current and future observations in terms of microscopic nuclear forces.