Green's Function Techniques for Infinite Nuclear Systems

I review the application of self-consistent Green's function methods to study the properties of infinite nuclear systems. Improvements over the last decade, including the consistent treatment of three-nucleon forces and the development of extrapolation methods from finite to zero temperature, have allowed for realistic predictions of the equation of state of infinite symmetric, asymmetric and neutron matter based on chiral interactions. Microscopic properties, like momentum distributions or spectral functions, are also accessible. Using an indicative set of results based on a subset of chiral interactions, I summarize here the first-principles description of infinite nuclear systems provided by Green's function techniques, in the context of several issues of relevance for nuclear theory including, but not limited to, the role of short-range correlations in nuclear systems, nuclear phase transitions and the isospin dependence of nuclear observables.


INTRODUCTION
The recent discoveries of neutron-star binaries GW170817 [1] and GW190425 [2] are formidable feats in gravitational-wave (GW) and multimessenger astronomy. These events provide unique insight into a great variety of astrophysical questions, from stellar physics to cosmology. In dense matter physics, the analysis of multimessenger data is providing a unique way forward in our understanding of matter in its most neutron-rich form. Gravitational waves from the inspiral carry imprints of the equation of state (EoS) of neutron-star matter [3]. The initial analysis of the source GW170817 has already provided significant constraints on the masses, radii and EoSs of neutron stars [4]. The expectation is that, as more of these GW sources are detected in the next few years, they will provide more and more accurate observational constraints. Simultaneous developments in X-ray burst observations are also providing equally important and independent constraints on the EoS and the physical properties of neutron stars (NSs) [5]. Soft X-ray waveforms of rotation-powered pulsars also provide promising NS constraints [6,7].
Theoretical nuclear physics tools provide relevant input for dense matter physics. Among others, there have been several significant efforts trying to put constraints on the symmetry energy [8,9] or the EoS itself [10]. The data in the new generation of astrophysical observations will tighten these constraints, and may even allow for further relevant conclusions on macroscopic observables, like the EoS [11]. Connecting these bulk observable constraints to the strong interaction among neutrons requires significant efforts on the nuclear theory side. There is a stringent need for the development of realistic predictions with systematic errors based on our understanding of the physics of dense systems. This necessarily requires the use of ab initio nuclear theory techniques that can provide a solid, quantifiable link between macro-and micro-physics [12,13]. For the case of GW inspirals, this link also demands for an understanding of dense matter properties not only in the zero-temperature domain, but at temperature scales which are significant enough to modify some nuclear properties (e.g., temperatures in the range of tens of MeV) [14][15][16][17][18]. There is a similar need for microphysics inputs in the case of core-collapse supernovae, where the formation process of proto-NSs [19] and the subsequent development of a neutrino-heated shock wave require an understanding of neutrino emission from compact objects [20,21].
The accuracy of theoretical nuclear physics predictions is hampered by a variety of different factors. First, the strong interaction among nucleons is difficult to pin down. Scattering experiments provide extremely useful insight at the two-body level, but three-neutron forces remain relatively unconstrained [22]. For neutrons, four-body and higher level forces are expected to play a small role for the EoS [23]. Chiral effective field theory (χEFT) [24,25] and Lattice QCD [26,27] efforts in the next decade will be crucial to provide clear limits on these terms. As a first approximation, one can estimate the uncertainty in the strong force by considering different input Hamiltonians that fit scattering phase-shifts and few body data [28,29].
In addition, even if the interaction is perfectly understood, one still needs to solve a many-body multicomponent 1 problem to link the NN force to the macroscopic properties related to the EoS. While several methods have been devised to attack this problem, a clear-cut solution is difficult, if not impossible. The treatment of many-body correlations differs depending on the method that is used to describe dense neutron-rich systems [30,31]. In particular, the structure of the strong force itself-particularly its short-range core-may preclude low-order perturbative treatments. Self-consistent Green's function (SCGF) techniques can be used to obtain both perturbative and nonperturbative results from these interactions of relevance for nuclear physics [32][33][34]. This approach offers the promise to simultaneously tackle issues that are relevant for neutron-star astrophysics, while providing insight into problems that are appropriate for nuclear structure (particularly in terms of shortrange correlations).
My aim with this review is to provide a summary of the recent advances of the SCGF treatment of infinite nuclear systems. I will only explore the uncertainty associated to the Hamiltonian by looking at calculations performed with three different chiral interactions. This represents a small subset of possible forces, and is by no means representative of the overall uncertainty due to the strong interaction in infinite systems. Even then, these interactions predict relatively different micro-and macroscopic properties below, around and above saturation density. A benchmark analysis with different methods, while necessary to fully understand the limitations of nuclear theory, lies beyond the scope of this contribution. I will instead provide a qualitative, and somewhat limited discussion of the many-body scheme dependence of the results, by providing comparisons to manybody perturbation theory results based on similar Hamiltonians.
I start this review with a brief summary of the formal tools associated to the SCGF formalism in the following section. I 1 Multicomponent in the sense that it requires both neutrons and protons. note that there are monographs [32,34] and books [33,35,36] that provide many more details, and I refer the interested reader to those. Section 3 is devoted to a discussion of recent numerical results that concern single-particle properties in dense matter. These include spectral functions, self-energies and momentum distributions. In section 4, the discussion turns into the macrophysics of infinite nuclear systems, with a particular emphasis on the EoS. I finally provide some conclusions and an outlook in section 5.

