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 non-perturbatively 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, fully-fledged infinite nuclear matter studies are usually performed at finite temperature [43–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, thermally-induced 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 Ak(ω). 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 Ak(ω), 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 Ak(ω) 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 a^_k (a^k†) are destruction (creation) operators of single-particle states with well-defined momentum k. A many-body system with N_n particles can have different eigenstates of the many-body Hamiltonian with energy E_n, labeled with quantum number n. Z=∑ne-β(En-μNn) 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ωAk(ω)/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 f (ω) = [1 + exp[(ω − μ)/T]]⁻¹ is such that the normalization condition

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-k2/2m-Σk(z)]Gk(z)=1 provides a connection between the one-body propagator G and the self-energy Σ. 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 quasi-particle 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 momentum—a 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 self-energy, the imaginary component is negative, and is directly related to the available phase space. Indeed, for an energy-independent 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(ω=μ)≈ak(ω-μ)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–56].

In the ladder approximation, the imaginary part of the self-energy is given by the expression:

where b(Ω) = [exp[(Ω − 2μ)/T] − 1]⁻¹ is a Bose-Einstein 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 so-called mean-field or Hartree-Fock approximation, the self-energy 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² 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=kF, 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 isospin-imbalanced systems [66, 67]. In asymmetric matter, all quantities become functions of isospin—i.e., the spectral functions for neutrons Akn(ω) and protons Akp(ω) 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), neutron-proton (np), proton-neutron (pn), and proton-proton (pp) in-medium T-matrices are explicitly different. Extensions of the normal ordering procedure to compute effective one and two-body 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.

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 ³H-³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 normal-ordering 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→12(q2+q′2). This procedure also involves a 3NF regulator, which is chosen to be non-local f(p,q)=exp[-(q2+p2/3Λ2)n], 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].

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 zero-temperature 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 zero-temperature 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 self-energies 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 Zk=[1-∂ωReΣk(ω)|ω=εk]-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⁻³. 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 Ak=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 high-energy 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² [37]. The phase-space 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 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.

Figure 3 provides an example of the structure of the zero-temperature spectral function in nuclear systems. The three panels illustrate the energy dependence of the spectral function at three different momenta: k = 0, below the Fermi surface; k = k_F, at the Fermi surface; and k = 2k_F, above the Fermi surface. Results are shown for three interactions: Λ = 414 MeV (red dotted lines); Λ = 450 MeV (yellow dashed lines); and Λ = 500 MeV (blue solid lines).

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 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 ω = ε_kF [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 quasi-particle 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 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⁻³. Results are presented for the Λ = 500 MeV N3LO potential supplemented with 3NFs². 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 μⁿ 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 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=kFp), 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>kFp), 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 high-momentum 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–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–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].

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⁻³ 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. Quasi-particle 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⁻³, 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 ω < μⁿ, 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Σkp disappear steadily. This reflects the decrease in density, and hence of available phase space of protons. The corresponding low-momentum (k<kFp) 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).

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.

A unique advantage of SCGF simulations is the access to both off-shell, 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 quasi-particle 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 Ukτ=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⁻³. 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 Ukτ 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 Ukτ 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–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 momentum, the results of the 414 MeV and 450 MeV force go to zero relatively quickly, whereas the potential for the 500 MeV interaction shows a characteristic shallow minimum around k = 800 MeV.

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 low-momentum 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, neutron-rich 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 Ukτ 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=kF should be zero, thus signaling an infinite lifetime of quasi-particles 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]. For neutrons, this hole component is relatively independent of the isospin asymmetry. In fact Γk=0n increases from values around 30 − 40 MeV in symmetric nuclear matter to 40 − 50 MeV around η = 0.6, and decreases again to values around 20 − 30 MeV in neutron matter. In contrast, for protons the hole component decreases steadily with asymmetry. This reflects the decrease in phase space of protons as asymmetry increases. For both neutrons and protons, the cutoff dependence of the hole component is of the order of 10 MeV and relatively independent of isospin asymmetry.

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=kFτ. 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>kFτ) of Γkτ is remarkably more cutoff dependent. The momentum dependence immediately above kFτ and up to about 500 − 600 MeV is very similar in all cases, and shows an increase from zero up to an isospin-dependent 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<kFτ. Single-particle 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].

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–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 in-medium 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].