Chiral Effective Field Theory and the High-Density Nuclear Equation of State

Recent advances in neutron star observations have the potential to constrain the properties of strongly interacting matter at extreme densities and temperatures that are otherwise difficult to access through direct experimental investigation. At the same time, chiral effective field theory has developed into a powerful theoretical framework to study nuclear interactions and nuclear matter properties with quantified uncertainties in the regime of astrophysical interest for modeling neutron stars. In this article, we review recent developments in the chiral effective field theory approach to constructing microscopic nuclear forces and focus on many-body perturbation theory as a computationally efficient tool for calculating the structure, phases, and linear-response properties of hot and dense nuclear matter. We also demonstrate how effective field theory combined with Bayesian methods enables statistically meaningful comparisons between nuclear theory predictions, nuclear experiments, and observational constraints on the nuclear equation of state. 1 ar X iv :2 10 1. 01 70 9v 1 [ nu cl -t h] 5 J an 2 02 1

Recent advances in neutron star observations have the potential to constrain the properties of strongly interacting matter at extreme densities and temperatures that are otherwise difficult to access through direct experimental investigation. At the same time, chiral effective field theory has developed into a powerful theoretical framework to study nuclear interactions and nuclear matter properties with quantified uncertainties in the regime of astrophysical interest for modeling neutron stars. In this article, we review recent developments in the chiral effective field theory approach to constructing microscopic nuclear forces and focus on many-body perturbation theory as a computationally efficient tool for calculating the structure, phases, and linear-response properties of hot and dense nuclear matter. We also demonstrate how effective field theory combined with Bayesian methods enables statistically meaningful comparisons between nuclear theory predictions, nuclear experiments, and observational constraints on the nuclear equation of state.

Introduction
Neutron stars are one of Nature's most intriguing astronomical objects and provide a unique laboratory for studying strongly interacting, neutron-rich matter under extreme conditions. With masses ≈ 1 − 2 times that of the Sun and radii of only ≈ 10 km, neutron stars contain the densest form of matter in the observable Universe and lie just at the threshold for collapse to a black hole. Much has already been learned about neutron stars through mass and radius measurements, pulsar timing, x-ray observations, and gravitational-wave measurements of binary mergers in the new era of multimessenger astronomy (see, e.g., Refs. (1, 2, 3) for reviews). But many interesting questions remain to be answered, especially regarding the nature of ultra-compressed matter located in the inner cores of heavy neutron stars where a variety of exotic new states of matter have been theorized to exist.
The mass of the Sun is M ≈ 2 × 10 30 kg.
While neutron stars are bound together by gravity acting over macroscopic length scales, strong short-ranged nuclear interactions provide the essential pressure support to counteract gravitational collapse. The central densities in the heaviest neutron stars may reach up to 5−10 n0, where n0 ≈ 0.16 fm −3 is the nucleon number density typical of heavy atomic nuclei (the associated mass density is ρ0 ≈ 2.7 × 10 14 g cm −3 ). Although the strong interaction is in principle described by quantum chromodynamics (QCD) over all relevant energy scales, at present no systematic computational method is available to calculate the properties of the high-density matter in the inner cores of heavy neutron stars.

QCD: quantum chromodynamics
ChEFT: chiral effective field theory With chiral effective field theory (ChEFT), however, a powerful tool has emerged to carry out microscopic calculations of nuclear matter properties at densities up to around 2n0. Instead of QCD's quarks and gluons, ChEFT is formulated in terms of nucleons and pions (and delta isobars), which are the effective strong interaction degrees of freedom present throughout most of the neutron star interior. In its range of validity, ChEFT provides a systematic expansion for two-and multi-nucleon interactions consistent with the symmetries of low-energy QCD. The unresolved short-distance physics is parameterized in terms of contact interactions whose low-energy couplings are fitted to experimental data. An essential advantage over phenomenological approaches is that theoretical uncertainties can be quantified by analyzing the order-by-order convergence of the ChEFT expansion. In the last few years, the combination of systematic nuclear matter predictions from ChEFT, uncertainty quantification, and neutron star observations has developed into a new avenue for constraining the high-density regime of the nuclear equation of state (EOS).

EOS: equation of state
MBPT: many-body perturbation theory In this review our aim is to describe recent advances in microscopic ChEFT calculations of the nuclear EOS and their application to neutron stars. We highlight many-body perturbation theory (MBPT) as an efficient framework for nuclear matter calculations at zero and finite temperature based on chiral two-and multi-nucleon interactions. We also discuss Bayesian methods for quantifying and propagating statistically meaningful theoretical uncertainties. Together with nuclear experiments, astrophysical simulations, and neutron star observations, next-generation ChEFT calculations will be crucial to infer the nature of the extreme matter hidden deep beneath the surface of neutron stars.
The review is organized as follows. In Section 2 we focus on recent progress in deriving nuclear forces from ChEFT and renormalization group (RG) methods to improve the manybody convergence in nuclear matter calculations. We then dedicate Section 3 to recent high-order MBPT calculations of the moderate-density nuclear EOS at zero temperature and advances in the Bayesian quantification of EFT truncation errors. We also discuss finite-temperature calculations and nuclear thermodynamics. In Section 4 we review the present status of the high-density nuclear EOS constrained by nuclear theory, experiment, and observation in the era of multimessenger astronomy, emphasizing the importance of ChEFT. Section 5 ends the review with our summary and perspectives on future advances in nuclear matter calculations and their applications to astrophysics.

From microscopic interactions to the nuclear equation of state
In this section, we briefly review delta-less ChEFT and the construction of chiral nuclear interactions as microscopic input for many-body calculations. Applying RG methods allows one to systematically generate (perturbative) low-momentum interactions, for which the nuclear EOS and related observables can be efficiently calculated using MBPT. We discuss both zero-and finite-temperature MBPT, complementary many-body approaches, and the implementation of 3N interactions in nuclear matter calculations.