Formalism
Many-body Green's function are a natural starting point in the discussion of quantum many-body systems [33,37]. Green's function (or propagators) arise naturally in diagrammatic treatments, and provide a description of several real or virtual excitation processes within the many-body system. Formally, these propagators take the form of time-ordered correlation functions associated to the creation and annihilation of excitations and can be expressed in terms of Feynman diagrams [38,39]. Broadly speaking, Green's function techniques attempt to sum series of the most physically relevant diagrams to describe the system under consideration. For instance, in the case of thermodynamics of correlated infinite systems, one often resorts to an expansion that involves renormalized interactions treating the short-range core nonperturbatively via infinitely repeated scattering processes [40]. When it comes to the response of finite systems, one instead considers diagrams that describe relevant excitations within the system [32]. This allows SCGF practitioners to look at many-body systems from different perspectives, providing a characterization of the system that includes microscopic and macroscopic properties.
In SCGFs methods, internal propagators of all diagrams are consistently dressed. In other words, renormalization effects are introduced from the onset in the formalism -and hence the monicker "self-consistent." At a diagrammatic level, this means that only skeleton (or 1-particle-irreducible) diagrams are required for single-particle self-energies. In practice, one pays the price of having to run simulations iteratively for a given system which, in infinite matter, translates into a fixed density (or temperature) [33]. A key advantage of SCGF methods is the direct access to the propagators themselves. These contain relevant information for the system's microphysics, in terms of spectroscopy and strength distributions [32]. Moreover, propagators can also be used to compute bulk properties, including total energies, thermodynamical properties and pressures. Within a single SCGF simulation, one therefore obtains a significantly broad scope in terms of physical information. This goes beyond the reach of other many-body approaches which are mostly tailored to address the energy of the system, either by minimizing it using quantum Monte-Carlo [31,41] techniques, or addressing it by means of diagrams, like the Brueckner-Hartree-Fock (BHF) approach [30].
Another strategic advantage of the SCGF formalism is its generality. Within the very same formalism, one can tackle both finite and infinite systems-see [42] for a contribution in this volume that specifically discusses recent applications of SCGF techniques to finite nuclei. More importantly in the context of astrophysics, SCGF methods can be formulated consistently at finite temperature [33,37]. In fact, fullyfledged infinite nuclear matter studies are usually performed at finite temperature [43][44][45]. For sufficiently large temperatures (typically T 2 − 5 MeV), this treatment avoids the complications associated to nuclear superfluidity [43,46] that plagued early attempts of SCGF calculations in the 1990s [47]. Furthermore, a finite temperature formulation allows for access to astrophysically-relevant temperature dependences of observables [18], and provides insight into liquid-gas, thermallyinduced phase transitions in dense systems [48,49]. It is now possible to use these non-superfluid finite temperature calculations to perform meaningful extrapolations to zero temperatures for both microscopic [50] and bulk properties [51], as we shall see in the following.
A key quantity in the SCGF formalism is the so-called spectral function A k (ω). This characterizes fully the one-body propagator G, which is generally defined as a time-ordered product of a creation and a destruction operator [33,36]. Upon a Hilbert transform with a complex energy variable z, the spectral function and the one-body propagator are related via a spectral decomposition, Different time orderings give rise to different components of the propagator, including the retarded and advanced components that can be addressed by taking z = ω ± iη, with η infinitessimal [35]. All time-ordered components can be accessed through the spectral function A k (ω), which also has the advantage of having a probabilistic interpretation. The spectral function provides the probability distribution associated to either extracting or attaching a particle with well-defined momentum k, and energy ω, to the infinite system. In the finite-temperature Matsubara formalism, the Lehmann representation of A k (ω) takes into account not only the thermal excitations on the ground state of the system, but also the presence of excited states that are populated according to a thermal distribution, In this equationâ k (â † k ) are destruction (creation) operators of single-particle states with well-defined momentum k. A manybody system with N n particles can have different eigenstates of the many-body Hamiltonian with energy E n , labeled with quantum number n. Z = n e −β(E n −µN n ) is the corresponding partition function and, in the grand-canonical ensemble, the system is described in terms of a given temperature T = 1/β, and a chemical potential µ.
As a well-defined probability distribution, the spectral function is normalized dωA k (ω)/2π = 1. When this integral is weighted by the corresponding thermal population of hole states, one obtains the momentum distribution of the system, In turn, even though the formulation of the problem is grand-canonical in nature, the momentum distribution can be used to fix the density of the system by requiring that the chemical potential µ in the Fermi-Dirac function is fulfilled. Here ν = 2 represents the degeneracy in spin for neutron or asymmetric matter, and ν = 4 is the degeneracy in spin and isospin for symmetric matter. The Dyson equation [z − k 2 /2m − k (z)]G k (z) = 1 provides a connection between the one-body propagator G and the selfenergy . The latter is traditionally decomposed into a real and an imaginary part, and the spectral function is then connected to the retarded self-energy by the equation: The real part of the self-energy Re k (ω), is often interpreted as an energy-dependent mean-field potential [52]. Systems with a well-defined quasi-particle structure have spectral functions with strong peaks around the energy for which the denominator of Equation (5) is smallest. This happens at the so-called quasiparticle energy, which can be accessed from Re k (ω) by solving the self-consistent equation at each momentum k. This on-shell condition provides a one-to-one correspondence between energy and momentuma dispersion relation that is often used to characterize relevant single-particle properties. Other many-body approaches, like many-body perturbation theory or the BHF diagrammatic method, make use of such on-shell dispersion relations throughout [39]. The SCGF method, in contrast, considers all energy-dependent, offshell components in the self-energy and the spectral function. These effects are relevant to describe the fragmentation of single-particle strength [33]. The imaginary part of the self-energy is usually discussed in terms of quasi-particle lifetimes [33]. For the retarded selfenergy, the imaginary component is negative, and is directly related to the available phase space. Indeed, for an energyindependent self-energy, the Fourier transform of Equation (5) to the time domain provides a propagator that decays exponentially in time with a characteristic timescale, The energy and momentum dependence of the self-energy can be characterized in some domains by a model-independent analysis [33,52]. According to Luttinger's theorem for normal fermionic systems, at zero temperature the imaginary part of the self-energy should be exactly zero at the Fermi energy and, close to this energy, it should behave quadratically, so that Im k (ω = µ) ≈ a k (ω−µ) 2 [53]. When temperature sets in, Im k at ω = µ grows quadratically with T [37]. Expressions for the self-energy can be derived at different orders in many-body perturbation theory [33,36]. I will present results here based on the ladder approximation [33]. A summary in terms of diagrams of this approximation is provided in Figure 1. The key ingredient in this approach is the T−matrix shown in diagram (c), which is obtained from a Lippman-Schwinger equation where the intermediate kernel involves a particle-particle and hole-hole two-body propagator at finite temperature. This energy-dependent in-medium interaction describes re-scattering effects in the medium and can therefore capture short-range effects. Three-nucleon forces are included in the ladder approximation by constructing an effective two-body force, diagram (a) of Figure 1. This is the result of adding up the bare two-nucleon force (dashed line), V, and a three-body force (dotted line) that has been appropriately averaged over a third nucleon [54][55][56].
In the ladder approximation, the imaginary part of the selfenergy is given by the expression: distribution and ImT is the imaginary part of the T−matrix.
In SCGF methods, one generally computes the imaginary part first, which tends to be simpler to evaluate thanks to energy and momentum conservation. Whereas explicit expressions for the real part are normally available too, one can avoid a series of complex integrals over momentum and energy phase-space by computing the real part from the dispersion relation, where P denotes the principal value. The energy-independent (but momentum-dependent) component ∞ k , is akin to a Hartree-Fock self-energy. This effective one-body interaction is built by averaging the two-nucleon and the three-nucleon interaction over the Fermi sea of one or two particles, as depicted in diagram (b) of Figure 1 [54]. This average includes propagator renormalization effects via a correlated momentum distribution, In this expression, the matrix elements correspond to antisymmetrized two-nucleon (V) and three-nucleon (W) interactions, respectively. The chiral interactions that are used in this work are discussed in more detail in the following section. The self-energy operator can be defined diagrammatically and encodes different many-body processes depending on the level of approximation [33]. At the lowest order, the socalled mean-field or Hartree-Fock approximation, the selfenergy is energy-independent and formally given by the same expression as in Equation (10), but with internal momentum distributions n k , computed from on-shell Fermi-Dirac distributions n k = f (ε k ). Second-order self-energies are instead energy-dependent, proportional to V 2 and have already a non-trivial imaginary part [33,51]. These account for 2p-1h and 1p-2h excitations within the system. In infinite matter applications, non-perturbative approximations beyond the second order are typically considered. A self-energy arising from a ladder resummation of the in-medium interaction deals effectively with the strong short-range and tensor components of the NN force, and leads to stable numerical results even for hard core interactions-for details see [57,58].
When a SCGF calculation converges for a given set of external parameters, one typically stores the complex self-energy k (ω). All other microscopic properties, including spectral functions, can be derived from it according to Equation (5). The energy per particle of the system is then accessed by an energy-weighted integral over the spectral function, the so-called Galitskii-Migdal-Koltun sum-rule [33], Here, we also incorporate a correction term proportional to the expectation value of the 3NF W . This is computed by using an additional integral over momentum of the second term in Equation (10). In other words, we keep the lowest order approximation to the expectation value of the 3NF, in agreement with techniques typically employed in finite nuclei [59]. The formal procedure to go beyond this approximation exists, as discussed in [54], but has not been implemented in infinite systems yet. At finite temperature, the system is described in terms of thermodynamical potentials like the Helmholtz free energy F = E − TS. This requires an explicit calculation for the entropy S, based on Green's function methods. A proposal to compute S from spectral functions in dense matter was put forward in [58,60] based on the formal techniques derived by Pethick and Carneiro [61,62]. With access to the entropy, one can then derive all thermodynamical properties, including the pressure p = ρ(F/A − µ). Formally, the SCGF method is known to give rise to thermodynamically consistent results [33,63]. For instance, a zero-temperature calculation of the chemical potential µ, from the Fermi energy of the system ε k=k F , or from the density derivative of the free energy density, provides the same result. This formal agreement has also been demonstrated numerically in the past [45,60].
In terms of formalism, there is no inconvenience to extend the SCGF method to asymmetric, multicomponent systems. In fact, the SCGF method has already been used to describe isospin asymmetric matter in the past [64,65]. One can therefore characterize correlations and fragmentation in isospinimbalanced systems [66,67]. In asymmetric matter, all quantities become functions of isospin-i.e., the spectral functions for neutrons A n k (ω) and protons A p k (ω) are different, as are the self-energies and all other one-body properties. In addition to the isospin splitting, one must also consider changes in the effective interactions, so that neutron-neutron (nn), neutronproton (np), proton-neutron (pn), and proton-proton (pp) inmedium T-matrices are explicitly different. Extensions of the normal ordering procedure to compute effective one and twobody forces from 3NF matrix elements are also necessary [56], but relatively straightforward [68]. The results for asymmetric matter presented here within the SCGF formalism include explicitly the effect of isospin asymmetry in the normal-ordering of 3NFs and in the extrapolation procedure from finite to zero temperature.

