In the previous section, we mentioned the linked-cluster perturbation series for the energy of a many-body system [1–6]. To facilitate convergence (otherwise problematic in view of the strong repulsive core of the NN force), the linked-cluster expansion for the energy per particle in nuclear matter [3] is written in terms of the reaction matrix or “G-matrix,” which itself is solution of the Bethe-Goldstone equation. Schematically, the Bethe-Goldstone equation can be written as

where V is the NN potential, Q is the Pauli operator, and E₀ the starting energy of the two nucleons. The second term on the RHS of Equation (1) represents the infinite ladder sum which builds short-range correlations (SRC) into the wave function. The correlated (ψ) and the uncorrelated (ϕ) wave functions satisfy

from which one can write

At large distances, the correlated wave function is expected to approach the uncorrelated one (a behavior known as the “healing” property), whereas the two can be very different at short range. Hence, the difference between the correlated and the uncorrelated wave functions, or “defect function” f = ψ − ϕ, can be associated to the degree of SRC.

Usually, its momentum-dependent Bessel transform is considered instead, so as to bring out the dependence on specific partial waves. For each angular momentum state [5], we then have

where the angle-averaging has been applied to the Pauli operator, Q̄. Equation (4) is related to the probability of exciting two nucleons having relative momentum k₀ and relative orbital angular momentum L to a state with relative momentum k and relative orbital angular momentum L′. The integral of the probability amplitude squared is known as the “wound integral” and defined, for each partial wave at some density ρ, as

Thus, both f and κ contain information on correlations present in the wave function and the G-matrix. The degree of SRC has been traditionally associated with the “strength” of a given potential, as indicated, for instance, by the deuteron D-state probability [26].

The topic of SRC deserves a review by itself and will not be covered here. However, we have taken the opportunity to recall how one may obtain, through Equations (4–5), some information about SRC in nuclear matter. The latter is complementary to studies of SRC in nuclei, which are currently the object of intense experimental investigations through high momentum-transfer (inclusive or exclusive) electron scattering measurements. (For a review on this topic, see [27] and references therein). Two-nucleon dynamics at short distances is mostly determined by the presence of short-range repulsion in the two-nucleon force, which is one of the reasons why a mean-field picture of the nucleus has strong limitations. Short-range correlations, particularly two-nucleon correlations, are therefore fundamentally important and open intriguing questions concerning momentum distributions in nuclei as a tool to probe the off-shell nature of the NN potential. For a recent work of our group on SRC in A=2,3 nuclei see [28].

Back to Equation (1), we solve it self-consistently to obtain the G-matrix together with the single-particle potential, which we define for (anti-symmetrized) states below and above the Fermi level according to the so-called “continuous choice”:

The starting energy is written as

in terms of on-shell single-particle energies

where T is the kinetic energy. The average energy per particle in nuclear matter is then obtained from

Equation (9) as a function of density is the nuclear EoS. Next, we will address how the NN potential V in Equation (1) is constructed.

It is our philosophy that constructing the EoS microscopically from state-of-the-art few-body interactions is the right way to gain insight into effective nuclear forces in the medium. High-precision meson-theoretic interactions [29–31] are often utilized in contemporary calculations of nuclear matter, structure, and reactions. However, in the meson-theoretic approach it is difficult, if not impossible, to maintain a strong connection between the 3NF, or more generally the A-nucleon forces with A > 2, and the associated 2NF [32]. On the other hand, χEFT [18, 20, 22] provides a systematic way to construct nuclear many-body forces consistently [17] with two-body forces, as well as to assess theoretical uncertainties through a systematic expansion controlled by a counting scheme [15]. Furthermore, and perhaps most importantly, χEFT maintains consistency with the symmetries and symmetry breaking pattern of the fundamental theory of strong interactions, QCD.

Because of the strengths described above, χEFT has become the authoritative approach for developing nuclear forces. Applications include few-nucleon reactions [33–38], nuclear structure, especially of light- and medium-mass nuclei [39–55], cold infinite matter [22, 53, 56–64], infinite matter at finite temperature [65, 66], and various aspects of nuclear dynamics [67–73].