Chiral effective field theory for nuclear forces
The interactions among nucleons arise as an effective low-energy phenomenon of QCD, the theory of the strong interaction. At the momentum scales relevant for nuclear physics, p ∼ mπ, QCD is strongly coupled and features nonperturbative effects such as spontaneous chiral symmetry breaking and the confinement of quarks and gluons into hadrons. Direct applications of QCD to hadronic physics at finite density, where lattice QCD faces a formidable sign problem, are therefore extremely challenging and not feasible at present and in the near future. However, one can construct a systematic description of nuclear physics in terms of the effective degrees of freedom at low energies: nucleons and pions (and delta isobars). This effective description is given by ChEFT (8,9,10,11,12).
The starting point of ChEFT is to write down the most general Lagrangian consistent with the symmetries of low-energy QCD, in particular, the spontaneously broken chiral symmetry, for which pions are the (pseudo) Nambu-Goldstone bosons. This naturally sets a limit for the applicability of ChEFT, i.e., the breakdown scale Λ b will be of order of the chiral symmetry breaking scale Λχ ∼ 1 GeV. A truncation scheme, known as power counting, is then needed to organize the infinite number of operators in the effective Lagrangian in a systematic expansion. This expansion is governed by the separation of scales inherent in ChEFT, i.e., the power counting is according to powers of a typical momentum p (or the pion mass) over the ChEFT breakdown scale, Q = max(p, mπ)/Λ b .
A Bayesian analysis of free-space NN scattering with several (not-too-soft) chiral NN potentials in Weinberg power counting estimated Λ b ≈ 600 MeV (13).
In perturbative EFT, both power counting and ultraviolet renormalization are essentially unambiguous and straightforward. The situation is different for applications of ChEFT in nuclear physics, where the calculational framework must be able to account for nonperturbative effects such as bound states (atomic nuclei) and large S-wave scattering lengths in NN scattering. While there has been some controversy in the literature as to how the ChEFT expansion should be set up precisely (see, e.g., Refs. (10,12) and references therein), the prevalent and most successful power counting for ChEFT (in particular regarding many-body applications) is the one first suggested by Weinberg.
NN: nucleon-nucleon 3N: three-nucleon 4N: four-nucleon Within Weinberg power counting, chiral nuclear interactions (and currents) are organized according to naive (i.e., perturbative) dimensional analysis. The nuclear potentials constructed at a given truncation order in the ChEFT expansion are then used for computing observables. Renormalization in this approach is approximative, and carried out by equipping the potentials with regulator functions that suppress contributions above a cutoff scale Λ Λ b , typically chosen in the range 450 − 600 MeV. That is, the cutoff independence of the observables will be achieved only approximatively through Λ-dependent low-energy constants (LECs), which have to be fit to experimental data at a given scale. The residual cutoff dependence can then be attributed to higher-order terms in the expansion, so results are expected to become less cutoff dependent with increasing truncation order.
LECs: low-energy constants are in practice optimized for a given value of Λ to reproduce low-energy NN scattering data and few-nucleon observables, see, e.g., Refs. (12,14). Figure 1 depicts the hierarchy of nuclear interactions up to fifth order (or N 4 LO) in the chiral expansion without delta isobars. At each order the interactions are composed of short-range contact interactions as well as one-and multi-pion exchanges at long-and intermediate distances, respectively. The LECs associated with pion exchanges have recently been determined with high precision through an analysis of pion-nucleon scattering within the framework of Roy-Steiner equations (15). The short-range LECs corresponding to NN couplings are generally fixed by matching to NN scattering data. Figure 1 shows that ChEFT naturally predicts the observed hierarchy of two-and multi-nucleon interactions, i.e., VNN > V3N > V4N, etc. The first nonvanishing 3N forces appear at N 2 LO in three topologies; from left to right: the long-range two-pion exchange (involving the pionnucleon LECs c1, c3, and c4), intermediate-range one-pion exchange-contact (∝ cD), and short-range 3N contact interaction (∝ cE). At N 3 LO the 3N forces are significantly more involved and operator-rich, and also 4N interactions start to contribute. Apart from the two N 2 LO 3N LECs, cD and cE, chiral interactions up to N 3 LO are completely determined by the πN and NN system. While N 4 LO NN forces have already been worked out, partly even at N 5 LO, the derivation of N 4 LO 3N interactions has not been finished yet. The 3N LECs cD and cE can be fit to (uncorrelated) few-body observables; for instance, the 3 H binding energy combined with, e.g., the charge radius of 4 He, the 3 H β-decay half-life, or the nucleon-deuteron scattering cross section. Also heavier nuclei and even saturation properties in infinite nuclear matter have been used to constrain 3N forces.
ChEFT with explicit delta isobars is currently less developed than the delta-less version we focus on here. For recent work on delta-full ChEFT, see, e.g., Refs. (16,17).
Although the residual regulator and cutoff dependence of obervables at a given chiral order is expected to decrease at higher orders, actual calculations show significant influence of these so-called regulator artifacts on the ChEFT convergence depending on the specific regularization scheme and computational framework. These issues have resulted in the development of a flurry of chiral potentials with nonlocal, local as well as semilocal regulators for a range of cutoff values; see, e.g., Table I of Ref. (18). Moreover, as discussed in Section 2.2, RG methods allow one to modify a given set of two-and multi-nucleon potentials such that observables are left invariant (up to RG truncations) but the convergence of manybody calculations is optimized. These RG transformations are most suitably formulated at the operator (i.e., Hamiltonian) level. The nuclear Hamiltonian constructed at a given order in the ChEFT expansion reads H = T kin + VNN(Λ, ci) + V3N(Λ, ci) + V4N(Λ, ci) + . . . , where Λ stands for the (initial) cutoff or resolution scale, and ci for the set of LECs inferred from fits to experimental data. The nuclear Hamiltonian is not an observable, and the basic idea of the RG is to exploit this feature to generate more perturbative Hamiltonians.

Perturbative chiral nuclear interactions
The strong short-range repulsion ("hard core") and tensor force found in nuclear potentials constructed at cutoff scales Λ 500 MeV question the applicability of perturbation theory for many-body calculations. In fact, nuclear many-body calculations were historically considered a nonperturbative problem (see also Section 2.3). Both features give rise to strong couplings between high-and low-momentum states, i.e., large off-diagonal matrix elements, which enhance the intermediate-state summations in perturbation theory. RG methods allow one to amend this feature while preserving nonperturbative few-body results. The initial application (19) of RG methods to study the scale dependence of nuclear forces was based on T -matrix equivalence, but in recent years the similarity renormalization www.annualreviews.org • Chiral EFT and the High-Density Nuclear EOS group (SRG) has been the standard RG method for "softening" nuclear interactions. The SRG decouples high-and low-momentum states through continuous infinitesimal unitary transformations, Hs = UsHU † s , described by a differential flow equation in the evolution parameter s. As the SRG flow progresses, the matrix-elements of the NN potential are driven toward a band-diagonal (or block-diagonal) form in momentum space, see Ref. (20) for illustrations. While NN observables are by construction invariant under any RG evolution of the NN potential, A-body observables will remain so only if one also consistently evolves the multi-nucleon part of the nuclear Hamiltonian. The SRG allows one to implement this in principle exactly, in terms of so-called "induced" many-body forces. However, for practical applications, a truncation of the consistent evolution of multi-nucleon interactions is required, e.g., at the 3N level. Neutron matter calculations with SRG-evolved chiral interactions truncated at the 3N level indicate that induced 4N forces are (in that case) negligible within uncertainties for a wide range of SRG resolution scales (21,22).

SRG: similarity renormalization group
There have been many developments recently in the application of ChEFT and SRG technology for the construction of high-precision nuclear potentials. Hebeler et al. (23) explored a set of low-momentum N 3 LO NN potentials combined with unevolved N 2 LO 3N forces where the two 3N LECs were fit to reproduce few-body data (assuming that the 3N contact interactions capture dominant contributions from induced 3N forces). For the softest of these potentials (with λ = 1.8 and Λ3N = 2.0 fm −1 ), which was found to predict nuclear saturation properties (24) and ground-state energies of light-to medium-mass nuclei in agreement with experiment (25), Stroberg et al. (26) have computed ground-state and separation energies of nearly 700 isotopes up to iron. Moreover, Hüther et al. (27) have constructed a family of SRG-evolved NN and 3N potentials up to N 3 LO. Also, Reinert et al. (14) have developed the first chiral NN potentials up to N 4 LO with semilocal regulators in momentum space. They showed that several N 3 LO contact terms present in previous generations of chiral NN potentials can be eliminated using unitary transformations, leading to considerably softer potentials (even without SRG evolution). The Weinberg eigenvalue analysis (28,18) is a powerful tool for quantifying and monitoring the perturbativeness of nuclear forces at different resolution scales. Given an NN potential, the Weinberg eigenvalues ην (W ) of the operator G0(W )VNN determine the (rate of) convergence of the Born series for NN scattering. Here, G0(W ) is the (free-space or in-medium) propagator as a function of the complex energy W . The Born series converges if and only if all eigenvalues satisfy |ην (W )| < 1. Bound states of the potential (such as the deuteron) correspond to ην (W ) = 1 at energies W < 0, so the free-space Born series diverges even for soft potentials. In nuclear matter at sufficiently high densities, however, Pauli blocking suppresses the (in-medium) eigenvalues associated with bound (or nearlybound) states. For potentials constructed at Λ 550 MeV, other sources of nonperturbative behavior (such as the hard core) are suppressed as well, both in free-space and in-medium (see, e.g., Ref. (20) for details). This implies that a nonperturbative treatment of in-medium NN scattering in the particle-particle channel (see Section 2.3) is not mandatory for these interactions. Instead, order-by-order MBPT calculations can be used for systematically studying the many-body convergence of (low-momentum) chiral nuclear interactions.
The free propagator is given by G 0 (W ) = (W − H 0 ) −1 , with H 0 the kinetic energy operator. The in-medium propagator involves also Pauli-blocking.