Interactions
For consistency in the presentation, all the results that follow have been generated with the same set of interactions. I will consider three sets of χEFT-inspired forces. At the NN level, these three interactions have been introduced in [69]: they have been fitted to the same NN scattering data and incorporate different cutoffs in the relative momentum: = 414 MeV, 450 MeV and 500 MeV. These forces also differ slightly in the shape (particularly the sharpness) of the regulator in relative momentum. 3NFs arise naturally in the χEFT formalism [24,25]. Reference [70] introduced a set of N2LO 3NF forces based on the previous three NN interactions. The associated 3NF low-energy constants c D and c E were fitted to the binding energies of A = 3 nuclei and to the 3 H-3 He Gamow-Teller transition matrix element [70]. Calculations are performed for two-body matrix elements with partial waves up to J = 9.
The 3NFs are included in the SCGF simulations by building effective one-and two-body density-dependent forces using a procedure that is akin to normal ordering [71,72], but includes a renormalization of occupation probability in n k . Diagrams (a) and (b) in Figure 1 show the corresponding one-and two-body effective interactions. Diagrams (c) and (d) indicate how these are used in the T−matrix calculation and the corresponding self-energy, respectively. This method is rooted on an extension of the SCGF formalism to multi-body forces [54], and has been implemented in symmetric nuclear matter and neutron matter simulations in [55,56]. We resort to some numerical approximations in this averaging procedure. First, the normalordering procedure ignores the center-of-mass dependence of the 3NF, which is set to P = 0. Second, matrix elements off the diagonal in relative momentum q = q ′ , are extrapolated from diagonal matrix elements using the prescription q → 1 2 (q 2 + q ′2 ). This procedure also involves a 3NF regulator, which is chosen , and is based on the two-body relative momentum q and the internal, integrated single-particle momentum p. The exponent n is chosen in accordance to the regulator in the 2NF [70]. Within the P = 0 approximation, the off-shell prescription provides results of very similar quality to many-body perturbation theory calculations that incorporate fully-fledged 3NFs from the outset [73].

Extrapolation Procedure to Zero Temperature
Finite-temperature calculations are useful on their own, but benchmark comparisons with other methods are typically performed in the zero-temperature limit [30,31]. A zerotemperature extrapolation procedure for SCGF simulations was developed in the context of beyond-BCS pairing [50,74]. This procedure is now available generically and can be used to provide zero-temperature data for symmetric and asymmetric systems. For a given fixed density, the extrapolation requires as input a few (typically 4 − 10) finite-temperature simulations. One typically performs simulations from high temperatures of order 20 MeV, and subsequently dials down the temperature to values that are closer to zero. Too close to zero, though, the simulations will not converge as one enters into the superfluid regime -so an intermediate range must be found. It is also important to note that the dependence on temperature typically scales with the dimensionless degeneracy parameter T/ε F , where ε F is the Fermi energy. Because ε F typically increases with density, the expansion is closer to the zerotemperature result (and therefore more accurate) at higher than at lower densities.
The extrapolation from finite temperature to zero temperature works in three independent steps. In a first step, the selfenergies from different temperatures are remeshed to the same values of momentum and energy. In the second step, the temperature dependence of the real and imaginary parts of the self-energy (for each value of k and ω) are fitted to a polynomial function in temperature. The third step is to use not one, but different polynomial fits to extrapolate to T = 0. A series of tests are run on these different extrapolations, and a quality measure is built based on a series of consistency checks. For instance, the renormalization factor ] −1 is compared to the size of the discontinuity in the momentum distribution at k = k F . In a consistent theory, the two quantities should provide the same value [33]. The quality measure penalizes fits where these values are significantly different. Importantly, some of these fits have a built-in parabolic energy dependence of Im around ω = µ. Additional details of this extrapolation procedure are discussed in [50].
I provide an example of this extrapolation procedure in Figure 2. This shows the temperature dependence of three properties that are representative of the system's behavior as a function of temperature at a fixed density of ρ = 0.16 fm −3 .

Calculations have been performed with the N3LO
= 500 MeV NN force, supplemented with 3NFs at the N2LO level. The averaging procedure for these 3NFs consistently takes into account the density and temperature dependence of the internal propagators. Qualitatively similar results are obtained for the = 450 and 414 MeV interactions. The energy dependence of the spectral function A k=0 (ω) and the imaginary part of the self-energy k=0 (ω) are displayed in the left and central panels, respectively. I focus on the k = 0 value for simplicity, but the calculation and extrapolation procedure are performed for all available momenta. The effect of temperature in the spectral function is concentrated in the region near ω = µ. The quasi-particle peak and the positive and negative highenergy tails are insensitive to thermal effects. For the self-energy, there is also a relatively strong temperature dependence near ω = µ, which is known to scale with T 2 [37]. The phasespace dominated peaks in the particle and hole domains at energies about 200 MeV away from µ are also sensitive to the temperature. The extrapolated T = 0 result is able to capture the temperature dependence of these quantities, while providing consistent results.
The right panel of Figure 2 shows the associated momentum distribution n k , obtained under the same conditions. Thermal effects provide additional depletion at low momenta and populate moderately the vicinity of the Fermi surface. In fact, temperature modulates substantially the momentum distribution around k = k F even for relatively low values of T ≈ 4 − 5 MeV. The extrapolation procedure, however, captures the energy-and FIGURE 3 | Spectral function A k (ω) as a function of energy for three different momenta: k = 0 (left) k = k F (central) and k = 2k F (right). SCGF calculations are performed at a density ρ = 0.16 fm −3 and have been extrapolated to zero temperature. Results for N3LO NN plus N2LO 3NFs have been performed with three different cutoffs. The vertical line and solid dot indicate the position of the quasi-particle peak ω = ε k , for each of these momenta. The height of this peak has been chosen arbitrarily.
momentum-information sufficiently well to produce a sharp transition at the Fermi surface, as expected from Luttinger's theorem [37]. The size of this discontinuity is also in agreement with the renormalization parameter associated to the self-energy itself. This extrapolation procedure has been implemented not only in the symmetric nuclear matter and pure neutron matter [74], but also in asymmetric infinite matter, as I shall discuss in the following sections.

MICROSCOPIC PROPERTIES
In the following, I provide an overview of the microscopic properties predicted by the SCGF approach for symmetric, asymmetric and pure neutron matter. The aim of this section is to show the broad scope of information that can be generated in SCGF simulations. I focus for simplicity on a single value of the density ρ = 0.16 fm −3 , which is representative of both nuclear systems and of neutron stars. Note however that the method is customarily used at other densities. The majority of the results are presented in their zero-temperature extrapolated form, although the temperature dependence of most quantities is available and has been used for extrapolation purposes. I leave a full analysis of the density, isospin asymmetry and temperature dependence of these results for future studies. The general features of the spectral strength in infinite systems are well-understood [44,47,57,75]. At all momenta other than k F , the spectral strength has a prominent quasi-particle peak. In infinite nuclear matter, the peak is very close to the quasi-particle energy ω = ε k obtained from Equation (6). The quasi-particle energy for the = 500 MeV interaction is displayed in the three panels of Figure 3 with vertical lines. The solid dot is shown for illustrative purposes, and its height is chosen arbitrarily. The quasi-particle peak is well-defined in infinite nuclear systems, in the sense that the peak of A k (ω) and ω = ε k agree relatively well. At this density, the quasi-particle peak is independent of the interaction under consideration. In fact, the cutoff dependence of the quasi-particle peak is not resolvable on the scale of this graph. I note, however, that while the peaks are well-defined, they are also relatively broad. Typical widths at k = 0 and k = 2k F are of the order of Ŵ k ≈ 40 MeV.