In regard to the connection between nuclear matter properties and finite nuclei, it is interesting to point out a persistent problem encountered in structure calculations and related to the bulk properties of medium-mass nuclei. Typically charge radii are underpredicted [74] while the opposite is true for binding energies [75]. Including the desired properties of medium-mass nuclei directly into the fitting protocol for the low-energy constants (LECs) which parametrize short-distance physics in chiral nuclear forces has resulted in improved predictions [76]. However, for a truly microscopic approach the 2NF should be constrained by two-nucleon data and the 3NF by three-nucleon data, without additional adjustments. Applications to A > 3 systems would then be actual predictions, although they may carry substantial uncertainties.

Two recent studies [54, 55] provide indications for how the overbinding problem may be overcome. In these studies, a rather soft nucleon-nucleon (NN) potential (due to renormalization group evolution) along with 3NFs fitted to the binding energy of ³H and the charge radius of ⁴He were employed to calculate the ground-state properties of closed shell nuclei from ⁴He to the light Tin isotopes [54, 55]. Predictions of the ground-state energies were accurate, whereas the radii were somewhat underpredicted, although still in fairly good agreement with experiment. These features can be linked to the good nuclear matter saturation properties of the employed 2NF + 3NF combination [57]. In the above example, the 2NF was soft and alone would lead to substantial overbinding in nuclear matter, whereas the addition of a repulsive 3NF contribution leads to a much better description of the nuclear matter saturation point [57]. As we mentioned earlier, the first quantitative explanation of nuclear matter saturation was achieved in this way within the framework of Dirac-Brueckner-Hartree-Fock theory [12, 14, 77–79]. As an alternative, one could begin with a relatively repulsive 2NF and then add an attractive, density-dependent 3NF contribution. An example of such combination is provided by the Argonne v₁₈ (AV18) 2NF [31] together with the Urbana IX 3NF [80]. However, in this way satisfactory predictions for both the nuclear matter saturation energy and density cannot be obtained [81] and the binding energies of medium-mass nuclei are seriously underpredicted [82]. A similar scenario presents itself when the AV18 2NF is used in combination with the Illinois-7 3NF [82, 83]. Efforts to treat the 3NF microscopically were reported in Zuo et al. [84] and Li et al. [85]. In Li et al. [85], in particular, a 3NF including the Δ, Roper, and nucleon-antinucleon excitations was proposed, based on the Bonn [86] and the Nijmegen [30] one-boson-exchange potentials.

The predictions reviewed in this work are based on the high-quality soft chiral NN potentials from leading order to fifth order of the chiral expansion constructed in Entem et al. [87]. More details are provided below.

The properties of isospin-polarized matter have relevance for a number of open questions in nuclear physics and nuclear astrophysics. For instance, the existance of the neutron drip lines, the thickness of neutron skins, and the properties of neutron stars all have in common a strong sensitivity to the EoS of neutron-rich matter. The symmetry energy determines to a good approximation the energy per particle in homogeneous nuclear matter with any degree of isospin asymmetry (cf. Equation 14 below). The symmetry energy and its density dependence are therefore a key focus of contemporary theoretical and experimental investigations, and much effort has been devoted to identifying nuclear observables which correlate with this important property of infinite matter [25, 45, 53, 56–64, 110–120].

The isospin asymmetry parameter is a measure of the relative densities of neutrons and protons and is defined as

where ρ_n and ρ_p are the neutron and proton densities. It is useful to write the energy per particle in isospin asymmetric matter at some density as an expansion with respect to α:

Frequently, the expansion above is truncated at the term quadratic in α, resulting in the popular parabolic approximation:

where e₀(ρ) = e(ρ, α = 0). [Equation (14) has been verified to be valid up to fairly high densities [119].] Within the assumption of Equation (14), the symmetry energy, e_sym, is the difference between the energy per particle in neutron matter and the one in symmetric matter. We can expand the symmetry energy about the saturation density, ρ₀,

The slope parameter, L, is defined as