Many-body perturbation theory at zero and finite temperature
MBPT starts with partitioning the nuclear Hamiltonian H into a reference one-body part H0 = T kin + U and a perturbation H1 = V − U , where T kin is the kinetic-energy operator and U is an effective single-particle potential. We consider here NN-only potentials and discuss the implementation of 3N interactions in Section 2.5. The standard choice for U is the Hartree-Fock potential given by U (HF) i = j V ij,ij fj, with the antisymmetrized NN matrix elements V ij,ab = kikj|(1 − P12)VNN|kak b , the Pauli exchange operator P12, the momentum integral j = d 3 kj /(2π) 3 , and the zero-temperature distribution function fj = θ(kF − kj). (For simplicity, we assume here a single-species system and neglect spinisospin degrees of freedom.) In zero-temperature MBPT, the ground-state energy density E is obtained by expanding Hint about its reference value E0. Truncating the many-body expansion at a finite order L then leads to the approximation E(kF) E0(kF)+ L l=1 E l (kF), where the Fermi momentum kF is in one-to-one correspondence with the particle number density via n(kF) = i fi(kF).
The MBPT series is in fact a divergent asymptotic series, but the divergent behavior is expected to appear only for high truncation orders L 20 (29).
The first-order correction is determined by the expectation value of U (HF) i (30). At higher orders it is useful to represent the contributions diagrammatically, e.g., as Hugenholtz diagrams. The diagram and expression for the second-order contribution E2 are given by with the distribution functions fij = fifj (holes) andf ab = (1−fa)(1−f b ) (particles), energy denominator D ab,ij = εa + ε b − εi − εj, and single-particle energies εi = k 2 i /(2M ) + U (HF) i . Writing down the expression associated with a diagram follows these simple rules: In Hartree-Fock MBPT, the −U part of H 1 = V − U cancels all diagrams involving single-vertex loops (30,24,22).
• each vertex gives a factor V ij,ab , with i and j (a and b) being the lines directed towards (away from) the vertex, • downwards lines give factors of fi while upwards lines give (1 − fi), corresponding to hole and particle excitations of the reference ground-state, respectively, and • for adjacent vertices there is an energy denominator given by subtracting the energy of the reference ground-state from the excited state corresponding to the particle and hole lines that are crossed by a virtual horizontal line between the two vertices.
Each diagram's overall factor can be inferred from the diagrammatic structure as well (30). For instance, the expression of the third-order particle-particle (pp) diagram reads Finding all valid diagrams (and associated expressions) at a given MBPT order has been formalized using graph-theory methods (31). Together with automated code generation for the efficient Monte Carlo integration of arbitrary MBPT diagrams developed in Ref. (24), a fully automated approach to MBPT calculations has become available.
In the traditional Brueckner (or G-matrix) approach (32), the pp ladder diagrams are resummed to all orders, motivated by the large high-momentum components of traditional NN potentials to which the pp ladders are particularly sensitive. The pp bubbles in these diagrams are even ultraviolet divergent if the potential is not sufficiently suppressed at high momenta. For modern low-momentum potentials, however, the pp ladders no longer play a distinguished role in the many-body expansion, and explicit MBPT calculations at third and fourth order have shown that they are not enhanced compared to other diagrams at the same order (24). Nevertheless, partial diagrammatic resummations are still pertinent for consistent calculations of in-medium single-particle properties and response functions as performed in the self-consistent Green's functions method (for more details see Section 2.4).
The consistent generalization of MBPT to finite temperatures (T > 0) is a nontrivial issue. From the standard finite-T perturbation series for the grand-canonical potential the free energy density F(T, µ) is obtained via the thermodynamic relation F(T, µ) = Ω(T, µ) + µ n(T, µ). Here, the density is given by n(T, µ) = −∂Ω(T, µ)/∂µ. The issue is now that the relations between (F, n) and (E, n) obtained in finite-and zero-T MBPT, respectively, do not match in the limit T → 0 (as discussed further below).
The use of the grand-canonical ensemble is required for the evaluation of quantum-statistical averages in the thermodynamic limit.
Regarding this, we first consider the finite-T expression for the second-order diagram, Equation (4) differs only slightly from E2(kF) in Eq. (1). First, the energy denominator is replaced by G2 = (1 − e −D ab,ij /T )/(2D ab,ij ). The numerator in this expression vanishes at any zero of the denominator, i.e., there are no poles at finite T . In the T → 0 limit the integration regions corresponding to the two terms in the numerator of G2 separate into two equivalent parts (with integrable poles at the integral boundaries), i.e., G2 → 1/D ab,ij for T → 0. These features pertain for higher-order diagrams (33). The second difference at finite T compared to T = 0 is that the fi = fi(T, µ) are Fermi-Dirac distributions instead of step functions centered at the Fermi energy ε k F .
The poles (at the integral boundary) at T = 0 lead to nonanalyticities in the asymmetry dependence of the nuclear EOS, see Section 3.2.
Similar to the free Fermi gas (i.e., MBPT with U = 0 and L = 0), for Hartree-Fock MBPT at L = 1, the chemical potential µ at T = 0 matches the reference Fermi energy ε k F , with n(T, µ) = i fi(T, µ) (33). But these relations cease to be valid at higher orders due to higher-order contributions in the expression for n(T, µ). Note that these contributions involve factors ∂fi/∂µ = fi(1 − fi)/T , which become δ(εi − µ) at T = 0 (so there is a nonvanishing contribution at T = 0). Contributions involving factors fi(1 − fi)/T are also present in certain perturbative contributions to Ω, starting at fourth order for Hartree-Fock MBPT (33). [For U = 0, they appear already at second order.] These contributions can be associated with the presence of additional so-called anomalous diagrams in finite-T MBPT, see Refs. (34, 33) for more details. As evident from the discussion above (i.e., below Eq. (4)), the T → 0 limit of the finite-T expressions for normal (i.e., not anomalous) contributions Ω l matches the corresponding zero-T contributions E l , except that the reference Fermi energy is replaced by the (true) chemical potential. Therefore, a consistent finite-T version of Hartree-Fock MBPT for L 3 would be given by where F l = Ω l (for l = 1, 2, 3) and the auxiliary "chemical potential"μ is related to the density via n(T,μ) = i fi(T,μ), implyingμ → ε k F in the T → 0 limit.
The true chemical potential is obtained from F (T,μ) via µ = ∂F /∂n.
In the U = 0 case, the method for constructing a finite-T perturbation series of the form of Eq. (5) for any L is well known (34): one expands each contribution to F(T, µ) aboutμ according to µ =μ + L l=1 µ l (T,μ) while neglecting all terms beyond the truncation order L. This process can also be applied to Hartree-Fock MBPT (22), with the caveat that the single-particle potential has to be evaluated atμ, i.e., no derivatives of U (HF) i (T,μ) inμ appear. In both cases, U = 0 and Hartree-Fock MBPT, the resulting perturbation series for the free energy reproduces zero-T MBPT at each truncation order L, even though the terms F l contain anomalous contributions for l 4 (l 2, for U = 0). The fact that Eq. (5) results from a truncated re-expansion shows explicitly that the original grand-canonical series is not consistent with zero-T MBPT. For general arguments why the free-energy series is expected to give improved results compared to grand-canonical MBPT, see Refs. (33,35).
The terms µ l are determined by the requirement that the truncated expansion of n(T, µ) aboutμ reproduces n(T,μ) = Altogether, MBPT as formulated in the free-energy series [Eq. (5)] provides a consistent framework for nuclear matter calculations at zero-and finite-temperature, where many-body uncertainties can be systematically assessed by increasing the truncation order L. Although the number of MBPT diagrams increases rapidly with L, the technologies recently developed for automated diagram generation and evaluation (31,24) enable calculations at highenough orders to probe in detail the many-body convergence for chiral low-momentum NN and 3N interactions. Furthermore, exploring MBPT with single-particle potentials beyond the Hartree-Fock level is an important task for future research. In particular, the singleparticle potential U can be chosen at each truncation order such that the grand-canonical and free-energy series are also equivalent for L > 1 (33); for Hartree-Fock MBPT, they are only equivalent for L = 1. First investigations of this order-by-order renormalization of the single-particle potential have shown that higher-order contributions to U can have a significant effect on low-order MBPT results and the many-body convergence (36).