Spectral Functions
Spectral functions in infinite nuclear matter display slowly decaying high (positive and negative) energy tails. These are particularly prominent at the Fermi surface k = k F , where the zero-temperature spectral function should be a combination of a broad energy-dependent background and a zero-width δ function centered at ω = ε k F [33]. In practice, the T = 0 results at k = k F do show some minor structures around ω = µ as the non-analytical behavior of the δ function is not built in the extrapolation procedure. The results of Figure 3 clearly show a substantial cutoff dependence in these high-energy components. This is relatively unsurprising, since these tails are a direct consequence of the short-range core in the NN force [32]. As expected, the interaction with the largest momentum cutoff ( = 500 MeV), and hence the largest resolution in real space, provides the largest high-energy components. These include a shoulder at positive energies in the region 100−500 MeV which is characteristic of a relatively strong short-range core. In contrast, the lower cutoff interactions provide much faster (exponential) decaying spectral functions beyond the quasi-particle peak. For positive energies, they provide very little single-particle strength beyond 500 MeV. The = 414 MeV interaction has faster decaying tails than the 450 MeV force, which agrees with intuition. Note that this sharply decaying spectral strength is at odds with calculations employing traditional phenomenological NN forces, which typically populate these offshell regions with higher probabilities [57,76]. Off-shell high-energy components are observed in electron-scattering experiments [77] and may provide a way to quantify the short-range component of NN forces [78].
The results presented so far concerned symmetric nuclear matter, where by construction one assumes the same fraction of neutrons and protons. Isospin asymmetric systems, however, are interesting on their own. The deep interior of heavy stable and unstable isotopes is expected to be a neutron-rich environment, characterized by an isospin asymmetry parameter, of order η ≈ 0.2 [79]. In contrast to this relatively mild isospin asymmetry, the core of neutron stars is expected to be extremely asymmetric with η ≈ 0.8 [80]. Finite-temperature SCGF simulations including explicitly isospin asymmetry have been available for over a decade [64]. In these simulations η is an external, tuneable parameter. SCGF calculations then give isospin-dependent self-energies, spectral functions and quasiparticle properties, as well as the corresponding bulk properties. These calculations provided insight into the isospin-asymmetry dependence of the fragmentation of the single-particle strength [66,67], but were somewhat limited by two factors. First, finite temperature effects have a larger importance in the minority component as asymmetry is turned on, because the corresponding degeneracy parameter T/ε τ F decreases. Second, the calculations were performed with two-body interactions only Calculations with the N3LO = 500 MeV and an N2LO 3NF at a density ρ = 0.16 fm −3 have been extrapolated to zero temperature for different isospin asymmetries η. The vertical dashed line is a guide to separate hole (ω < µ τ ) and particle (ω > µ τ ) components. and thus missed the effects of 3NFs. None of these factors is particularly critical. Temperature effects can be controlled, and the effect of 3NFs on the spectroscopic strength is small below saturation density for both symmetric nuclear and pure neutron matter [55,56].
Having said that, recent technical developments have allowed us to overcome these two shortcomings recently. Figure 4 shows the spectral functions of neutrons (top panels) and protons (bottom panels) as a function of energy for different values of the isospin asymmetry η, at a fixed density of ρ = 0.16 fm −3 . Results are presented for the = 500 MeV N3LO potential supplemented with 3NFs 2 . The results for the isospin dependence of these spectral functions are qualitatively similar to previously obtained results which lacked 3NFs. The evolution of the neutron spectral function from symmetric (η = 0) to neutron (η = 1) matter is relatively mild. The quasi-particle peaks with respect to the corresponding chemical potential µ n become slightly more (attractive) repulsive for the (hole) particle components. The high energy tails are relatively insensitive to the isospin asymmetry, and are mostly affected for the negative energy components at momenta above the Fermi surface. All in all, this agrees with the generic picture of a neutron majority component evolving from a symmetric situation, where tensor correlations provide larger correlations, to neutron matter, where only the short-range component is responsible for fragmentation.
In contrast, the proton spectral function (bottom panels) is much more affected by isospin asymmetry. In asymmetric conditions, the relative importance of np pairs increases for the minority proton component. Since np interactions are mediated by the tensor component of the NN force, one naively expects 2 I note here for completeness that the η = 0 and η = 1 calculations are performed with independent codes that do not include isospin asymmetry explicitly, and hence provide a good numerical test of the arbitrary isospin simulations. enhanced correlation effects for protons in neutron-rich matter. Simultaneously, for neutron-rich conditions the overall proton density decreases, so the proton component should somehow become "less interacting." SCGF simulations allow us to provide quantitative predictions to distinguish between these scenarios.
As the isospin asymmetry increases, the single-particle peak for protons becomes more repulsive with respect to µ p . Importantly, the width of the quasi-particle peak for components below the Fermi surface (see k = 0 panel) also decreases substantially. As a consequence, positive high-energy tails increase with isospin. Interestingly, the central panel (k = k p F ), shows a background component that behaves differently for the addition and removal energy components. Whereas for ω < µ p , the background of the spectral function decreases with isospin asymmetry, at ω > µ p it increases. Above the Fermi surface (k > k p F ), the quasi-particle peak is more compressed with increasing isospin asymmetry, both in terms of a lower peak and a significantly decreased width. The limit case of a proton impurity in an η = 1 system cannot be tackled with the present set of SCGF technology.
While a full analysis is beyond the scope of this initial presentation of results, the isospin dependence shown here is relevant for both nuclear experiments and astrophysics. Hadronic and electron two-nucleon knock-out reactions from the early 2000s up to the present have shown that np pairs, rather than pp or nn pairs, are more likely to be knocked off isospin symmetric targets [81,82]. Further electron-scattering results in isospin asymmetric nuclei have provided additional evidence for this neutron-proton dominance [83]. It is likely that the tensor component preferentially acts in relatively highmomentum pairs and provides the physical mechanism behind these strong isophobic correlations [66,84]. Extending the np dominance picture to isospin imbalanced systems would naively indicate that protons should be more correlated than neutrons in neutron-rich systems, although quantitative statements depends sensitively on the definition of a correlation measure [85][86][87]. In the astrophysical context, the strong modification of the proton spectral function in isospin asymmetric nuclear matter predicted by SCGF simulations has been explicitly explored at the level of the symmetry energy in [88][89][90]. The fragmentation of strength in neutron matter has also been used to provide a description of the pairing properties in dense matter [50,74].

Self-Energies
Within the SCGF formalism, the spectral strength is directly related to the self-energy. The imaginary part of the retarded self-energy encodes information on both the phase space and the NN interaction, and can even be used to diagnose interactions by means of sum rules [76]. I present a density plot of the zero-temperature extrapolated imaginary part of the self-energy at ρ = 0.16 fm −3 in Figure 5. Left, central, and right panels correspond to chiral interactions with cutoffs = 414, 450, and 500 MeV, including in all cases 3NFs. As expected, Im k is zero at ω = µ throughout all momenta. For energies ω < µ, this function has a hole component with support below k ≈ k F = 263 MeV. The size and shape of this pocket is relatively interaction-independent. In contrast, the positive energy (ω > µ) components of Im k show a clearer cutoff dependence. The deep pocket around k = 0 in this component is sensitive to the underlying NN force. Typical phenomenological interactions have deep pockets, and these extend to very high positive energies (into the GeV domain). For this set of chiral interactions, the depth of the pocket is relatively insensitive to the cutoff and of order 50 − 70 MeV. The tails of Im k=0 fall off relatively quickly with energy. In fact, the = 414 and 450 MeV self-energies are almost zero beyond an energy of about 500 MeV. The cutoff also has a clear effect on the momentum dependence of Im k in the positive energy region. Generally speaking, higher-cutoff forces allow the imaginary part of the self-energy to extend to higher momenta.
The solid lines in Figure 5 show the trajectory of the quasi-particle peak in the energy-momentum plane. Quasiparticle lifetimes τ k in Equation (7), are computed along these trajectories. The plot indicates that the lifetimes are relatively insensitive to the differences in self-energies. In other words, largely different offshell self-energies can give rise to relatively similar on-shell properties. On-shell components will be discussed in detail in section 3.4.
The isospin asymmetry dependence of Im τ k (ω) is explored in Figure 6. Results are shown for the = 500 interaction only at a fixed density of ρ = 0.16 fm −3 , for a variety of isospin asymmetries. The leftmost panel corresponds to the symmetric system (η = 0), and is the same as the right panel of Figure 5. The top panels indicate that the effect of asymmetry in the neutron self-energy is relatively mild. The pocket at negative energies ω < µ n , associated to neutron holes, is relatively unchanged as η increases. In contrast, the positive energy (particle) components evolve with isospin in a substantially different way. The pocket that in the symmetric situation decayed relatively quickly with momentum, becomes an almost momentum-independent structure in the most asymmetric case explored here.
In contrast, proton self-energies show a strong dependence on isospin asymmetry. As the proton fraction decreases, the negative-energy (hole) components of Im p k disappear steadily. This reflects the decrease in density, and hence of available phase space of protons. The corresponding low-momentum (k < k p F ) inverse quasi-particle lifetime is small, which in turn gives rise to an increasingly narrow spectral function. This is clearly in agreement with the results in the bottom left panel of Figure 4. On the positive energy side of the self-energy ω > µ p , the deep pocket at low momenta becomes noticeably deeper as η increases but, unlike the neutron component, the momentum dependence  remains relatively unchanged. The increase in depth of Im k can be interpreted in terms of an increase in the interaction strength, as expected in the asymmetric case with an increase in the number of np pairs. In other words, the imaginary component of the self-energy for protons in asymmetric matter shows the competition of the two effects discussed earlier in two different energy regimes. On the one hand, proton holes, associated to the negative energy components, are dominated by phase-space and their self-energy decreases steadily as the proton density decreases. On the other hand, proton particles at positive energies show an increase in interaction strength as expected in a more imbalanced strongly interacting system. The real part of the self-energy is not shown here for brevity. I note, however, that this is connected to Im τ k by the dispersion relation in Equation (9).