and therefore is a measure of the density dependence of the symmetry energy around saturation density. We recall that L is an important quantity because of its significance for the skin thickness in neutron-rich nuclei. Experiments which plan to measure the neutron radius of ²⁰⁸Pb and ⁴⁸Ca using electroweak probes, such as PREX II [121] and CREX [122], respectively, are expected to provide accurate measurements of the neutron skin. As a consequence, one hopes for reliable constraints on the symmetry pressure, clearly related to the slope parameter (see section 3.3.1 below). Also, the radius of the average-mass neutron star is known to be sensitive to the pressure in neutron matter at normal density, P0NM, which is simply related to L (for fixed ρ₀) due to the vanishing of the density derivative of e(ρ, α = 0) at saturation. That is:

The reader is referred, for instance, to Sammarruca and Millerson [123] and the comprehensive list of citations therein.

The isovector incompressibility, K_sym, is associated with the next higher-order derivative, that is, it measures the curvature of the symmetry energy at saturation density. It is defined as:

Correlations between L and both K_sym and the symmetry energy at saturation, e_sym(ρ₀) [124–126] have been examined. Predictions for the isovector incompressibility carry large uncertainty, as is the case for the isoscalar one. Attempts to constrain the second derivative of the symmetry energy (that is, its curvature) are discussed in Vidaña et al. [127], Ducoin et al. [128], and Santos et al. [129].

Figure 6 displays the energy per particle in isospin asymmetric matter as a function of density and for increasing degree of asymmetry, cf. Equation (14), for one selected order and cutoff [130].

As we already noted in conjunction with Figure 2, the saturation properties of the chiral interactions we are considering are different from one another, with the saturation density varying between about 0.16 and 0.20 fm⁻³. Clearly, this will impact the expansion parameters contained in Equation (15), see definitions in Equations (16–18), differently than if the derivative were evaluated, in all cases, at some common, nominal saturation density ρ₀. On the other hand, analyses of correlations between the symmetry energy, its density slope, and the neutron skin thickness are typically done utilizing families of phenomenological models, such as large sets of Skyrme interactions or relativistic mean-field (RMF) models [131]. These models are constructed so as to have in common good saturation properties (usually by adjusting parameters to empirical properties of nuclei) while differing in the slope of the symmetry energy which, at saturation, is essentially a measure for the pressure in pure neutron matter (see Equation 17). Already almost two decades ago, Brown [132] considered a set of Skyrme interactions whose predictions of the density slope of the NM EoS around normal density differed dramatically and found a linear relation between such derivative and the neutron skin thickness in ²⁰⁸Pb. Similar investigations have been and continue to be done with RMF models, with families of interactions constructed so as to span a large range of L values. For instance, RMF models such as NL3 [133] and IU-FSU [134] give values of L equal to 118.2 and 47.2 MeV, respectively. (Not surprisingly, these models span a large range of both neutron skin values and stellar radii). In Roca-Maza et al. [24], the authors utilize a large set of RMF models all of which describe accurately the nuclear binding energies and charge radii across the periodic table (which should constrain tightly the binding energy and saturation density of SNM). On the other hand, the same models predict very different neutron root-mean-squared (r.m.s.) radii, since the isovector channel is poorly constrained [24].

Our EoS are microscopic and parameter-free and we are not in the practice of constructing families of parameterized EoS models to establish phenomenological correlations. Nevertheless, for the purpose of demonstration, next we wish to perform a study meant to highlight the role of neutron matter pressure for the neutron skin thickness once the uncertainty associated with the saturation point in SNM, cf. Figure 2, is removed. To that end, we will construct “semi-microscopic” models of asymmetric matter as follows: for the symmetric part, we will use an established phenomenological EoS, such as the one from Alam et al. [135]. For the neutron matter part, currently our focal point, we will continue to use the chiral EoS presented in section 2.3.2. We then proceed treating these six cases (three chiral orders and two cutoffs) as six EoS models differing in their NM components. Figure 7 shows the phenomenological EoS of SNM in comparison with our chiral predictions with cutoff equal to 450 MeV. Figure 8 displays the symmetry energy obtained with our microscopic NM EoS combined with SNM EoS represented by the black curve in Figure 7. We also include in the figures the results of several analyses and constraints [136–139]. The predictions based on our microscopic NM EoS are considerably softer than those constraints above normal density. Table 4 contains values for the parameters defined previosly through Equation (15), for the six EoS models which we have constructed as described.