Other many-body methods
The advances in ChEFT and RG methods have established MBPT as a central approach for studying the nuclear EOS at zero and finite temperature. While MBPT is the focus of this review, various other many-body methods have been applied in initial nuclear matter studies with chiral NN and 3N interactions. In particular, nonperturbative frameworks are important to benchmark the MBPT convergence and probe aspects of many-body physics beyond the nuclear EOS. Below we will briefly discuss the self-consistent Green's functions (SCGF) approach and Quantum Monte Carlo (QMC) methods. Other methods not discussed here for brevity are coupled-cluster (CC) theory (37), the in-medium SRG (4), and lattice EFT (38). Systematic comparisons between different many-body frameworks will provide a coherent picture of microscopic interactions and nuclear many-body properties.
The SCGF approach (39,40) is based on the self-consistent computation of in-medium propagators (or Green's functions) in Fourier (Matsubara) space, corresponding to the toall-orders resummation of some perturbative contributions to the propagators. SCGF calculations of the nuclear EOS at zero and finite temperature (41) have been implemented using the in-medium T -matrix approximation, where the ladder diagrams are resummed to all orders, providing a thermodynamically consistent generalization of Brueckner theory (40). Furthermore, SCGF calculations have been used to benchmark the order-by-order convergence of MBPT (up to third order) in neutron matter (42). The energy per particle obtained in SCGF and MBPT was found to agree well for a range of unevolved chiral NN and 3N interactions up to N 3 LO. The SCGF approach allows for fully consistent computations of response functions and transport properties, which will be vital for comparisons with MBPT calculations of these quantities. QMC refers to a family of stochastic methods that solve the many-body Schrödinger equation through random sampling (6). As such, QMC methods are truly nonperturbative and provide important benchmarks for many-body methods with basis expansions. However, apart from the fermion sign problem a caveat is that most QMC methods require local nuclear potentials to obtain low-variance results, restricting both the regularization scheme and the interaction operators that can be included in the ChEFT expansion. QMC calculations with local chiral NN and 3N potentials up to N 2 LO have been carried out in neutron matter (43,44) and recently also symmetric nuclear matter (45). The regulator artifacts are (due to Fierz-invariance breaking) significantly larger compared to MBPT calculations with nonlocal potentials. On the other hand, since QMC methods are not restricted to soft interactions, a much wider range of momentum cutoffs can be studied with QMC. Hence, QMC methods can provide important insights into the residual cutoff dependence of observables and the breakdown scale of ChEFT at high densities.

Implementing three-body forces
Three-nucleon forces are crucial for understanding properties of finite nuclei and nuclear matter (46), such as drip lines along isotopic chains and nuclear saturation in SNM. Even though partial-wave decomposed matrix elements of chiral 3N forces have become available recently up to N 3 LO (47), implementing 3N forces in many-body calculations remains computationally difficult and usually requires approximations (48). The large uncertainties due to 3N forces, e.g., in the nuclear EOS at densities n n0, emphasize the need for improving these approximations as well as developing novel chiral NN and 3N potentials in general.
Normal ordering allows one to include dominant 3N contributions in many-body frameworks using density-dependent effective two-body potentials (49). Through Wick's theorem the general three-body Hamiltonian can be exactly normal ordered with respect to a finitedensity reference state (e.g., the Fermi sea of noninteracting nucleons or the Hartree-Fock ground state) instead of the free-space vacuum (20). This shifts contributions from the three-body Hamiltonian operator to effective zero-body, one-body, and two-body operators plus a residual (reduced) three-body operator. A many-body framework built for NN interactions can then incorporate a density-dependent effective interaction V med NN derived from V3N as VNN → VNN + ξ V med NN . The combinatorial factor ξ is determined by Wick's theorem and depends on the many-body calculation of interest. The matrix elements of V med NN are obtained by summing one particle over the occupied states in the reference state: with the shorthand notation |i = |kiσiτi , antisymmetrized 3N interactionsV3N, and momentum distribution function of the reference state f1.
In contrast to the (Galilean-invariant) NN potential, the effective two-body potential (6) depends on the center-of-mass momentum P of the two remaining particles. Hence, both potentials cannot be straightforwardly combined in a partial-wave basis and different approximations for the P dependence have been used to enable applications to nuclear matter. Under the assumption that P = 0 first implementations evaluated Eq. (6) semi-analytically in symmetric nuclear matter and pure neutron matter starting from the N 2 LO 3N interactions (50,51,23). Extensions to asymmetric nuclear matter and finite temperature have followed (35,52,42,53), and a new method that allows for the construction of an effective two-body potential from any partial-wave decomposed 3N interaction in an improved P angle-averaging approximation has been developed (53). The latter approach is especially advantageous for studying 3N forces at N 3 LO (53), bare and SRG-evolved, and in different regularization schemes. Semi-analytic expressions along the lines of Ref. (50) have been derived up to N 3 LO and also partially to N 4 LO (54).
The three-body term in the normal-ordered Hamiltonian cannot be implemented using effective two-body potentials. In nuclear matter such residual 3N contributions have been studied in CC (55) and MBPT calculations (56,57,58,24). Explicit calculations of the residual 3N diagram in MBPT at second order (see the margin note) showed for a range of chiral interactions that its contribution is typically much smaller than both the overall EFT truncation error and the individual contributions from the other MBPT diagrams up to this order (24). While these findings give some justification for the commonly used approximation where residual 3N contributions are neglected, the automated approach introduced in Ref. (24) implements chiral NN, 3N, and 4N interactions exactly in nuclear matter calculations using a single-particle spin-isospin basis. Combined with high-performance computing, this method sets the stage for systematic studies of ChEFT interactions in MBPT up to high orders and without the mentioned approximations.
Residual 3N diagram at second order:

Nuclear equation of state at zero and finite temperature
In this Section we survey recent nuclear matter calculations up to n ≈ 2n0 in MBPT with chiral NN and 3N interactions. We discuss advances in the quantification and propagation of EFT truncation errors, confront different microscopic constraints on the nuclear symmetry energy with experiment, and examine contributions beyond the standard quadratic expansion of the EOS in the isospin asymmetry. We conclude the Section with results for the nuclear liquid-gas phase transition at finite temperature. interactions (23). Several general observations can be gleaned. Nuclear interactions are much stronger in SNM compared to PNM, which is closer to the free Fermi gas (FFG, solid lines). Consequently, the uncertainties are larger in SNM, especially for densities n n0. In PNM they are well controlled for n n0, and a wide range of chiral NN and 3N interactions leads to similar results for PNM (see, e.g., Refs. (60,46,61)). Increasing uncertainties toward higher densities are predominantly due to 3N interactions. Although the complexity of 3N interactions is much reduced in PNM (51), they provide at all values of δ important repulsive contributions that grow stronger with the density than those of NN interactions. The 3N interactions are therefore crucial for understanding the high-density EOS and the structure of neutron stars. In PNM all chiral interactions up to N 3 LO are completely determined by the πN and NN system. The intermediate-and short-range 3N interactions at N 2 LO that are proportional to the LECs cD and cE, respectively, vanish (for regulators symmetric in the particle labels) due to the coupling of pions to spin and the Pauli principle, respectively. Also the long-range two-pion exchange 3N forces at N 2 LO are simplified since the LEC c4 does not contribute. This allows for tight low-density constraints on the neutron-rich matter EOS from PNM calculations and systematic highdensity extrapolations (see Section 4).

FFG: free Fermi gas
Nuclear matter represents an ideal system for testing nuclear interactions at the densities accessible to laboratory experiments and their implementation in many-body methods. As illustrated in Fig. 2 (left), the nuclear EOS in the vicinity of n0 is (to good approximation) characterized by only a few experimentally accessible quantities. That is, the EOS of SNM can be expanded about its minimum n0 as E(n, δ = 0)/A ≈ E0 + (K/2) η 2 , with the saturation energy E0 = E(n0, 0)/A, incompressibility K, and η = (n − n0)/(3n0). Further, explicit ANM calculations with chiral NN and 3N interactions (53,62,63) have shown that the asymmetry dependence of the nuclear EOS is reasonably well reproduced by the standard quadratic approximation E(n, δ)/A = E(n, 0)/A + Esym(n) δ 2 , where the symmetry energy expanded in density reads Esym(n) ≈ Sv + L η. In this approximation one finds E(n, 1)/A ≈ (E0 + Sv) + L η for PNM. Microscopic predictions and empirical constraints for (n0, E0, , K) and (Sv, L) can then be confronted with one another. Only the latter fall into the empirical range for (n0, E0). However, judging the extent a nuclear potential reproduces empirical (saturation) properties can be quite misleading without taking meaningful uncertainties into account; especially, the truncation of the EFT expansion at a finite order can result in sizable EFT truncation errors (even at N 3 LO) that need to be quantified.
See Ref. (65) for a review of the link between the nuclear EOS and nuclear observables; e.g., from measurements of the isoscalar giant monopole resonance K ≈ 210 − 230 MeV was inferred.
Until a few years ago, the prevalent way of estimating theoretical uncertainties in nuclear matter calculations was parameter variation within some (arbitrary) range; that is, probing the observable's sensitivity to, e.g., the 3N LECs or momentum cutoff. Recently, the focus has been more on the systematic quantification of EFT truncation errors (66), which can be estimated by assuming that an observable's EFT convergence pattern at order k takes the form y k (n) = y ref (n) k m=0 cm(n)Q m (n) (67). Here, y ref (n) sets a dimensionful reference scale, Q(n) is the dimensionless expansion parameter, and the cn(n) are the dimensionless coefficients not to be confused with the LECs of the interaction (e.g., y4 = E/A at N 3 LO). Note that c1 = 0 in Weinberg power counting. For given choices of y ref (n) and Q(n), the c m k (n) are obtained from order-by-order calculations {y0, y1, . . . , y k } of the observable. Since y ref (n) and Q(n) factor in all physical scales, the cm(n) are expected to be of order one (i.e., natural ), unless the coefficients are fine-tuned. The standard EFT uncertainty, which assumes that the truncation error is dominated by the first omitted term, has been implemented by Epelbaum et al. (68) and applied to a wide range of observables in finite nuclei and infinite matter. This "EKM uncertainty" can be summarized at N j LO as δy(n) = y ref Q j+2 max(|c0|, |c1|, . . . , |cj+1|), whose point-wise estimates can be interpreted as Bayesian credibility regions under a particular choice of priors for cm (13).
Both y ref (n) and Q(n) need to be chosen; e.g., y ref y 0 and Q ∼ p/Λ b with typical momentum p ∝ kF(n) and some estimate of the EFT breakdown scale Λ b have been applied to the nuclear EOS.
The Bayesian Uncertainty Quantification: Errors in Your EFT (BUQEYE) collaboration has recently introduced a Bayesian framework for quantifying correlated EFT uncertainties in the nuclear EOS (69,59). In contrast to the standard EFT uncertainty, the new framework allows for the quantification and propagation of statistically meaningful uncertainties to derived quantities (e.g., the pressure) while accounting for correlations across densities and between observables. Without considering these correlations, uncertainties can be overestimated. The framework also includes Bayesian model checking tools (70) for diagnosing and testing whether the in-medium ChEFT expansion works as assumed (e.g., inference for Λ b ). Gaussian Processes (GPs) with physics-based hyperparameters are trained on the order-by-order calculations of the energy per particle under the assumption that all cm(n) are random curves drawn from a single GP (67). The Gaussian posterior for the cm(n) is then used to estimate the to-all-orders EFT truncation error δy k (n) = y ref (n) ∞ m=k+1 cm(n)Q m (n) and combined with additional (e.g., many-body) uncertainties. From the posterior also arbitrary derivatives in n can be obtained.   (25) found that the corresponding binding energies (charge radii) of medium-mass nuclei are predicted too small (too large) compared to experiment, in (dis)agreement with the expectations from SNM. Since both observables were also much more sensitive to the 3N LEC cD, SNM and medium-mass nuclei seem more intricately connected than one might naively expect (25). Figure 3 shows the order-by-order predictions for the energy per particle, pressure, and speed of sound squared in PNM at the 1σ confidence level based on "GP-B 500" (69). The observables show an order-by-order convergence pattern at n 0.1 fm −3 , whereas N 2 LO and N 3 LO have a markedly different density dependence at n n0 due to repulsive 3N contributions. This is also manifested in the Bayesian diagnostics (67). Assuming Q = kF/Λ b , the inferred breakdown scale Λ b ≈ 600 MeV is consistent with free-space NN scattering and could be associated with n > 2n0. The EFT truncation errors are strongly correlated in density and to those in SNM. A correlated approach is therefore necessary to propagate uncertainties reliably to derived quantities, although the standard EFT uncertainty for the energy per particle is broadly similar to the 1σ confidence level (69).

Nuclear symmetry energy and the isospin-asymmetry expansion
The nuclear symmetry energy is a key quantity to understand the structure of neutronrich nuclei and neutron stars. Although masses of heavy nuclei constrain the value of the symmetry energy well at nuclear densities, its density dependence is much less known (73). Studying the density-dependent symmetry energy from theory, experiment, and observation is therefore an important task in the era of multimessenger astronomy.  Overall, the constraints from ChEFT are consistent with each other, even at the highest densities shown, but the uncertainties in Esym(n) are generally sizable, e.g., 20.7 ± 1.1, 31.5 ± 3.0, and 49.0 ± 12.0 MeV at n0/2, n0, and 2 n0, respectively, for Lim et al. at the 1σ confidence level. Drawing general conclusions from comparing the sizes of these bands can be misleading since the underlying methods for estimating uncertainties are quite different. Order-by-order comparisons for a wider range of chiral NN and 3N interactions with EFT truncation errors quantified are called for to provide more insights in and stringent constraints on Esym(n). The Bayesian statistical tools introduced by the BUQEYE collaboration allow for such systematic studies.
Despite the large uncertainties in the SNM EOS (see Sec. 3.1), predictions for Esym(n) [as an energy difference] can be made with significantly smaller uncertainties than those in E(n, 1)/A and E(n, 0)/A individually, if correlations are properly considered. Reference (59) found that the EFT truncation errors associated with the PNM and SNM calculations in Ref. (24) are highly correlated, meaning that the uncertainty in Esym(n) is less than the usual in-quadrature sum of errors. Combined with order-by-order calculations up to N 3 LO this led to the narrow bands "GP-B 500" in Fig. 4 (left) Fig. 4 (left)) are in excellent agreement with the joint experimental constraint ["GP-B 450" is slightly shifted to higher (Sv, L)], indicating a stiffer neutron-rich matter EOS at n0 compared to the other theoretical constraints. This is, however, consistent with joint theory-agnostic posteriors from pulsar, gravitational-wave, and NICER observations (e.g., compare with Figure 1 in Ref. (77)). An important feature of the correlated GP approach is that the theoretical uncertainties in n0 (including truncation errors) are accounted for through marginalization over the Gaussian posterior for the saturation density predicted from the SNM calculations, pr(n0) ≈ 0.17 ± 0.01 fm −3 . Apart from "HK" allowing slightly lower (Sv, L), all shown theory calculations satisfy the constraint (solid black line) derived from the conjecture (78) that the unitary gas (UG) sets a lower bound for the PNM EOS. Overall, Fig. 4 (right) shows that current constraints from nuclear theory and experiment predict the symmetry energy parameters in the range Sv  Fig. 4 (right)).
While the standard quadratic approximation [E(n, δ)/A = E(n, 0)/A + Esym(n) δ 2 ] is in general sufficient to characterize the isospin-asymmetry dependence of the nuclear EOS, certain neutron-star properties, such as the crust-core transition (79) and the threshold for the direct URCA cooling process (80), are sensitive to nonquadratic contributions. Neglecting charge-symmetry breaking effects, the energy per particle may be assumed to have an expansion in the asymmetry δ of the form E(n, δ)/A ≈ E(n, 0)/A+ L l=1 S 2l (n) δ 2l , where the standard quadratic approximation corresponds to S 2l>2 (n) = 0. Note, however, that already the FFG contributes to the nonquadratic terms, e.g., S FFG Recently, however, Kaiser (82) has shown that MBPT at second order gives rise to additional logarithmic contributions ∼ δ 2l ln |δ| with l 2. Furthermore, Wellenhofer et al. (63) found that the analogous expansion of the free energy exhibits convergent behavior for δ 1 only at high temperature. That is, the expansion's radius of convergence decreases to zero in the limit T → 0 (with diverging S 2l>2 ), as implied by the logarithmic terms at T = 0. Nevertheless, Wen and Holt (83) demonstrated that the coefficients of the normal and logarithmic terms at T = 0 can be extracted up to O(δ 6 ) from high-precision MBPT calculations with chiral interactions. Such calculations allow for the improvement of existing parametrizations in δ at T = 0 and help motivate the investigation of alternative schemes, such as an expansion in terms of the proton fraction x = np/n = (1 − δ)/2 for neutron-rich matter.