Momentum Distributions
A useful characterization of beyond mean-field correlations is obtained with the momentum distribution in Equation (3). For a non-interacting Fermi gas or within the lowest-order Hartree-Fock approximation at zero temperature n k = 1 (0) for momenta below (above) the Fermi surface. This behavior is illustrated by the double-dotted-dashed line in the left panel of Figure 7. In contrast, the momentum distribution of correlated infinite nuclear matter shows a characteristic depletion below the Fermi surface. This depletion is of the order of 7 − 12%, with small quantitative differences between different interactions.
In turn, a depletion of momentum eigenstates below k F means that n k above the Fermi surface must be populated. The logarithmic scale in the right panel of Figure 7 illustrates the effects of such high-momentum components. These components are noticeably different depending on the cutoff. Smaller cutoff interactions give rise to momentum distributions that decay faster and start decaying at lower momenta, than larger cutoff forces. To illustrate the relatively small effect that 3NF have on correlations, I show in Figure 7 the momentum distribution obtained in a SCGF calculation without 3NFs (dashed-dotted line). Remarkably, the difference between the calculation with and without 3NFs is substantially smaller than the cutoff dependence of the interaction. This bodes well with previous observations indicating that the effect of 3NFs on fragmentation properties is quantitatively small, even though they have a substantial influence on the energetics and thermodynamics of the system [55].
The evolution of the momentum distribution with isospin also provides an illustrative understanding of the interplay between short-range and tensor correlations and isospin asymmetry. Figure 8 shows the momentum distribution as a function of momentum for different isospin asymmetries. These results have been obtained with the = 500 MeV interaction, but qualitatively similar plots are found for the other two forces. The left panel shows the results for neutrons. Isospin asymmetry has two main effects on the majority component. On the one hand, the discontinuity in n k moves to higher momenta, as the Fermi momentum changes from k F = 263 MeV in symmetric nuclear matter to k F = 331 MeV in pure neutron matter. This effect would be identical in the free Fermi gas (FFG). In contrast, the distinct increase of n k below the Fermi surface is specific to correlated systems. The depletion of strength decreases (or gets closer to 1) linearly with the increase of asymmetry. Whereas in symmetric nuclear matter it is around 10%, in pure neutron matter it goes down to a level of around 4 − 5%. In fact, this isospin dependence is rather universal and independent of the NN interaction under consideration [66,67]. A similar effect is expected for the size of the discontinuity across the Fermi surface, which decreases as isospin asymmetry increases.
The right panel of Figure 8 provides insight on the asymmetry evolution of the minority component (protons). Mirror effects are obtained in this case. The discontinuity of n k occurs at smaller and smaller momenta, as expected from FFG considerations. The depletion below the Fermi surface departs further away from 1 as the isospin asymmetry increases. In other words, by this measure, protons (or minority components) become more correlated in a neutron-rich environment. Close to the astrophysically relevant condition at η ≈ 0.8, protons are depleted by about 15% for this specific interaction. This depletion goes in hand with a change of high-momentum components for the proton, which become more important as asymmetry increases [67].
The picture obtained by these complete calculations including zero temperature extrapolations and 3NFs is qualitatively and quantitatively very similar to the one discussed in [66,67]. Overall, the effect of isospin asymmetry on occupation numbers is relatively small. For the characteristic isospin asymmetry in the interior of nuclei (η = 0.2), the change in occupation below the Fermi surface is of the order of only a few percent. Overall, this shallow asymmetry dependence is at odds with some of the past and recent experimental analysis in both electron-and hadron-induced knock-out reactions [85,87,91]. Importantly, however, these calculations at fixed density do not account for the density and isospin-asymmetry dependence within nuclei, which is likely to be relevant for the description of nuclear momentum distributions.