We now proceed to discuss the spread of our values for the symmetry energy, the L parameter, and the isovector incompressibility in the framework of chiral uncertainties of the NM EoS. We recall that one of the strengths of χEFT is the opportunity of order-by-order improvement of the predictions. Naturally, the truncation error at a given order should be a reasonable measure of the uncertainty which arises from omitting the next order contributions. If the value of the observable X has been calculated at order n+1, than the truncation error at order n can be estimated as the difference between the values at order n+1 and n:

which is a reasonable way to estimate what one is missing by retaining only terms up to order n. On the other hand, if X_n+1 is not known, then some alternative prescription must be used. We use the definition [87]

where Q is a momentum characteristic for the observable under consideration and Λ is the cutoff. For the fifth (and highest) order, we use Equation (20) and we find it reasonable to define Q as the r.m.s. value of the relative momentum of two neutrons in neutron matter at the given density [see [123] and references therein].

We wish to express our final results for the symmetry energy, the slope parameter, and the isovector incompressibility at N³LO. To that end, we average the predictions for the quantity X obtained with the two values of the cutoff separately at N³LO and N⁴LO, yielding X̄4 and X̄5, respectively. The truncation uncertainty at N³LO can then be estimated as ΔX=|X̄4-X̄5|. As an alternative, we choose to take the largest of the errors at the two cutoff values.

Applying that prescription, we obtain for the symmetry energy, the slope parameter, and the isovector incompressibility at N³LO (all numbers in MeV):

We see that K_sym shows large variations, which reflect the extreme sensitivity of the second derivative to the details of the interactions for each of the curves in Figure 8. We emphasize that variations among those curves are due entirely to the NM predictions.

A phenomenological study of the EoS based on Skyrme density functionals [135] reports the slope parameter to be L = 65.4±13.5 MeV, whereas the isovector incompressibility is found to be within the range K_sym = −22.9±73.2 MeV. Lattimer and Lim [140] determined L to be between 40.5 and 61.9 MeV. For the isovector incompressibility, they suggest a linear relation between K_sym and L, that is, K_sym ≈ aL − b, with a, b equal to 3.33 and 281 MeV, respecively [128], or 2.867 and 260 MeV [127]. More recent constraints obtained from tidal deformabilitie inferred from GW170817, report 30 < L < 86 MeV and −140 < K_sym < 16 MeV or 40 < L < 62 MeV and −112 < K_sym < −52 MeV [141].

Before closing this section, we take the opportunity to address the pressure in neutron matter at saturation density, which, for the EoS of SNM which we have chosen is equal to 0.155 fm⁻³. Using Equation (17) and the uncertainty on L, we find (in MeV/fm³):

The neutron skin is defined as the difference between the r.m.s. radii of the neutron and proton density distributions:

where

i = n, p and T_n, T_p = N, Z respectively.

As mentioned before, the neutron skin thickness, particularly for ²⁰⁸Pb, is of great contemporary interest due to the possibility of constraining the slope of the symmetry energy through skin measurements [24, 142–147].

It is remarkable that the relation between the mass and the radius of neutron stars is uniquely determined by the EoS together with the star's self-gravity through the Tolman-Oppenheimer-Volkoff (TOV) equations of General Relativity [152]. In fact, although the detailed structure of a neutron star is complex and varies as a function of density, the part of its core mostly composed of a uniform liquid of neutrons, protons, and leptons in β-equilibrium accounts for almost all the mass and the volume. Therefore, these compact systems are intriguing testing grounds for both nuclear physics [153–156] and General Relativity. Extensive effort has been and continue to be devoted to constraining properties of compact stars from astrophysical observations see, for instance [156–160].

The largest possible value for the mass of a neutron star was estimated by Rhoades and Ruffini [161] based on the following assumptions: (1) General Relativity is the appropriate theory of gravitation; (2) the EoS obeys the Le Chatelier's principle (∂P/∂ϵ ≥ 0) and is consistent with causality, ∂P/∂ϵ ≤ c²; and (3) the EoS is reliably known below some density. Subject to these conditions, it was determined that the maximum mass of a neutron star cannot exceed 3.2 solar masses. Note that, releasing the causality constraint, the limit can be as high as 5 solar masses [162, 163] due to the increased stiffness of the EoS at supernuclear densities.