Nuclear thermodynamics
While thermal effects are negligible in isolated neutron stars, they become important in neutron star mergers and core-collapse supernovae, where T 100 MeV can be reached. Dense matter at such high temperatures not only consists of nucleons and leptons but also of additional particles such as pions and hyperons. The consistent inclusion of these particles in medium is work in progress (85,86). In the nascent field of multimessenger astronomy, one of the immediate theoretical needs is consistent modeling of (i) cold neutron stars, (ii) hot hypermassive neutron stars formed in the aftermath of compact object mergers, and (iii) core-collapse supernovae so that observations and simulations in any one of these astrophysical regimes can be propagated to constrain the others. Finite-temperature MBPT with ChEFT interactions is a suitable framework for this purpose, and here we describe some of the results on nuclear thermodynamics in recent years (for reviews, see Refs. (87,88)).
The salient thermodynamic feature of homogeneous nuclear matter at sub-saturation densities is the presence of a liquid-gas type instability toward the formation of clustered structures. In neutron stars, this instability corresponds to the crust-core transition, involving such intricate features as a variety of "pasta" shapes (89). The nuclear liquid-gas instability is also connected to the observed multifragmentation events in intermediate-energy heavy-ion collisions. In the idealized case of (infinite) nuclear matter, there is a liquid-gas phase transition of van-der-Waals type. Nuclear matter calculations at finite temperature with chiral interactions have provided predictions for the properties of this phase transition, in particular the location of the critical point. Figure 5 (left) shows the second-order MBPT results for the boundary of the liquid-gas coexistence region (so-called binodal) of SNM obtained in Ref. (35). The predicted critical point, especially the associated temperature Tc ≈ 17 − 19 MeV, is consistent with estimates (e.g., Tc ≈ 15 − 20 MeV (90)) extracted from multifragmentation, nuclear fission, and compound nuclear decay experiments (90,84).
The binodal of SNM was recently computed using the SCGF method (91) and with lattice EFT (38); the results are similar to those of Fig. 5.
In the interior of the binodal a region where the homogeneous system is unstable with respect to infinitesimal density fluctuations can be found. The boundary of this region is called spinodal. Between the binodal and spinodal the uniform system is metastable. [The two boundaries coincide at the critical point.] For SNM, the unstable region is identified by a negative inverse isothermal compressibility, κ −1 T = n(∂P/∂n) < 0. An equivalent stability criterion is ∂µ/∂n > 0, corresponding to a strictly convex free energy density F(T, n) as a function of n. If charge-symmetry breaking effects are neglected, SNM can be treated as a pure substance with one particle species (nucleons), whereas ANM is a binary mixture with two thermodynamically distinct particles (neutrons and protons). This implies that the stability criteria are different in the two cases, and for ANM the region with κ −1 T < 0 is a subregion of the spinodal region. There are various equivalent stability criteria for binary mixtures (92). A useful criterion is that outside the spinodal the free energy F(T, nn, np) is a strictly convex function of nn and np (see Ref. (52) for details). The MBPT results for the asymmetry dependence of the critical temperature Tc(δ) from Ref. (52) are shown in Fig. 5 (right). For comparison, we also show the trajectory of the temperature Tκ T (δ) where the subregion with negative κ −1 T vanishes. The trajectory of Tc(δ) reaches its T = 0 endpoint at a small proton fraction x; i.e., while PNM is stable at all densities, already small x lead to a region where the system undergoes a phase separation (52,93).
While for SNM the binodal corresponds to the Maxwell construction, for ANM the more involved Gibbs construction is required (93).
A useful characteristic for the temperature dependence of the nuclear EOS is the thermal index Γ th (T, n, δ) = 1 + P th (T, n, δ)/E th (T, n, δ), where P th is the thermal part of the pressure, and E th is the thermal energy density. For a free gas of nucleons with effective masses m * n,p (n, δ) one obtains for Γ th the temperature-independent expression The thermal index of a free Fermi gas is Γ th,free = 5/3, [To be precise, for δ / ∈ {0, 1} the above expression is valid only in the classical limit, but it provides a good approximation to Γ th (n, x) for intermediate values of δ as Ref. (61) showed.] Recently, Refs. (22,94) showed that Γ th (n, δ) reproduces the exact Γ th with high accuracy. This implies that the temperature-dependence of the EOS can be characterized in terms of a temperature-independent effective mass (see Ref. (61) for a recent implementation), which is in particular useful for monitoring thermal effects in astrophysical applications (95, 96).