On-Shell Properties
A unique advantage of SCGF simulations is the access to both offshell, energy-dependent quantities, as well as the corresponding information on the energy shell-at the quasi-particle energy ω = ε k . On-shell properties are often used in interpreting the properties of the system and are the input to most astrophysical simulations. In particular, quasi-particle potentials in asymmetric matter at relatively low densities (neutrino-sphere) are important for charged-current neutrino opacities [20,92]. The quasiparticle spectrum of Equation (6) is a key quantity that allows for a characterization of single-particle excitations in the system. Since the kinetic term in the quasi-particle energy obscures the many-body effects at large momentum, I instead focus on the single-particle potential U τ k = Re τ k (ε τ k ). This is shown as a function of momentum in Figure 9. Top panels correspond to neutrons, whereas bottom panels correspond to protons. From left to right, the isospin asymmetry of the system changes from symmetric nuclear matter (η = 0) to pure neutron matter (η = 1) at a constant density of ρ = 0.16 fm −3 . I show results for the = 500 (solid), 450 (dashed) and 414 MeV (dotted lines) interactions.
As discussed previously, these results are obtained from a zero-temperature extrapolation of finite temperature SCGF simulations. As an illustration of the importance of thermal effects, I also show in the top left and top right panels the = 500 MeV results for U τ k at temperatures up to T = 20 MeV. The thermal effects are clearly small on the scale of this figure, particularly in comparison to the cutoff and the isospin dependence of the results. I note in passing that the temperature dependence is not monotonous. At momenta below the Fermi surface U τ k becomes more attractive as temperature increases. Above the Fermi surface and up to about 500 MeV, the potential instead becomes more repulsive with temperature. This is in agreement with previously reported SCGF results that did not incorporate 3NFs [44].
In terms of isospin dependence, the results also agree with previously reported outcomes in the SCGF approach [64], and in other many-body calculations [93][94][95]. In symmetric nuclear matter, the single-particle potential has a minimum at low momenta close to U k=0 ≈ −60 MeV, and increases steeply with momenta. At k = 0, the lowest cutoff interaction produces a potential that is about 5 MeV more attractive than the highest cutoff force. This cutoff-dependence is relatively constant as a function of momentum up to k ≈ 600 MeV, where regulator artifacts start to dominate. Above this value of As asymmetry increases, the neutron single-particle potential becomes more repulsive and shallower in terms of momentum dependence. In contrast, protons experience a larger attraction when they become the minority species. Their single-particle potentials become significantly more attractive, and in the astrophysically relevant region of η ≈ 0.8, they reach values of 80 − 90 MeV at k = 0. The isospin dependence is to a good approximation linear with isospin asymmetry, as expected on general grounds [94]. I note however that the momentum and isospin dependence of these single-particle potentials is significantly different from those predicted by phenomenological mean-field approaches [96].
Interestingly, the cutoff dependence of the results depends on isospin asymmetry and on the specific nucleon. As neutrons become the predominant species and their single-particle potentials become shallower, the cutoff dependence in the lowmomentum single-particle potential below 400 MeV decreases. In neutron matter, the differences between the potentials are of the order of 1 MeV. On the contrary, as the isospin asymmetry increases and the proton spectrum becomes more attractive, the differences between interactions with different cutoffs increase and become of the order of more than 10 MeV.
These differences can be explained by two interrelated factors. One is of a relatively straightforward nature. np interactions play no role in neutron matter, and hence the intrinsic differences in the NN interaction are less likely to appear as several partial waves channels are suppressed. In other words, if the fits to the nn channels have intrinsically smaller differences between forces than the corresponding np channels, neutronrich systems will show less cutoff dependence. Another factor that plays an important role is the strength (in the sense of "perturbativeness") in different channels. In terms of χEFT, the small cutoff dependence in neutron matter agrees with the idea that neutron matter is a more perturbative system than symmetric nuclear matter [69,97,98], because the strong np components are missing. It is important to emphasize here that the strong cutoff dependence of the impurity spectrum may not be identified in the energetics of the system, but will clearly have an impact on the dynamics of isospin asymmetric systems.
The momentum dependence of U τ k is directly related to the concept of the effective mass. SCGF simulations can be used to extract not only the mass associated to the single-particle spectrum, but also to the so-called ω− and k−effective masses. A detailed analysis of these properties within the SCGF approach is ongoing.
The on-shell imaginary part of the self-energy, Equation (7), is interpreted as an inverse quasi-particle lifetime. For a normal correlated many-fermion system at zero temperature Ŵ k=k F should be zero, thus signaling an infinite lifetime of quasiparticles at the Fermi surface. The lifetime should be non-zero everywhere else and, as discussed in the context of section 3.2, its magnitude is a combination of phase-space and interaction effects. Figure 10 provides an overview of the momentum dependence of Ŵ τ k for neutrons (top panels) and protons (bottom panels) as a function of isospin asymmetry. As in the previous figure, extrapolated zero-temperature results are shown.
The isospin and cutoff dependence of the results is very different above and below the Fermi surface. Below the Fermi surface, the SCGF simulations predict a non-zero width. Unlike in BHF calculations, intermediate hole states account explicitly for this non-zero width of hole states k < k F [94].   Around the Fermi surface, for all isospin asymmetries, the momentum dependence of Ŵ τ k is rather universal. In most cases, one cannot distinguish the curves from different interactions within about 100 MeV of the Fermi surface. I demonstrate again the importance of thermal effects for the = 500 MeV interaction at the two extreme asymmetries η = 0 (top left) and η = 1 (top right panel). Whereas regions away from the Fermi surface are temperature independent, these panels show that the temperature dependence in Ŵ τ k is prominent around k = k τ F . There, the quasi-particle width switches from zero at T = 0 to a finite (and increasing) value with temperature. This is to be expected, as the behavior around k F is governed by phase space arguments [33,53].
The particle component (at k > k τ F ) of Ŵ τ k is remarkably more cutoff dependent. The momentum dependence immediately above k τ F and up to about 500 − 600 MeV is very similar in all cases, and shows an increase from zero up to an isospindependent maximum for the = 450 MeV and 414 MeV cases. For neutrons (protons), the peak height decreases (increases) from around 30 MeV in symmetric nuclear matter η = 0, to around 25 MeV (50 MeV) in pure neutron matter η = 1. For these two forces Ŵ τ k becomes a decreasing function of momentum above k ≈ 500 MeV and eventually vanishes above k ≈ 1 GeV. This is stark contrast to the = 500 MeV interaction, where the peak of the particle component occurs at much higher momenta in the vicinity of k ≈ 1 GeV.
The peak is also relatively independent of the isospin asymmetry, and very similar for neutrons and protons. This peak and, in general, the relatively large values of Ŵ τ k in this region are indicative of the short lifetime of high-momentum states. The cutoff dependence clearly indicates that this component of Ŵ τ k depends strongly on the interaction, and the isospin dependence suggests that it is mediated by isoscalar short-range physics independently of the density of the individual components. These results match with the findings for the spectral functions reported in Figure 4. Whereas the width of the neutron spectral function is qualitatively unaffected by a change in isospin asymmetry, the width of the proton spectral function changes significantly with asymmetry for the hole states k < k τ F . Singleparticle widths provide an alternative way of characterizing beyond mean-field correlations, since in the Hartree-Fock approximation Ŵ k = 0. The results in Figure 10 indicate that some of the features of Ŵ k are robust against isospin asymmetry changes, whereas others are not. In particular, the relative strong dependence of proton single-particle widths with isospin asymmetry may be related to the isospin dependence of the absorptive component of optical potentials [99].

Other Microscopic Properties
The microscopic results presented so far provide a very complete picture of the dynamics of infinite nuclear systems. The versatility and wide range of applications is an advantage of SCGF techniques. On-shell and off-shell quantities provide valuable information which, when complemented with modern interactions based on χEFT, can help provide estimates of systematic theoretical errors in simulations. A useful characterization of in-medium nuclear properties can be obtained from energy-weighted moments (or sum-rules) of the spectral function. The first moment is connected to effective single-particle energies [100,101] and has been studied in the context of both symmetric and asymmetric nuclear mater at zero and finite temperatures [102][103][104]. The conclusions reached by these studies in the early 2000s indicated that the first moment is sensitive to NN interactions, and provides a useful characterization of correlations. The second moment has only been studied more recently [76]. It measures fluctuations around the single-particle peak, and is a connected to the energy dependence of the imaginary part of the self-energy. Because the latter is also particularly sensitive to the NN force, the second moment can provide a measure of the Hamiltonian properties in both nuclei and infinite nuclear matter.
The SCGF method can be used to study a variety of inmedium quantities which are relevant for nuclear physics experiments. For instance, in the late 1990s, SCGF methods were used to provide a formal redefinition of scattering phaseshifts in the medium [105]. Practical calculations also allowed for predictions of in-medium scattering cross sections in dense matter [106]. The finite lifetime of quasi-particles away from the Fermi surface can be formally translated into an in-medium nucleon mean-free path. This mean-free path is relevant for a variety of nuclear physics considerations, from the validity of the shell model [107] to the characterization of multi-fragmentation reactions [108]. The consistent description of the real and imaginary parts of the self-energy via dispersion relations makes SCGF techniques particularly appealing in this context. In [109], the mean-free path of nucleons in dense infinite matter was computed by means of an extension of the SCGF formalism to the complex plane. Calculations agree well with results obtained with other methods and, most importantly, with a variety of nuclear experiments at intermediate energies [109].
Superfluidity in nuclear systems can be detected by the appearance of in-medium bound states in the two-body scattering matrix, or the two-body propagator. This can be used to characterize the normal-superfluid phase diagram of dense matter. Extensions of the SCGF method to the superfluid domain are possible, by means of the symmetry-breaking Gorkov formalism [43,46]. Simulations in the normal phase, however, can already provide relevant information for superfluid properties, by extrapolating normal self-energies and spectral functions to zero temperature as discussed above. With this data, one can build effective gap equations that contain the effects associated to the fragmentation of strength [110]. Results with the three chiral interactions presented here indicate that short-range correlations reduce the pairing gap, and provide a gap closure at lower density compared to BCS solutions for both the singlet and the triplet channel [74]. In principle, SCGF calculations at different orders can also be used to estimate systematic uncertainties in this method. This approach may provide a better handle on the role of superfluidity at intermediate and high densities, where the role of long-range correlations is poorly understood [50,111].