While the maximum mass is mostly determined by the stiffness of the EoS at densities greater than a few times saturation density, the star radius is impacted mainly by the slope of the symmetry energy. More precisely, it is closely connected to the internal pressure (that is, the energy gradient) of matter at densities between about 1.5ρ₀ and 2-3ρ₀ [157]. The mass and the radius of the neutron star are predicted by the TOV equation as we review next.

The equations for a perfect fluid in hydrostatic equilibrium allow to determine the pressure and the total mass-energy density as a function of the radial distance from the center of the star. These coupled equations are

with

where ϵ is the total mass-energy density. The star gravitational mass is

with R the value of r where the pressure vanishes. It's worth recalling that no mass limit exists in Newtonian gravitation.

Recently, the LIGO/Virgo [164] detection of gravitational waves originating from two neutron stars spiraling inward and merging, the neutron star merger GW170817, has generated even more interest and excitement around these highly exotic systems.

The dimensionless tidal deformability is related to the neutron star response to the tidal field induced by the companion star and is defined as

where the Love number k₂ reflects the quadrupole component of the gravitational potential induced by the companion star at the surface [165]. It depends on the neutron star compactness, M/R, and the energy density and pressure profile of the star. The tidal deformability can be obtained by solving the appropriate equations together with the TOV equations which yield the M(R) relation [166]. Hence, the merger detection can provide constraints on the star radius based on the tidal deformabilities of the colliding system [167]. In fact, the August 2017 first direct detection of a binary neutron star merger helped establish new limits on the radius of a 1.4 M_⊙ neutron star. Additional references addressing the radius of a 1.4 M_⊙ neutron star include [166, 168–172].

The correlation between the neutron skin thickness (discussed in section 3.3) and the radius of a neutron star originates from the sensitivity of the star radius to the pressure at normal density. Note that such correlation weakens as the mass increases see, for instance [110], which is why the radii of lighter stars are good candidate to help constrain the neutron skin of ²⁰⁸Pb and, in turn, the slope of the symmetry energy around saturation density. Based on these considerations, an upper limit of 0.25 fm was found for the neutron skin thickness of ²⁰⁸Pb. Additional observations from the LIGO-Virgo collaboration scheduled for 2019 are likely to provide stronger constraints.

In the remainder of this section, after reviewing how the EoS of β-stable matter is obtained from conditions of charge neutrality and energy minimization (section 4.2), we will address (spherical) neutron star properties, with emphasis on the radius of a “typical” neutron star, namely one with a mass approximately equal to 1.4 M_⊙. The reasons for this choice have been given in the previous paragraph.

In this section, we review the basic equations which we use to obtain the EoS for stellar matter in β-equilibrium.

The total energy per particle, e_tot, related to the total energy density, ϵ_tot, by e_tot = ϵ_tot/ρ, for neutrons and protons in β equilibrium with leptons (electrons and muons) is given by:

where Y_i, i = n(p), stands for the neutron(proton) fraction. On the right-hand side are the baryon contributions including their rest energies (first three terms), and the relativistic electron and the muon energies per baryon (last two terms). Note that, in the equation above, e₀ and e_sym are the EoS of symmetric nuclear matter and the symmetry energy, respectively. All terms are functions of density.

The relativistic energy density for particle species “i” having Fermi momentum (k_F)_i is given by

where γ is an appropriate degeneracy factor. The partial densities are related to the respective Fermi momenta as

which gives, for spin-12 fermions (γ = 2),

with the corresponding particle fractions given by

The chemical potential for species “i” is defined as

where we have used Equations (38–39) to perform the derivatives with respect to the upper integration limit.

The standard procedure is to minimize the total energy per particle with the constraints of fixed baryon density and global charge neutrality:

and

The resulting set of equations allow to solve for the various lepton fractions from which one can easily obtain the corresponding energy densities.

For the purpose of applying the Lagrange multipliers method, we define the functional

where

and

and set

Equations (45–48) then yield

Thus,

and

The equations above allow to solve for the various lepton fractions and, through Equation (38), the corresponding energies are easily obtained. For electrons, we find the ultra-relativistic approximation to be appropriate and set their rest energy to zero in Equation (38).