Applications to neutron star physics
In this Section our goal is to emphasize the prominent role of nuclear theory in modeling neutron stars, core-collapse supernovae, and neutron star mergers. We begin by placing high-energy nuclear astrophysics in the more general context of the QCD phase diagram and discuss under what ambient conditions ChEFT can serve as a tool to constrain the properties of hot and dense matter. Specific applications include the neutron-star massradius relation, moment of inertia, and tidal deformability, as well as the nuclear EOS and neutrino opacity for astrophysical simulations.

Scales in hot and dense stellar matter
The extreme astrophysical environments found in core-collapse supernovae, neutron star interiors, and neutron-star mergers span baryon number densities nB ∼ 10 −7 − 10 1 n0, temperatures T ∼ 0 − 100 MeV, and isospin asymmetries δ ∼ 0 − 1 (corresponding to electron fractions Ye ∼ 0 − 0.5) (97). In Sections 2 and 3 we have shown that ChEFT provides a suitable framework to constrain the EOS, transport, and response properties of hadronic matter when the physical energy scale is well below the chiral symmetry breaking scale of Λχ ∼ 1 GeV. In practice, ChEFT descriptions of nuclear matter based on highprecision NN and 3N forces begin to break down at densities n ≈ 2 − 3 n0 and temperatures T 30 MeV. Therefore, additional modeling is needed at high densities and temperatures to cover all regions of astrophysical interest. For this purpose, high-energy heavy-ion collisions at RHIC, LHC, and especially FAIR aim to probe states of matter similar to those that exist naturally in neutron stars, but reaching sufficiently large proton-neutron asymmetries remains a significant challenge that may be addressed with next-generation radioactive ion beam facilities, such as FRIB. The interplay of microscopic ChEFT, whose convergence pattern is not especially sensitive to the isospin asymmetry, together with upcoming nuclear experiments that create and study hot, dense, and neutron-rich matter, will provide a direct line of inquiry probing neutron-star physics from low to high densities.
From the observational side, measurements of neutron star masses, radii, tidal deformabilities, and moments of inertia are expected to place constraints on the pressure of betaequilibrium matter at n 2n0 (98,74,99). In Fig. 6, we present a qualitative overview of the QCD phase diagram and highlight regions probed by nuclear experiments (RHIC, LHC, FAIR, and FRIB), theory (lattice QCD and ChEFT), and astrophysical simulations of neutron stars, supernovae, and neutron star mergers. We see that ChEFT intersects strongly with the region of FRIB experiments and nuclear astrophysics, providing a bridge between new discoveries in the laboratory and their implications for neutron stars. The next decade is expected to witness a strong interplay among all of these different fields, with nuclear theory predictions getting confronted with stringent empirical tests.