MACROSCOPIC PROPERTIES
Simulations based on SCGF methods can access the thermodynamics of the system. A wide range of other many-body methods have been extended to the finite temperature domain, including diagrammatic BHF approaches [112], variational methods [113], and many-body perturbation theory techniques [68,114,115]. All of these methods share common limitations in terms of the definition of quasi-particles [116]. This can lead to an inconsistency in terms of the chemical potential. On the one hand, the chemical potential is typically obtained from a normalization condition like Equation (4). On the other, one can find µ from the thermodynamical relation where ǫ is the energy density. Formally and numerically, these two values may be different depending on the theoretical framework under consideration.
A key advantage of the SCGF formulation at finite temperature is that it guarantees that the two quantities are the same, as long as the self-energy fulfills basic criteria [33]. These criteria are related to the so-called −derivability, and lead to a thermodynamically consistent theory [63]. In practice, this means that in SCGF calculations there is no need to perform numerically the derivative in Equation (13). This is particularly convenient in calculations of the pressure, which can directly be computed from the relation p = ρ(F/A − µ). I note, however, that differences between microscopic and macroscopic chemical potentials may arise at a numerical level, particularly when the approximated treatment of 3NFs is considered [51].
I show the pressure as a function of density for symmetric nuclear matter (left panel) and pure neutron matter (right panel) in Figure 11. These have been computed with the three chiral interactions discussed in this work and include the effect of 3NFs. The results are computed at two finite temperatures T = 20 and 10 MeV, and are also extrapolated to T = 0. Note that the extrapolation uses many more temperatures, which are not displayed here for simplicity. The discussion below is brief, and details of similar calculations can be found in [44,48,49].
For symmetric nuclear matter, the zero-temperature pressure shows a characteristic van der Waals shape, with an area of negative pressure that extends from zero up to the corresponding saturation density of each interaction. This structure signals the existence of a liquid-gas phase transition in nuclear matter [48,49,112,115]. At zero temperature, the pressure vanishes at zero density and at the saturation point. Figure 11 indicates that the three chiral interactions have significantly different saturation densities, in the region ρ 0 = (0.14 − 0.165) fm −3 . As temperature increases, the spinodal region (where dp/dρ < 0) shrinks, until it disappears. At that stage, the pressure can still have an inflection point, which will smoothen out only at the corresponding critical temperature T c of the liquid-gas phase transition. Different many-body approaches based on the same NN interaction predict different critical points [48]. Similarly, the same many-body set-up with different interactions can yield different critical points [115].
The recent SCGF analysis presented in [49] suggests that the critical temperature predicted by several different χEFT interactions can range from T c ≈ 11.0 MeV to T c ≈ 18 MeV. The cutoff dependence in Figure 11 falls within this rage, and T c lies between about T c ≈ 11 MeV for the = 500 MeV interaction and T c ≈ 16.5 MeV for the = 414 MeV interaction. There is a clear correlation between the saturation energy and the critical temperature, with more attractive saturation energies leading to higher critical temperatures [49]. This intuitive behavior is not supported by the analytical predictions of mean-field models [117]. The relatively wide range of critical temperatures dominates over the uncertainties associated to the many-body approximation. Indeed, in [49], the maximum difference obtained between BHF and SCGF predictions for the critical temperature was about 2 MeV. While further many-body benchmarks may be necessary to understand the full manybody dependence of this result, it is clear that the Hamiltonian uncertainty is relevant for finite temperature predictions [51]. In fact, the cutoff dependence of a sub-saturation density property like T c can be taken as a sign that finite-temperature many-body correlations are significant in this region.
The temperature dependence of the pressure in neutron matter is relevant for astrophysical considerations, particularly in the context of neutron star mergers [16] and proto-neutron star formation [19]. The right panel of Figure 11 shows the pressure on a semilogarithmic scale as a function of density. The same three temperatures (0, 10 and 20 MeV) are presented. At low densities, the pressure should be well-described in terms of the model-independent virial expansion, which only depends on neutron-neutron scattering phaseshifts [118]. The pressure for fugacities z < 0.5 at T = 10 MeV and 20 MeV is shown with two solid bands at low densities in the figure. Clearly, the low-density SCGF results agree well with virial predictions. The low-density pressure is in fact remarkably cutoff independent, and the first distinguishable cutoff-dependent features only show up above ρ ≈ 0.10 fm −3 . This is in stark contrast to the case of symmetric nuclear matter. Noticeably, this cutoff independence also holds for the low-density, zero-temperature extrapolated results. All in all, the cutoff and many-body dependence of low-density neutron matter predictions are well under control (independently of temperature), as already discussed elsewhere [23,28,29,69,97].
The cutoff dependence in neutron matter is only clearly distinguishable above saturation density. There, and toward the high density limit of applicability of chiral interactions, the = 500 MeV interaction is providing more pressure at a given density. This is followed by a decreasing pressure from the = 414 MeV and = 450 MeV forces. It is interesting to note that the last two switch orders with respect to symmetric nuclear matter, where for a given density the = 450 MeV interaction provides more pressure than the = 414 MeV interaction. The cutoff dependence is only one of several systematic uncertainties in these results. SCGF can also provide an estimate of the manybody uncertainties, by performing calculations with self-energies obtained not only at a ladder level, but also at the first-and second-order level. The results of a recent analysis indicate that the many-body uncertainty is smaller than that associated to the underlying NN force by a factor of 2 − 3 at twice saturation density [51].
With respect to the temperature dependence, one can distinguish two distinct regimes. At sub-saturation densities, cutoff independent results dominate. At higher densities, reaching the degenerate supra-saturation regime, the cutoff dependence overtakes temperature effects. A more quantitative analysis of the temperature dependence can be obtained by looking at the so-called thermal index (or adiabatic index) of the equation of state [14]. It characterizes the temperature dependence of the pressure (or energy) of the system. Carbone and Schwenk have recently provided an in-depth analysis of the density, temperature and Hamiltonian dependence of the thermal index obtained in SCGF simulations [18]. They find that some of the assumptions made in previous literature regarding the temperature and density independence of this quantity do not hold. In addition, the effect of 3NFs is important in the high-density regime.
The two extreme limits of isospin asymmetry are interesting on their own, and provide insight into different relevant phenomena. However, both in the study of the liquid-phase transition [119] and in neutron-star astrophysics [73], the actual FIGURE 12 | Energy per particle (top) and pressure (bottom) as a function of density for different isospin asymmetries. Results have been extrapolated to zero temperature. The thin lines have been generated from Skyrme functionals fitted to second-order calculations with the same underlying NN forces by Lim and Holt [120]. The bands are obtained from fits to many-body perturbation theory calculations by Drischler et al. [73].
isospin dependence of the results is relevant. Arbitrary isospinasymmetric systems can be explored with SCGF simulations, and the results in Figure 12 provide both the energy per particle (top) and the pressure (bottom panels) as a function of density for several asymmetries. Results are presented for the three interactions with different cutoffs.
For comparison, I show two different sets of results based on similar ab initio calculations. The narrow lines, labeled Sx450 and Sx414, correspond to Skyrme parametrizations derived in [120]. These density functionals were in turn fit to asymmetric matter calculations performed at second order in many-body perturbation theory for the = 450 MeV and 414 MeV chiral interactions. I do not show the corresponding Sx500 Skyrme results, which were fitted to data with free-particle (rather than self-consistently dressed) intermediate states. The internal dressing of fermion lines is closer in spirit to the SCGF results presented here, and one thus expects to find a better agreement. I note that if non-perturbative effects were small, one would expect an agreement between the fits and the SCGF data.
The bands in the different panels of Figure 12 have been obtained from fits to second-order many-body perturbation theory calculations for seven different chiral interactions in [73]. These fits reproduce the results obtained at second order in perturbation theory with internal Hartree-Fock propagators. These calculations are different to the SCGF approach in that they implement a different averaging procedure for the 3NFs, including an average value of the center-of-mass momentum.
The results from other methods are provided here for reference and only represent a crude estimate of the manybody uncertainties. The discussion below is therefore qualitative, rather than quantitative. The top left panel of Figure 12 clearly indicates that all these interactions and approximations are able to provide saturation for symmetric nuclear matter. The SCGF = 500 MeV results saturate at a too low density and a too high energy per particle. In contrast, the 450 MeV and 414 MeV interactions predict better saturation energies, but are slightly off in terms of density. This may be due to deficiencies of the interactions, the many-body scheme or the implementation of the 3NFs into the SCGF approach [29]. I note however that the saturation energy itself is an extrapolated quantity and a recent reanalysis based on nuclear data suggests it may be a bit more repulsive than anticipated [121]. Comparing the energy per particle to the results of the Sx parametrizations, I find that the SCGF results for the 414 MeV and 450 MeV chiral interactions are qualitatively similar to those obtained in perturbation theory calculations in [70] (for symmetric nuclear matter) and [120] (for asymmetric nuclear matter). The ladder resummation in the SCGF calculations seems to bring in repulsion at intermediate densities for isospinsymmetric systems. This repulsive effect, of the order of a couple of MeV in symmetric nuclear matter, reduces as the isospin asymmetry increases.
As it is well-known, the energy per particle becomes more repulsive as η increases. For pure neutron matter, and in the subsaturation region shown in the figure, the cutoff dependence of the SCGF results is very small, certainly within 1 MeV. The size of the cutoff dependence is in agreement with the initial findings reported in [69], although the latter were obtained using thirdorder perturbation-theory calculations. In fact, as the isospin asymmetry of the system increases, the differences in energy between the Sx414 and Sx450 results and the SCGF calculations become much smaller. For densities above 0.06 fm −3 , all calculations agree in the density dependence of the energy per particle of neutron matter. This is to be expected if one assumes that neutron matter is more perturbative than symmetric nuclear matter. The SCGF and the Sx414 and Sx450 results lie in the lower range of the band provided by [73], which indicates that the cutoff dependence associated to the three chiral interactions presented in this work is not necessarily representative of the systematic uncertainty due to different Hamiltonians.
The reduced cutoff dependence of the results, as the isospin asymmetry evolves from the left to the right panels, has clear implications on the symmetry energy. Because symmetric nuclear matter is more bound for the = 414 MeV and 450 MeV forces, and the neutron matter energy is basically the same, one expects the symmetry energies of these two forces to be relatively larger than that of the 500 MeV interaction. Indeed, computing the symmetry energy as the difference between the two extremes of isospin at ρ = 0.16 fm −3 , one finds values that range from S 0 = 31.5 MeV for the 414 MeV force up to 29.8 MeV and 25.4 MeV for the 450 MeV and 500 MeV forces. These values are in agreement with those previously reported for these forces in [122], and are somewhat smaller than empirical determinations [8].
The bottom panels of Figure 12 provide a relatively similar picture in terms of the pressure. In symmetric nuclear matter (left panel), the pressure is negative at sub-saturation densities and, other than the = 500 MeV SCGF results, all results agree qualitatively. The picture evolves as isospin asymmetry increases. The spinodal, negative pressure zone decreases as η becomes larger, as expected in general grounds [119]. For most interactions, the critical asymmetry where the spinodal area vanishes is in the range η = 0.6 to 0.8. At these relatively large asymmetries, all the different many-body predictions are contained with the bands obtained by Drischler et al. [73]. Noticeably, the results of the = 500 MeV interaction that were off the band at small asymmetries, fall into the band as isospin asymmetry increases and neutron matter is approached. This indicates that the discrepancies associated to this force with respect to other simulations lie in the isovector component of the interaction.
The cutoff dependence of the SCGF results offers only an initial indication of the size of uncertainties associated to the NN force. In qualitative terms, however, the behavior is similar to that predicted by the bands in [73] which could be explored with non-perturbative methods in the future. At the level of the interaction, there are other sources of uncertainty (like the low-energy constants or the regulator dependence) that are not really considered by these exploratory considerations [51]. The dependence on the many-body method is also relevant as the system becomes more and more isospin symmetric, and quantities like the thermal index or the symmetry energy may be substantially affected by the many-body method under consideration. As such, benchmark comparisons in asymmetric systems would provide an interesting testing ground for newly developed forces.