In Figure 10, we show the predicted fractions for the various species (neutrons, protons, and leptons) at the three highest orders of χEFT which we consider. We note that the proton fraction goes up to just above 10% at the highest densities being considered. This is a rather low value, most likely related to the relatively soft nature of the symmetry energy displayed in Figure 8. It implies that neutron stars (with central densities up to those included in the figure) will not cool down via direct Urca processes.

One is now in the position to calculate the pressure in β-stable matter. The pressure is related to the energy density through

In order to continue the discussion started in section 3.2 and extend it to neutron star radii in a consistent manner, in this section we will use the same interactions constructed in section 3 in terms of an empirical SNM EoS and the microscopic NM chiral EoS. Note, from section 4.2, that the symmetry energy, and thus the EoS of both SNM and NM are needed to obtain the various particle fractions. As mentioned earlier, those fractions tend to be rather small (see Figure 10), and thus the results shown in this and the next sections are to a large extent determined by our chiral predictions in NM.

We close this section by showing in Figure 11 the calculated pressure in β-stable matter at the third, fourth, and fifth orders of the 2NF together with the 3NF constructed as described in section 2.2.2. Of course we are referring to the chiral 3NF appropriate for NM, where the LECs c₄, c_D, and c_E vanish [56].

The BHF approach to nuclear matter is appropriate for the description of homogeneous matter, such as a homogeneous fluid of nucleons. Below nuclear densities, the chiral EoS are joined with the crustal equations of state from Harrison and Wheeler [173] and Negele and Vautherin [174], performing a smooth interpolation between the two EoS. The crust has crystal-like composition, and contains light [173] or heavy [174] metals together with a gas of electrons.

We now proceed to discuss specifically neutron star predictions. As mentioned above, we will focus on the radius of a star with mass equal to 1.4 M_⊙. We note, in passing, that the increased population of neutron stars observed around the mass range of 1.4 M_⊙ may be related to the physics of white dwarfs, atomic stars supported by electron (rather than neutron) degeneracy pressure. Since the Chandrasekhar limit of white dwarfs is approximately 1.4 M_⊙ [175], their collapse is likely to generate neutron stars in that mass range.

As stated in section 1, χEFT is a low-energy theory and thus limitations to its domain of applicability must be carefully considered. To begin with, the chiral symmetry breaking scale, Λ_χ ≈ 1 GeV, imposes clear limitations on the momentum or energy ranges where pions and nucleons can be taken as suitable degrees of freedom [15, 17]. Furthermore, the cutoff parameter Λ appearing in the regulator function (cf. Equation 10), has the purpose to remove high momentum components. Naturally, the strength of the cutoff determines to which degree such high-momentum components are suppressed. On the other hand, central densities of compact astrophysical systems can reach as high as several times the density of normal saturated matter, resulting, of course, in the presence of Fermi momenta which are beyond the reach of χEFT. Therefore, if one wishes to make predictions based, to some extent, on χEFT, methods to extend those predictions must be devised.

It has been observed that the pressure as a function of baryon density (or mass density) for a very large number of existing EoS can be fitted by piecewise polytropes, namely functions of the form P(ρ) = αρ^Γ [176]. (Note that, in our definitions, ρ denotes the baryon density). Guided by this observation, we find it reasonable to extrapolate the pressure predictions obtained from the EoS shown in Figure 11 using polytropes, as we have done in Sammarruca and Millerson [123]. More precisely, we employ our semi-microscopic predictions up to about 2ρ₀, where ρ₀ is defined to be 0.16 fm⁻³, approximately the density of saturated matter. The reason for choosing 2ρ₀ as a matching density is as follows: since we are dealing with a perturbative expansion in the parameter Q/Λ, we base our arguments on the size of the expansion parameter for typical momenta of the system under consideration. The highest momentum for pure neutron matter around twice normal density is approximately 420 MeV, as obtained from the usual relation ρ=kF,n33π2, with k_{F, n} the neutron Fermi momentum. And of course, the highest momentum in β-stable matter is slightly lower due to the presence of a small proton fraction. In conclusion, we are still below (although getting close to) Λ ~ 450 − 500 MeV. Furthermore, the r.m.s. value of the relative momentum for two nucleons in infinite matter is lower than their maximum momentum, and in fact it can be estimated to be about 60% of the Fermi momentum [177]. Thus, on statistical grounds, we should be safe from “cutoff artifacts,” even in the presence of smooth regulators.