Neutron star structure
The mass-radius relation of non-rotating neutron stars is determined from the EOS by the general relativistic equations for hydrostatic equilibrium, the Tolmann-Oppenheimer-Volkoff (TOV) equations: where r is the radial distance from the center of the star, M (r) is the mass enclosed within r, ε is the energy density, and p is the pressure. Analysis of spectral data from neutron stars in quiescent low-mass x-ray binaries and x-ray bursters (100,101) have resulted in radius measurements R1.5 = 10−13 km for typical 1.5 M neutron stars. More recently, the NICER FRIB Figure 6: Schematic view of the QCD phase diagram. We highlight regions probed by experiments (RHIC, LHC, FAIR, and FRIB), regions of validity for lattice QCD and ChEFT, and environments reached in neutron stars, supernovae, and neutron star mergers.
x-ray telescope has observed hot spot emissions from the accretion-powered x-ray pulsar PSR J0030+045. Pulse profile modeling of the x-ray spectrum from two independent groups have yielded consistent results for the neutron star's mass M = 1.44 +0.15 −0.14 M (102) and M = 1.34 +0. 15 −0.16 M (103) and radius R = 13.02 +1.24 −1.06 km (102) and R = 12.71 +1.14 −1.19 km (103) at the 68% credibility level. Future large area x-ray timing instruments, such as STROBE-X and eXTP, have the potential to reduce uncertainties in the neutron-star mass-radius relation to ∼ 2% at a given value of the mass. This would significantly constrain the neutron-rich matter EOS at n ≈ 2n0 and when combined with mass and radius measurements of the heaviest neutron stars could give hints about the composition of the inner core (104). In recent years numerous works have studied constraints on the neutron star EOS from ChEFT. In Ref. (105) the EOS of neutron-rich matter was calculated up to saturation density with MBPT using chiral NN and 3N interactions. To extrapolate to higher densities, a series of piecewise polytropes was used to parameterize the EOS. It was found that ChEFT generically gives rise to soft EOSs that lead to 1.4 M neutron stars with radii in the range R1.4 = 10 − 14 km. Subsequent studies (e.g., Refs. (106,107,61)) have employed a wider range of chiral forces, increased the assumed range of validity for ChEFT calculations to 2n0, and explored other high-density EOSs, including smooth extrapolations and speed of sound parameterizations. The choice of transition density at which ChEFT predictions are replaced by model-dependent high-density parameterizations has a particularly large influence on neutron-star radius constraints. For instance, when the transition density was raised to nt = 2n0, Ref.  to R1.4 > 13 km. To demonstrate how a precise neutron-star mass and radius measurement can constrain the EOS of beta-equilibrium matter at n = 2n0, in Fig. 7 we show the correlated probability distribution (74) for the radius of a 1.4 M neutron star and the nuclear symmetry energy at twice saturation density Esym(2n0). In the inset we show the conditional probability distribution for Esym(2n0) assuming a precise measurement of R1.4 = 12.38 km. For the specific EOS modeling used in Ref. (74), such a precise radius constraint determines Esym(2n0) with an uncertainty of approximately 10%.
In addition to a high-density extrapolation, a uniform-matter EOS from ChEFT needs to be supplemented with a neutron-star crust model, e.g., the BPS crust model (109).
In addition to radius measurements, there has long been the possibility (110,111) to obtain a neutron-star moment of inertia measurement based on long-term radio timing of PSR J0737-3039, a binary pulsar system in which the periastron advance receives a small correction from relativistic spin-orbit coupling. A recent analysis (112) has shown that by 2030 a moment of inertia measurement of PSR J0737-3039A to 11% precision at the 68% confidence level is achievable. The moment of inertia for a uniformly rotating neutron star of radius R and angular velocity Ω can be calculated in the slow-rotation approximation, valid for most millisecond pulsars, by solving the TOV equations together with where λ and ν are metric functions andω is the rotational drag. In Refs. (113,114) the moment of inertia of PSR J0737-3039A, which has a very well measured mass of M = 1.338 M , was calculated from EOSs based on ChEFT. In Ref. (113) it was found that at the 95% credibility level, the moment of inertia of J0737-3039A lies in the range 0.98 × 10 45 g cm 2 < I < 1.48 × 10 45 g cm 2 , while Ref. (114) found a consistent but somewhat larger range of 1.06 × 10 45 g cm 2 < I < 1.70 × 10 45 g cm 2 . The moment of inertia is strongly correlated with the neutron star radius, and it has been shown (115) that measurements of the PSR J0737-3039A moment of inertia can constrain its radius to within ±1 km.
The radius and moment of inertia of a typical 1.4 M neutron star is strongly correlated with the value of Esym(n ≈ 2n 0 ) (98).
In the past ten years, several neutron stars (116,117,118) with well measured masses of M 2 M have been observed. The maximum mass (M TOV max ) of a non-rotating neutron star is a key quantity to probe the composition of the inner core, which must have a sufficiently stiff EOS to support the enormous pressure due to the outer layers. To date the strongest candidate for the heaviest measured neutron star is PSR J0740+6620, with a mass of M = 2.14 +0.20 −0.18 M at the 95% credibility level (118). As mentioned previously, ChEFT generically gives rise to relatively soft EOSs just above nuclear saturation density. The existence of a very massive neutron star with M = 2.14 M would require a stiff EOS at high densities, revealing a slight tension with ChEFT (119). However, even smooth extrapolations (74,61) of EOSs from ChEFT can produce maximum neutron star masses in the range 2.0 M M TOV max 2.4 M , and therefore more precise radius measurements (or the observation of heavier neutron stars) are needed to make strong inferences about the EOS in the ChEFT validity region n 2n0.
Beyond M TOV max , additional stable branches (120), such as hybrid quark-hadron stars or pure quark stars, may appear before the ultimate collapse to a black hole.

Neutron star mergers
The advent of gravitational wave astronomy has opened a new window into the visible Universe. Current gravitational wave detectors (LIGO and Virgo) are sensitive to frequencies 10 Hz < f < 10 kHz, which is the prime range for compact object mergers and supernovae. Gravitational wave astronomy therefore has major implications for the field of nuclear astrophysics (3). In particular, during the late-inspiral phase of binary neutron star coalescence, a pre-merger neutron star will deform with induced quadrupole moment Q under the large tidal gravitational field E: Qij = −λEij, where λ is the dimensionful tidal deformability parameter. Tidal deformations enhance gravitational radiation and increase the rate of inspiral. Gravitational wave detectors are sensitive to such phase differences and hence the dense matter EOS, but such corrections enter formally at fifth order (121) in a post-Newtonian expansion of the waveform phase and are therefore difficult to extract. The tidal deformability is an important observable in its own right, but this quantity is also strongly correlated with both the neutron star radius (122), since more compact stars experience a smaller deformation under a given tidal field, and especially the moment of inertia through the celebrated I-Love-Q relations (123). The post-merger gravitational wave signal from binary neutron star coalescence can also carry important information on the nuclear EOS. It has been shown (124) that the peak oscillation frequency f peak of a neutron-star merger remnant is strongly correlated with neutron star radii. Moreover, a strong first-order phase transition can show up as a deviation in the empirical correlation band between f peak and Λ (125).