CONCLUSIONS
In this work, I have briefly reviewed a set of SCGF techniques that have been used over the last couple of decades to study infinite nuclear systems. SCGF approaches are unique in terms of their flexibility, which allows for both perturbative and nonperturbative studies. The SCGF formalism can be formulated at finite temperature and provide thermodynamically consistent results. The same approach can be used to study the properties of finite nuclei with similar many-body approximations and equivalent interactions [42]. A single SCGF calculation yields not only bulk nuclear properties, including binding energies and thermodynamics, but also a wide range of microphysics information, ranging from predictions associated to the fragmentation of the single-particle strength to the characterization of quasi-particle lifetimes. With an appropriate account of 3NFs, this method can now competitively tackle issues in infinite nuclear matter physics with the same level of predictive power as other many-body approaches. The predictions in terms of the density and isospin dependence of fragmentation strengths are unique and provide an insight that is also relevant for the understanding of the evolution of correlations with isospin.
In this review, I have focused on the properties of symmetric and asymmetric nuclear matter and neutron matter as obtained from a set of three χEFT interactions based on N3LO NN interactions and N2LO 3NFs. This provides a limited but already illustrative set of results that spans a wide range in terms of short-range physics. Of course, one would ideally perform these simulations with the largest possible set of interactions to get a more quantitative estimate of the dependence of the results on the Hamiltonian. A more sophisticated treatment using Bayesian analysis techniques could eventually improve the estimates of systematic uncertainties [9].
The results indicate that the N3LO = 500 MeV interaction produces significantly different strength distributions than the = 450 MeV and 414 MeV forces. As expected, the 500 MeV interaction provides high momentum components in the momentum distribution, and also large additional components in the self-energy. The two softer interactions provide results which are significantly more perturbative. In other words, their momentum distributions die out quicker in momenta and the energy dependence of the associated self-energies is much more limited. These generic results hold independently of the isospin asymmetry. In fact, the results for the 450 MeV and 414 MeV interactions are very similar across a wide range of density, isospin asymmetry and temperatures. In pure neutron matter, the 500 MeV interaction results are also very close to the 450 MeV and 414 MeV interactions. This can be translated into a strong isospin dependence of the overall cutoff dependence. This point has already been discussed in the past in [69,122], and our SCGF simulations confirm these trends from a different many-body perspective.
SCGF provide access to quasi-particle properties that are relevant in various contexts. Single-particle potentials for neutrons and protons in neutron-rich matter show interesting asymmetry dependences. It is well-known that the single-particle potential for neutrons becomes less attractive with increasing isospin asymmetry. At the same time, it becomes less sensitive to the cutoff. In contrast, the proton potential becomes more attractive in neutron-rich matter and, as the isospin asymmetry increases, the cutoff dependence is enhanced. This behavior should be reflected in the corresponding effective masses. The asymmetry dependence of the single-particle width is also interesting, and again suggests substantial differences in the minority and majority components.
When it comes to bulk properties, I find results that are qualitatively close to those established with many-body perturbation theory methods. Again, the predictions of the = 500 MeV interaction are noticeably different than the lower cutoff results in the isospin symmetric case. In fact, the saturation energy of this interaction is about 4 − 5 MeV more repulsive than empirical estimates, in contrast to the lower cutoff simulations. Results for the 450 MeV-414 MeV forces are more consistent with each other. These also agree well with previous many-body perturbation theory calculations based on these interactions in isospin-symmetric, isospin-asymmetric nuclear matter and pure neutron matter [68].
One can draw several conclusions from this analysis. The isospin dependence of short-range correlations can be predicted with SCGF techniques and appears to be less pronounced in bulk systems than it is in finite nuclei [66,67,85,87]. The cutoff dependence of the results in isospin asymmetric systems affects more the single-particle properties of protons than those of neutrons. In this sense, the simulations of proton impurities in pure neutron matter could provide some insight into the modeldependence of neutron-proton forces. I also note that ab initio simulations of single-particle potentials are significantly different from those predicted by phenomenological approaches [96].
Looking forward, the SCGF method is at an advantage with respect to some other methods in its handling of finite temperature and isospin-asymmetry effects. For supra-saturation densities, this is relevant for neutron star mergers, as well as for neutron-star formation. There, the non-perturbative nature of SCGF methods can handle short-range correlation effects and hence provide meaningful results up to presumably higher densities than perturbative calculations [51]. The effects associated to the fragmentation of strength are also relevant in a variety of high-density astrophysical settings, where they have often been ignored. These include predictions for the equation of state, but also for pairing properties and response functions. Broadly speaking, previous SCGF simulations that have tackled these issues have not relied on systematically improvable chiral interactions that are also helpful in providing estimates of theoretical uncertainties. Interestingly, the low-density high-temperature regime of the neutrinosphere in supernova explosions is also relevant for astrophysical simulations, and it can be tackled with SCGF simulations. There, SCGF methods can provide insight into quantities like singleparticle potential shifts [20,92] and response functions needed for neutrino processes [21]. Our results indicate that these should be interaction-independent, and the dependence on the manybody method will provide a good handle on systematic errors.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.