We then proceed to match polytropes with diverse adiabatic indices, preserving continuity of the pressure. The range of the polytropic index was taken to be between 1.5 and 4.5 [123] (based on guidelines from the literature [176]), and these extensions are calculated up to about 3ρ₀. At this density, every polytrope is again joined continuously with another set of polytropes spanning the same range in values of Γ. In this way, we are able to cover a large set of possible EoS continuations, simulating scenarios where the EoS displays different degrees of “softness” or “stiffness” in different density regions, and thus we can estimate a realistic uncertainty. We stress again that this procedure is a way to simulate the uncertainty arising from reasonable parameterizations of the EoS as determined by phenomenological studies in the literature, and is not to be understood as a replacement for a theoretical model. A demonstration of this procedure is shown in Figure 12 (top row) for the case of N³LO with Λ = 450 MeV. (Note that cgs units are adopted in Figure 12 as those are popular in astronomy and may facilitate comparison by other authors).

In Hebeler et al. [178], the β-equilibrated EoS based on microscopic chiral interactions in neutron matter, is extended to high densities employing a general piecewise polytropic extrapolation, which leads to a very large number of EoS. Applying causality and the requirement that the EoS must be able to support a mass of 1.97 M_⊙, the author select a range of possible EoS, ranging from “soft” to “stiff.” We show the resulting uncertainty band in the bottom row of Figure 12, noting that our uncertainty band from Figure 12 is consistent with it.

Having built up the EoS at all needed densities, we are now in the position to solve the TOV [152] star structure equations and obtain the mass as a function of the radius for a sequence of stars differing in their central densities, up to several times normal density. The M(R) relations we obtain are shown on the bottom row of Figure 13. Note that only combinations of polytropic indices which can support a maximum mass of at least 1.97 M_⊙ have been retained for the purpose of Figure 13, to account for the observation of a pulsar with a mass of 2.01 ± 0.04 M_⊙ [179]. It is appropriate to point out here that most recent observations [180] are consistent with the even higher value of 2.14 M_⊙. In future work, we will apply this new constraint, which will result in a more limited set of acceptable EoS.

The causality constraint imposes limitations and those are applied in Figure 13. That is, one must require that the speed of sound in stellar matter is less than the speed of light, a condition which can be expressed as dPdϵ< 1, where ϵ is the total energy density.

Table 6 shows the the radius and the central density of the 1.4 M_⊙ neutron star when the pressure curves at the fourth and fifth orders from Figure 11 are extrapolated via piecewise polytropes with adiabatic indices Γ₁ and Γ₂ as shown. The speed of sound at central density is also included. Confirming what we found in Sammarruca and Millerson [123], Table 6 demonstrates quite clearly that the radius is practically insensitive to how the continuation is done. In particular, no changes are observed due to variations of the polytrope attached at 3ρ₀, and changes by less that one kilometer occurr in response to varying the first polytropic index. Note that the central densities we predict for the canonical-mass star are typically in the order of, and can exceed 3ρ₀. These densities are at or above those where the spreading of the pressure can be quite large (see Figure 12). Evidently, the radius of a star with this kind of mass responds to pressures at much lower than central densities, in line with earlier observations (see, for instance [178, 181]), where the insensitivity of the radius to the higher densities was pointed out. Tables similar to Table 6 with Λ=500 MeV and at N⁴LO with changing value of Λ are not included but do lead to very similar observations.

At the same time, the very small spreading of the pressure at normal to moderately high densities (see Figure 11), would suggest similarity of the radius in all cases (differing in chiral order and/or cutoff). This is in fact the case. Taking into consideration both the truncation error and the uncertainty from the polytropes, one may state, estimating the error pessimistically, that RN3LO≈ 11.8 ± 1 km for M=1.4M_⊙.

For completeness, in Figure 13 the full M(R) relation is also displayed. However, we stress that, at the high central densities probed by the heaviest stars, it is not possible to make reliable statements at this time. Predictions are no longer constrained by the chiral theory and are mostly phenomenology.