LIGO: Laser Interferometer Gravitational-Wave Observatory
The first observation (126), GW170817, of a neutron star merger through its gravitational wave emissions was accompanied by a short gamma-ray burst and optical counterpart (127). The combined multi-messenger observations of this single event have resulted in a wealth of new insights about the origin of the elements and the properties of neutron stars. Analysis of the gravitational waveform resulted in a prediction (128) (126) to deviate more than 20% from the canonical value of M 1.4 M , assuming low neutron star spins. In summary, GW170817 data were found to strongly favor the soft EOSs predicted from ChEFT, though many other models (131) with generically stiffer EOSs were consistent with the upper bounds on Λ and R1.4 from GW170817.
Current gravitational wave interferometers do not have large signal-to-noise ratios at the high frequencies expected during the post-merger ringdown phase and therefore GW170817 provided no clues about the fate of the merger remnant. Nevertheless, an analysis of the spectral and temporal properties of the kilonova (132) optical counterpart to GW170817 have been used (133,134,135,136,137) to infer the lifetime of the merger remnant. Depending on the component neutron star masses prior to merger (primarily the total mass Mtot) as well as the maximum mass for a nonrotating neutron star M TOV max , the merger remnant can (i) undergo immediate collapse to a black hole, (ii) exist as a short-lived hypermassive neutron star supported against collapse by differential rotation, (iii) persist as a longer-lived supermassive neutron star supported against collapse by rigid-body rotation, or (iv) form a stable massive neutron star. While there is still some uncertainty about what ranges of Mtot will lead to each of the above four scenarios, it has been suggested (134) that prompt collapse will occur when Mtot (1.3 − 1.6)M TOV max , hypermassive neutron stars will be created when 1.2 M TOV max Mtot (1.3 − 1.6)M TOV max , and supermassive neutron stars will result when Mtot 1.2M TOV max . Each merger outcome is expected to have a qualitatively different optical counterpart and total mass ejection, since longer remnant lifetimes generically give rise to more and faster moving disk wind ejecta.
Observations of the GW170817 kilonova suggest that the most likely outcome of the neutron star merger was the formation of a hypermassive neutron star, which would imply a value of M TOV max = 2.15 − 2.35 M (135,134,137). Eliminating the possibility of prompt black hole formation in GW170817 also rules out compact neutron stars with small radii and tidal deformabilities. In Ref. (133) such arguments were used to infer that the radius of a 1.6 M neutron star must be larger than R1.6 10.7 km, while in Ref. (136) it was found that the binary tidal deformability parameter for the GW170817 event must satisfyΛ 400. Both of these inferred constraints are compatible with predictions (74,129,107,106) from ChEFT. However, the constraint on the binary tidal deformabilityΛ > 400 can rule out a significant set of soft EOSs (138), roughly half of those allowed in the analysis of Ref. (74). Combined gravitational wave and electromagnetic observations of binary neutron star mergers together with more precise radius measurements therefore have the possibility to strongly constrain the dense matter EOS and related neutron star properties in the regime of validity of ChEFT (138,139). As a demonstration, in Fig. 8 we show the joint probability distributions (74) for the pressure of nuclear matter at n = 2n0 with (i) the radius of a 1.4 M neutron star, (ii) the tidal deformability of a 1.4 M neutron star, and (iii) the moment of inertia of PSR J0737-3039A with a mass of 1.338 M .

Core-collapse supernovae
Neutron stars are born following the gravitational collapse and ensuing supernova explosion of massive stars (M 8 M ). The core bounce probes densities only slightly above normal nuclear saturation density (97) and leaves behind a hot (T ∼ 20 − 50 MeV) nascent protoneutron star. During the subsequent Kelvin-Helmholtz phase that lasts tens of seconds, the proto-neutron star emits neutrinos, cools to temperatures T < 5 MeV, de-leptonizes, and contracts to reach supra-saturation densities in the innermost core. The success or failure of the supernova itself (140), the thermal and chemical evolution during the Kelvin-Helmholtz phase (141), and the possibility of novel nucleosynthesis in the neutrino-driven wind (142) depend on details of the nuclear EOS and weak reaction rates.
Investigating the qualitative impact of specific EOS properties, such as the incompressibility or the symmetry energy, on the fate of supernova explosions is often challenging due to "Mazurek's Law", a colloquial observation that feedback effects tend to wash out any fine tuning of parameters in core-collapse supernovae (143). Nevertheless, several recent systematic investigations (144,95) of EOS parameters have supported the idea that a high density of states, linked to a large value of the in-medium nucleon effective mass M * , reduces thermal pressure and leads to enhanced contraction of the initial proto-neutron star. This results in the emission of higher-energy neutrinos that support the explosion through the neutrino reheating mechanism (140). Since microscopic calculations based on ChEFT tend to predict larger values of the effective mass than many mean-field models (145), these observations help motivate recent efforts (146) to include thermal constraints from ChEFT directly into supernova EOS tables.
Neutrino reactions also affect the nucleosynthesis outcome in neutrino-driven wind outflows and the late-time neutrino signal that will be measured with unprecedented detail during the next galactic supernova. Charged-current neutrino-absorption reactions νe + n → p + e − andνe + p → n + e + , which can be calculated from the imaginary part of vector and axial vector response functions, are especially sensitive (147,148) to nuclear interactions and in particular the difference ∆U = Un − Up between proton and neutron mean fields. The isovector mean field is especially important in the neutrinosphere, the region of warm and dense matter where neutrinos decouple from the exploding star.
Recently, the calculation (149,145) of nuclear response functions that include mean-field effects from ChEFT interactions have shown that terms beyond the Hartree-Fock approximation are needed for accurate modeling. In particular, resummed particle-particle ladder diagrams were shown (149) to produce larger isovector mean fields due to resonant, nonperturbative effects in the NN interaction. Moreover, for neutral-current neutrino reactions, such as neutrino pair bremsstrahlung and absorption, resonant NN interactions were shown to significantly enhance reaction rates at low densities compared to the traditional one-pion exchange approximation (150).

Summary and outlook
In this article, we have reviewed recent progress in ChEFT calculations of nuclear matter properties (with quantified uncertainties) and their implications in the field of nuclear astrophysics. Combined with observational and experimental constraints, these microscopic calculations provide the basis for improved modeling of supernovae, neutron stars, and neutron star mergers. In particular, we have highlighted MBPT as an efficient framework for studying the nuclear EOS and transport properties across a wide range of densities, isospin asymmetries, and temperatures. We have also shown how advances in high-performance computing have enabled the implementation of two-and multi-nucleon forces in MBPT up to high orders in the chiral and many-body expansions. Finally, we have described new tools for quantifying theoretical uncertainties (especially EFT truncation errors) to confront microscopic calculations of the nuclear EOS with empirical constraints. Such systematic studies are particularly important in view of EOS constraints anticipated in the new era of multimessenger astronomy; e.g., from gravitational wave detection, mass-radius measurements of neutron stars, and experiments with neutron-rich nuclei.
In the following we briefly summarize several open research directions at the interface of nuclear EFT and high-energy nuclear astrophysics. (i) EFT truncation errors and their correlations in density and across observables need to be studied with different many-body frameworks and nuclear interactions at arbitrary isospin asymmetry and finite temperature. (ii) Together with EFT truncation errors, the uncertainties in the LECs parametrizing the interactions need to be quantified and propagated to nuclear matter properties using a comprehensive Bayesian statistical analysis. (iii) The full uncertainty quantification of the nuclear EOS will be aided by the development of improved order-by-order chiral NN and 3N potentials and the study of different regularization schemes as well as delta-full chiral interactions. (iv) Many-body calculations of nuclear matter properties beyond the nuclear EOS (e.g., linear response and transport coefficients) with chiral NN and 3N interactions are required for more accurate numerical simulations of supernovae and neutron star mergers.
(v) Future neutron star observations will provide stringent tests of nuclear forces and nuclear many-body methods in a regime that is presently largely unconstrained. The interplay between observation, experiment, and theory in the next decade can be expected to result in many further advances in our understanding of strongly interacting matter.

DISCLOSURE STATEMENT
The authors are not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this review.