To broaden the scopes of this discussion, we will include next a set of predictions obtained in a more “traditional” way rather than with χEFT. In particular, we will use a meson-exchange potential (the Bonn B potential [86]) and the DBHF approach mentioned in the Introduction and used extensively in the past by one of the authors of this review [14]. We recall that the characteristic feature of the DBHF method is its ability to effectively take into account an important class of 3NF generated by the so-called “nucleon-antinucleon Z-diagrams” (see [14] and references therein).

However, one of the problems with the traditional approaches based either on meson-theoretic potentials (such as Bonn B or CD-Bonn [29]) or entirely phenomenological ones (such as AV18 [31]), is the absence of guidelines to select the 3NF contributions to be included (among the infinitely many possibilities). Typically, a particular diagram or set of 3NF diagrams are chosen to accompany the 2NF, but no well-defined link exists between the 2NF and the associated 3NF. On the other hand, the chiral approach, through the order-by-order scheme, prescribes exactly which 2NF, 3NF, and higher-body force must be retained at each order.

In Figure 14, we show the pressure in β-stable matter from the DBHF EoS. Comparison between the blue and the red curves demonstrate that choosing a phenomenological SNM EoS (red curve) as compared to the microscopic one (blue curve) has only a minor impact on the EoS for β-stable matter (which is comprised mostly of neutrons), particularly from low to medium densities.

We then proceed to compare the pressure predictions based on the meson-exchange model with the predictions from Figure 11. As to be expected, differences become larger with increasing density, with the chiral EoS being substantially softer at the higher densities.

Once again, we place our focus on the radius of the average-mass star, which we find to be approximately equal to 12.5 km for the DBHF calculations. This value is reasonably close to our previous, chirally based predictions, which makes sense based on our earlier discussions and the fact that the DBHF pressure around normal density is not very different from the one in the chirally-based models (see Figure 15).

Finally, we employ the DBHF EoS in the TOV equations and calculate the M(R) relation, which is shown on the RHS of Figure 14. The blue and red curves correspond to the blue and red curves in Figure 14. Differences become noticible for the heavier stars and are consistent with those seen in Figure 14. In other words, the model with the larger pressure at the higher densities generates the larger maximum mass.

We end this exercise with an important comment: even if a theory (of nucleons and mesons) can be formally taken to high densities, as we have done with the DBHF predictions, the composition and thus the EoS of stellar matter in the inner core is simply unknown. At densities as high as those typical of compact stars, hyperons are expected to exist on simple energetic grounds. Similarly, other non-nucleonic degrees of freedom, such as quark degrees of freedom, can exist as the result of phase transitions. These possibilities have been explored by several groups (see, for instance [182–185]). Such investigations are not within the scope or the reach of χEFT. The polytropic extrapolations we have performed do indeed simulate a broad set of possible EoS consistent with current constraints but whose specific composition remains unknown.

Of course, we also take note of some other works aimed at incorporating aspects of chiral dynamics in the development of EoS suitable for astrophysical phenomena, such as Rapaj et al. [186]. In the latter reference, the authors calculate neutron star masses and radii with mean-field models whose parameters are made consistent with a chiral EoS at low to moderate densities. Constraints from χEFT on neutron star tidal deformabilities were investigated in Lim and Holt [166]. General relativistic simulations of neutron star mergers based on the EoS of Bombaci and Logoteta [111] have been reported in Endrizzi et al. [187].

We close this section with a few final remarks addressing predictions vs. constraints. Masses of neutron stars can be and have been measured with high precision. However, simultaneous measurements of radii are much more problematic. Some techniques do exist, such as those based on photospheric radius expansion [188]. Current observations have begun to determine the M(R) relation. In Steiner et al. [159], the authors determine the radius of a 1.4 M_⊙ neutron star to be between 10.4 and 12.9 km. Furthermore, from their Bayesian analysis of several EoS parameterized so as to be consistent with a baseline data set (see [159] and references therein), they are able to determine the M(R) relation within a range of masses. Our predictions fall within those constraints, shown in Figure 13 as the shaded purple area. Recent LIGO/Virgo measurements have constrained the radius of a 1.4 M_⊙ neutron star to be between 11.1 and 13.4 km [164, 167]. The predictions from our group are well within these new constraints.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.