Modeling Neutron Star Matter in the Age of Multimessenger Astrophysics

The interpretation of the available and forthcoming data obtained from multimessenger astrophysical observations -- potentially providing unprecedented access to neutron star properties -- will require the development of novel, accurate theoretical models of dense matter. Of great importance, in this context, will be the capability to devise a description of thermal effects applicable to the study of quantities other than the equation of state, such as the transport coefficients and the neutrino mean free path in the nuclear medium. The formalism based on correlated basis states and the cluster expansion technique has been previously employed to derive a well-behaved effective interaction -- suitable for use in standard perturbation theory -- from a state-of-the-art nuclear Hamiltonian, including phenomenological two- and three-nucleon potentials. Here, we provide a comprehensive and self-contained account of the extension of this approach to the treatment of finite-temperature effects, and report the results of numerical calculations of a number of properties of nuclear matter with arbitrary neutron excess and temperature up to 50 MeV.


INTRODUCTION
The first detection of a gravitational wave signal consistent with emission from a coalescing binary neutron-star system (Abbott et al. 2017a), supplemented by the later observation of electromagnetic radiation by space-and ground-based telescopes (Abbott et al. 2017b), arguably opened up a new age for both astrophysics and nuclear physics research.
In years to come, multimessenger observations are expected to provide unprecedented information on neutron star structure and dynamics, which will allow to shed light not only on bulk properties of nuclear matter-such as the Equation of State (EOS) determining the neutron-star mass and radius (see, e.g., Bauswein 2019)-but also on the underlying nuclear dynamics. The possibility to exploit the available data to infer direct information on nucleon interactions at microscopic level has been recently analysed in the pioneering study of Maselli et al. (2021).
To meet the challenges posed by the interpretation of upcoming data, theoretical models must be capable to provide a consistent description of both equilibrium and dynamical properties of neutron-star matter in the temperature regime corresponding to T m π -m π ≈ 140 MeV being the pion mass-where nucleons are believed to be the dominant degrees of freedom. In addition to the EOS, these include the transport coefficients describing matter viscosity (Benhar & Valli 2007;Camelio et al. 2022a,b) and heat transfer (Benhar et al. 2010), as well as the neutrino mean free path in the nuclear medium (Lovato et al. 2013(Lovato et al. , 2014, which play a critical role in the evolution of proto neutron stars (Camelio et al. 2017), as well as in the post-merger phase of neutron-star coalescence (Figura et al. 2020(Figura et al. , 2021. Owing to the complexity and non perturbative nature of nuclear forces, however, the achievement of the above goal involves non trivial conceptual and computational issues. The equation of state of nuclear matter can be obtained from highly accurate ab initio calculations, performed using phenomenological Hamiltonians-strongly constrained by the observed properties of two-and three-nucleon systems-and advanced theoretical approaches for the solution of the quantum-mechanical many-body problem (for a recent review, see, e.g., Benhar & Fantoni 2020). On the other hand, the present development of computational techniques does not allow to accurately describe transport phenomena or neutrino reaction rates using the same Hamiltonians.
In nuclear many-body theory, the problem of the occurrence of non perturbative interactions is circumvented through a renormalisation of the bare Hamiltonian, leading to its replacement with a density-dependent effective Hamiltonian suitable for use in perturbation theory. This scheme has been followed to derive effective interactions within the G-matrix approach (for a review, see Baldo 1990) or using the formalism of Correlated Basis Functions and the cluster expansion technique (Clark 1979;Fantoni & Pandharipande 1988). The resulting potentials are well behaved, and have been used to perform calculations of a variety of nuclear matter properties relevant to astrophysical processes (see Baldo et al. 1997;Akmal et al. 1998). More recently, it has been suggested that effective interactions suitable for perturbative calculations can also be obtained combining potentials derived within chiral effective field theory (χEFT) and renormalisation group evolution to low momentum, (see, e.g., Drischler et al. 2016;Keller et al. 2021;Drischler et al. 2021). As clearly stated by Tews et al. (2018), χEFT appears to be inherently inadequate for applications in the high-density region relevant to neutron stars; see also (Benhar 2021). However, studies of the convergence of the chiral expansion up to about twice nuclear-saturation density provide a quantitative estimate of the theoretical uncertainty associated with the derivation of the nuclear Hamiltonian.
It has to be pointed out that the effective interactions obtained from realistic dynamical models are conceptually different from those designed to merely explain the empirical information on average properties of atomic nuclei and isospin-symmetric nuclear matter near equilibrium density, such as the Skyrme interactions employed, e.g., in the studies of Chabanat et al. (1997) and Rikovska Stone et al. (2003). Because they lack a connection with microscopic dynamics, these models are in fact unable to describe nucleon-nucleon scattering in the nuclear medium, and their applications to the study of neutron star properties are severely limited (Benhar et al. 2010).
Starting from the early 2000s, the Correlated Basis Function, or CBF, formalism has been exploited to derive effective interactions from phenomenological Hamiltonians comprising both two-and three-nucleon potentials (Benhar & Valli 2007;Lovato et al. 2011Lovato et al. , 2013. The capability to embody into a density-dependent NN potential the effects of irreducible three-nucleon interactions-which are known to become large, or even dominant, at high densities-is a critical feature of the CBF approach. The effective interaction obtained from the latest and most advanced implementation of this approach, thoroughly described by Lovato et al. (2013), has been applied to the determination of a variety of properties of nuclear matter, at both zero and nonzero temperature, within a unified scheme (Benhar & Lovato 2017;Camelio et al. 2017).
This article is aimed at providing a comprehensive and self-contained account of the extension of the work of Benhar & Lovato (2017) to the case of non vanishing temperature, including a detailed derivation of the procedure employed to achieve thermodynamic consistency. The conceptual analogy between the CBF effective interaction, obtained from renormalisation in coordinate space, and the low-momentum interactions obtained from renormalisation group evolution is also discussed.
The body of the paper is organised as follows. The dynamical model underlying nuclear many-body theory and the derivation of the microscopic effective interaction using the CBF formalism and the cluster expansion technique are outlined in Section 2. Section 3 is devoted to a discussion of the issues associated with the use of many-body perturbation theory at nonzero temperature, while the results of numerical calculations of selected properties of hot nuclear matter-including the EOS, the single nucleon spectra, and the effective masses-and neutron stars are reported in Sections 4 and 5, respectively. Finally, our main findings and the prospects for future developments of our work are summarised in Section 6.

NUCLEAR DYNAMICS
Ideally, a model of nuclear dynamics should be capable to provide a unified description of all nucleon systems, from the deuteron to neutron stars (Wiringa 1993;Benhar & Fantoni 2020). In this section, we give an overview of the prominent issues associated with the treatment of nuclear interactions at supranuclear densities, and outline the key elements of the approach employed in our work.

The nuclear Hamiltonian
Nuclear many-body theory is founded on the hypothesis that nuclear systems can be treated as collections of pointlike protons and neutrons, whose dynamics are dictated by a non relativistic Hamiltonian of the form 1 In the above equation, A is the number of nucleons, p i and m denote the momentum of the i-th nucleon and its mass, and the potentials v ij and V ijk describe two-and three-nucleon interactions, respectively. The nucleon-nucleon (NN) potential is designed to reproduce the measured properties of the two-nucleon system, in both bound and scattering states, and reduces to the Yukawa one-pion exchange (OPE) potential at large distances. It can be conveniently written in the form where the functions v p only depend on the distance between the interacting particles, r ij = |r i − r j |, while the operators O p ij account for the strong spin-isospin dependence of nuclear forces, as well as for the occurrence of noncentral interactions. The most important contributions to the sum appearing in the right-hand side of Eq. (2) are those associated with the six operators where σ i and τ i are Pauli matrices acting in spin and isospin space, respectively, and the tensor operator S ij is defined as Note that the OPE potential can be written in terms of the O p≤6 ij defined by Eqs. (3) and (4). State-of-the-art phenomenological models of v ij , such as the Argonne v 18 (AV18) potential )-providing an accurate fit of the NN scattering phase shifts, the low-energy NN scattering parameters, and deuteron properties-include twelve additional terms. The operators corresponding to p = 7, . . . , 14 are associated with non-static components of the potential, notably spin-orbit terms, while those corresponding to p = 15, . . . , 18 take into account small violations of charge symmetry and charge independence.
The results discussed in this article have been obtained using the Argonne v 6 (AV6P) potential, constructed projecting the full AV18 on the basis of the six operators of Eqs. (3)-(4) (Wiringa & Pieper 2002). This potential reproduces the deuteron binding energy and electric quadrupole moment with accuracy of 1%, and 4%, respectively, and provides an excellent fit of the phase shifts in the 1 S 0 channel 2 .
As an example, the energy dependence of the 1 S 0 and 1 P 1 phase shifts is illustrated in Fig. 1. It is apparent that the results obtained using the AV6P potential, represented by the solid lines, provide an accurate description of the data reported by Bergervoet et al. (1990), Stoks et al. (1993), and Arndt et al. (2007) up to beam energies ∼ 600 MeV, well beyond the pion production threshold, E thr ≈ 280 MeV. For comparison, the 1 S 0 phase shifts calculated using the full AV18 potential are also shown by the dashed line of panel (A).
It is very important to realise that the capacity to explain scattering data at large energy is critical to assess the ability of a potential model to describe the properties of nuclear matter in the high-density region, relevant to neutron star physics. To see this, consider a scattering process involving two nucleons embedded in the nuclear medium at baryon number density B . Owing to Pauli's exclusion principle, in the strongly degenerate regime typical of neutron stars collisions predominantly involve particles with energies close to the Fermi energy, E F . As a consequence, a relation can easily be established between the kinetic energy of the beam particle in the laboratory frame, E lab , and the Fermi energy, which in turn is simply related to the density. In the case of head-on collisions, the resulting expression reads where the degeneracy of the momentum eigenstates, ν, equals to 2 and 4 in pure neutron matter (PNM) and isospin-symmetric matter (SNM), respectively. The densities of PNM corresponding to beam energies in the range 0 ≤ E lab ≤ 630 MeV, expressed in units of the equilibrium density of SNM, 0 = 0.16 fm −3 , are given in the top axes of Fig. 1. Three-body forces are long known to be required to model the interactions of extended composite bodies, such as protons and neutrons, without considering their internal structure explicitly (see, e.g., Friar 1986). In nuclear manybody theory, the inclusion of irreducible three-nucleon (NNN) forces, described by the potential V ijk , is needed to explain both the observed binding energies of the three-nucleon systems and saturation-that is, the occurrence of a minimum of the energy per particle at non-vanishing density 0 -of SNM.
The nature of NNN interactions has been first highlighted in the seminal paper of Fujita & Miyazawa (1957). These authors argued that the most prominent mechanism is the two-pion exchange process, in which one of the nucleons participating in a NN interaction is excited to a ∆ resonance, that then decays in the aftermath of the interaction with a third nucleon. Commonly used phenomenological models of the NNN force, such as the Urbana IX (UIX) potential adopted in this work (Pudliner et al. 1995), are written in the form where V 2π ijk is the attractive Fujita-Miyazawa term, while V N ijk is a purely phenomenological repulsive term. Like most phenomenological potentials, the UIX model involves two parameters. The strength of the two-pion exchange contribution, V 2π ijk , is fixed in such a way as to reproduce the observed ground-state energies of 3 He and 4 He, while that of the isoscalar repulsive contribution, V R ijk , is adjusted to obtain the saturation density of SNM inferred from extrapolation of nuclear data. It is remarkable, however, that saturation of SNM at density B ∼ 0 is also predicted by a potential tuned only to reproduce the properties of the few-nucleon systems (Coon & Ha 2001;Lovato et al. 2012). Nuclear Hamiltonians constructed combining the AV18 NN potential and a phenomenological NNN potential, such as the UIX model, have been shown to possess a remarkable predictive power. The results of Quantum Monte Carlo (QMC) calculations, extensively reviewed by Carlson et al. (2015), reproduce the measured energies of the ground and low-lying excited states of nuclei with mass number A ≤ 12 to few percent accuracy. The results reported in the present work have been obtained using the AV6P+UIX Hamiltonian.
Over the past two decades, a great deal of effort has been devoted to the derivation of nuclear potentials within the framework of χEFT (see, e.g., Epelbaum et al. 2009;Machleidt & Entem 2011, and references therein). This formalism, originally proposed by Weinberg (1990), is based on the use of effective Lagrangians involving pions and low-momentum nucleons, constrained by the broken chiral symmetry of strong interactions.
The approach based on χEFT provides an elegant and systematic scheme, in which the nuclear interaction is expanded in powers of a small parameter, e.g. the ratio between the pion mass or the typical nucleon momentum and the scale of chiral symmetry breaking, Λ χ ∼ 0.8 -1 GeV. In addition, it allows to obtain two-, three-, and many-nucleon potentials in a fully consistent fashion. Local coordinate-space representations of potentials derived from χEFT, suitable to carry out nuclear matter calculations using advanced many-body techniques, have been first derived by Gezerlis et al. (2013Gezerlis et al. ( , 2014. More general coordinate-space potential models, in which the appearance of ∆ resonances in intermediate states is explicitly taken into account, have been developed by Piarulli et al. (2015).
In principle, nuclear Hamiltonians obtained from χEFT may be used to carry out calculations of nuclear matter properties relevant to neutron stars. However, it has to be kept in mind that, being based on a low-momentum expansion, χEFT is inherently limited to the description of matter at densities 2 0 . This problem clearly emerges from the phase-shift analysis of Gezerlis et al. (2013), showing that their next-to-next-to-leading order (N 2 LO) potential only describes the data at E lab 150 MeV.
The potentials of Piarulli et al. (2015) have been obtained by fitting NN scattering data up to E lab = 125 or 200 MeV, and splitting the strong interaction component into short-and long-range contributions involving different coordinatespace cutoffs, R L and R T . The results of this analysis show that a suitable combination of R L and R T allows to obtain a description of scattering data comparable to that provided by the AV18 potential. It should be pointed out, however, that this procedure, while being fully justified on phenomenological grounds, involves a departure from the original formulation of the approach based on χEFT.

Renormalisation in coordinate space: the CBF effective interaction
The observation that the central density of atomic nuclei obtained from elastic electron-nucleus scattering data (Frois & Papanicolas 1987) becomes largely A-independent for A 12, its value being ∼ 0 , indicates that NN forces are strongly repulsive at short distance. This feature has been qualitatively confirmed by the results of pioneering lattice calculations based on the fundamental theory of strong interactions: Quantum Chromo-Dynamics, or QCD (Ishii et al. 2007).
Owing to the presence of the repulsive core, the matrix elements of the NN potential between eigenstates of the non interacting system are large, and standard many-body perturbation theory cannot be used to carry out calculations of nuclear properties.
A time-honored theoretical approach to overcome the above problem-known as G-matrix perturbation theory-is based on the replacement of the bare NN potential with a well-behaved operator describing NN scattering in the nuclear medium. The lowest order approximation of the resulting expansion has been extensively employed in early studies of cold nuclear matter (see, e.g., Day 1967;Bethe 1971). More recent developments, allowing the inclusion of higher order terms and the treatment of matter at finite temperature, are thoroughly reviewed in the volume edited by M. Baldo (1990).

The CBF effective interaction, that can be written in the form
with the O p ij given by Eqs. (3) and (4), is defined by the equation Here, |Φ 0 and T F denote the ground state of the non interacting Fermi gas at density B and the corresponding energy, respectively, while H is the nuclear Hamiltonian of Eq.
(1). The correlated ground state, |Ψ 0 , is obtained from the corresponding Fermi gas state |Φ 0 through the transformation where the operator F is a product of two-body correlation operators, whose structure is chosen in such a way as to reflect the complexity of NN interactions. The resulting expression, to be compared to Eq.
(2), reads with where the f p are NN correlation functions, to be discussed below. Note that the inclusion of the operator S, appearing in the right-hand side of Eq. (10), is needed to symmetrise the product under particle exchange, because, in general, The determination of the effective interaction from Eq. (8) is based on the cluster expansion of the left-hand side, leading to where (∆E) n denotes the contribution to the Hamiltonian expectation value arising from n-nucleon clusters. In order to explicitly take into account three-nucleon forces, the effective interaction employed in this work has been derived including terms with n = 2 and 3. The radial dependence of the correlation functions is obtained solving a set of Euler-Lagrange equations derived from functional minimisation of the two-body cluster approximation to H (Pandharipande & Wiringa 1979). The range of the f p (r ij ) is fixed in such a way as to simultaneously reproduce the ground-state energies of PNM and SNM obtained from highly accurate many-body calculations, carried out using the Auxiliary-Field-Diffusion-Monte-Carlo (AFDMC) technique or the variational approach referred to as Fermi-Hyper-Netted-Chain/Single-Operator-Chain (FHNC/SOC) (Benhar & Lovato 2017).
Note that the CBF effective interaction depends on density through both dynamical correlations, described by the operator F ij , and statistical correlations, arising from the antisymmetric nature of the state |Φ 0 . The radial dependence of v eff in the S = 0 and T = 1 channel at baryon density B = 0.04, 0.32 and 0.48 fm −3 is displayed in Fig. 2.
The short-distance behaviour of the f p -shaped by the strongly repulsive core of the NN potential-brings about a strong suppression of the probability of finding two nucleons at relative distance r 1 fm. To see this, consider a nucleon pair in the state of total spin and isospin S = 0 and T = 1, angular momentum = 0 and momentum k, embedded in in nuclear matter at equilibrium density. Its relative motion is described by the wave function  A great deal of effort has been also devoted to the derivation of soft effective interactions-suitable to carry out perturbative calculations of nuclear matter in the Fermi gas basis-from renormalisation-group (RG) evolution of potentials obtained within χEFT (see, e.g., Bogner et al. 2010;Furnstahl 2012). It has to be emphasised that the RG approach is conceptually equivalent to the one based on the CBF formalism. Within the RG scheme, screening of the repulsive core of the NN interaction is obtained by integrating out high-momentum components, and the evolution is driven by the value of the momentum cutoff. Using the CBF formalism, on the other hand, the same effect is realised in coordinate space through the action of the correlation functions, as illustrated in Fig. 3, and evolution is driven by nuclear matter density, determining the average distance between matter constituents. However, it has to be kept in mind that-in view of the limitations discussed in Section 2.1-the RG approach, being based on potentials derived from χEFT, is expected to be applicable in a narrow density region.

FINITE-TEMPERATURE PERTURBATION THEORY
The basic assumption underlying our treatment of nuclear matter at T = 0 is that at low-to-moderate temperatures-typically T m π , m π ≈ 140 MeV being the pion mass-nuclear dynamics, described by the potentials appearing in the Hamiltonian of Eq. (1), is largely unaffected by thermal effects. In principle, the CBF effective interaction involves an additional temperature dependence associated with the correlation functions, since the Fermi distribution appears in the Euler-Lagrange equations determining their shape. However, the results of detailed numerical calculations have shown that thermal modifications of the f p of Eq. (11) turn out to be negligibly small up to T ∼ 50 MeV (Valli 2007). The results reported in this article have been obtained using the zero-temperature effective interaction, involving correlation functions computed at T = 0.

Thermodynamics
All thermodynamic functions of a system in thermal equilibrium at temperature T can be derived from the grand canonical potential, defined as where β = 1/T . In the case of a one-component system, the partition function, Z, is given by In the above equation, µ is the chemical potential, while H and N denote the Hamiltonian and the particle number operator, respectively 3 . The basis for the derivation of finite-temperature perturbation theory is provided by the Bloch equation (see, e.g., Thouless 1961) to be solved with the boundary condition Φ(0) = 1. The perturbative expansion of the grand canonical partition function is easily obtained exploiting the formal similarity between Eq. (17) and the time-dependent Schrödinger equation of quantum mechanics, and rewriting the Hamiltonian in the form Substitution of Eq. (18) into the right-hand side of the Bloch equation, leading to shows that the formalism of time-dependent perturbation theory can be readily generalised by replacing t → −iβ, and using the operator H 0 to define the appropriate interaction picture. The fundamental relation (see, e.g., Landau & Lifshitz 1969) provides a link between the grand canonical potential, the pressure P , and the free energy F = E − T S, with E and S being the energy and entropy of the system, respectively. From Eq. (20) if follows that In the following, we will discuss the application of the above results to a system described by the Hamiltonian (18), with where, in general and Here, the label k specifies both the particle momentum and the discrete quantum numbers corresponding to oneparticle states, a † k and a k denote creation and annihilation operators, respectively, and v is the potential describing interparticle interactions. The single-particle potential U k , which in principle does not affect the results of calculations of physical quantities, is chosen in such a way as to improve the convergence of the perturbative expansion, or to fulfill specific conditions (see, e.g., Baldo 1990) .
Note that, through the use of the CBF effective interaction discussed in Section 2.2, the formalism based on the Hamiltonian defined by Eqs. (18) and (22)-(24) allows to take into account two-and three-nucleon interactions in a consistent fashion.
It has to be pointed out that, according to Eq. (20), the pressure can be written in the form with B = N/V , implying that at equilibrium, that is, for P = 0, µ = F/N . This result can be seen as the generalisation of the Hugenholtz-Van Hove theorem (Hugenholtz & Van Hove 1958) to the case of non vanishing temperature.

Perturbation Theory
At first order in H I , the grand canonical potential is given by (see, e.g., Lejeune et al. 1986) with where |kk A = |kk − |k k denotes an antisymmetrised two-particle state, and n k is the Fermi distribution, given by From Eqs. (27) and (28) it follows that the free energy per particle can be cast in the form In principle, for any assigned values of temperature and chemical potential, the above equations provide a scheme for the determination of the equation of state of nuclear matter at finite temperature, P = P (µ, T ). In view of the fact that baryon number is conserved by all known interactions, however, in nuclear matter it is convenient to use baryon density as an independent variable, and determine the chemical potential from the relation In the T → 0 limit and the chemical potential is given by µ = e k F , with the Fermi momentum being defined as k F = 6π 2 B /ν 1/3 .
For T = 0 and density-dependent potentials, thermodynamic consistency is not trivially achieved in perturbative calculations. A clear manifestation of this difficulty is the mismatch between the value of pressure obtained from Eq. (25) and the one resulting from the alternative-although in principle equivalent-thermodynamic expression A procedure fulfilling the requirement of thermodynamic consistency by construction can be derived from a variational approach, based on minimisation of the trial grand canonical potential (Heyer et al. 1988) with respect to the form of n k . Note that the above expression-the use of which is fully legitimate in the variational context-can also be obtained in first order perturbation theory neglecting terms involving ∂Ω 1 /∂T and ∂Ω 1 /∂µ (Lejeune et al. 1986). The condition δ Ω δn k = 0 , turns out to be satisfied by the distribution function with and Within this scheme, that reduces to the standard Hartee-Fock approximation in the case of density-independent potentials, all thermodynamic functions at given temperature and baryon density can be consistently obtained using the distribution n k of Eq. (36). Note, however, that, because both U k and δe depend on n k , see Eqs. (37) and (38), calculations must be carried out self-consistently, applying an iterative procedure.

EQUATION OF STATE OF HOT NUCLEAR MATTER
We will consider nuclear matter at fixed baryon density where the index λ = 1, 2, 3, 4 labels spin-up protons, spin-down protons, spin-up neutrons and spin-down neutrons, respectively, and the corresponding densities are denoted λ = x λ B . Even though our formalism is completely general, in the following we will restrict ourselves to the case of unpolarised matter, in which x 1 = x 2 and x 3 = x 4 . 4.1. Free energy per nucleon and single particle spectrum Following the procedure described in Sections 3.1 and 3.2, the energy per nucleon at temperature T can be obtained from In the above equation with the functions v eff,p and the operators O p defined in Eqs. (2) and (3), respectively, whereas |λλ denotes a two-nucleon state in spin-isospin space. The temperature dependence of the interaction contributions is contained in the generalised Slater functions L λ (x), defined as with the Fermi distribution given by where and j 0 (z) = sin z/z. In the T → 0 limit, µ λ = e k F λ , with k F λ = (6π 2 ρ λ ) 1/3 , and L λ (x) reduces to the ordinary Slater function (k F λ x), with (z) = 3(sin z − z cos z)/z 3 . At T = 0, however, E/N and e k,λ must be computed self consistently, and the chemical potentials µ λ are determined by the conditions The free energy per nucleon  (48), with the CBF effective interaction described in Section 2.2.
is obtained combining the above results with the corresponding expression of the entropy per nucleon Figure 4 shows the density and temperature dependence of the free energy per nucleon of SNM and PNM, corresponding to proton fraction Y p = x 1 + x 2 = 0.5 and 0, respectively, obtained from the procedure described above using the CBF effective interaction.
The formalism employed in this work allows to carry out calculations of the properties of isospin-asymmetric matter for any values of the proton fraction Y p . The smooth Y p -dependence of the free energy per nucleon is illustrated in Fig. 5, showing results corresponding to Y p = 0.05 and 0.25.
The energy per nucleon of asymmetric matter at zero temperature is often derived from an expansion in powers of the neutron excess δ = 1 − 2Y p . The resulting expression with the symmetry energy being given by E sym ( B ) = [E( B , 1) − E( B , 0)]/N , has been shown to provide a remarkably accurate approximation to the results of calculations in which the dependence on Y p is taken into account at microscopic level (see, e.g., Benhar & Lovato 2017). To test the validity of the generalisation of Eq. (49) to finite temperatures-previously discussed by, e.g., Camelio (2018) and Li et al. (2021)-we have compared the Y p dependence of the free energy per nucleon at different values of T to the predictions of the quadratic interpolation formula with The results shown in Fig. 6, corresponding to a representative baryon density B = 0.32 fm −1 , indicate that the excellent agreement observed at T = 0 tends to slowly deteriorate as the temperature increases. However, at T as high as 50 MeV the deviations-which are vanishing by construction at Y p = 0 and 0.5-turn out to be still limited to ∼ 20% at Y p 0.1.

Pressure and thermodynamic consistency
As pointed out in Sec. 3.2, in the case of density-dependent effective interactions thermodynamic consistency, expressed by the condition where µ n and µ p are the neutron and proton chemical potentials, requires the inclusion of a correction to the Hartree-Fock single-particle spectrum, given by the second line in the right-hand side of Eq. (45). The thermodynamic consistency of our formalism is illustrated in Fig. 7, showing the pressure of PNM at T = 0 and 50 MeV. Solid lines and triangles correspond to results obtained from the left-and right-hand side of Eq. (52), respectively. For comparison, the pressure of PNM at T = 50 MeV computed omitting the self-consistency correction to the single-particle spectrum is shown by the dashed line. It should also be pointed out that the effect of this correction on pressure turns out to be nearly independent of temperature.
A comparison between the curves corresponding to T= 0 and 50 MeV illustrates the density dependence of thermal effects. It is apparent that, while providing large to significant contributions at densities 2 0 , the thermal pressure plays a negligible role at higher densities, where statistical and dynamical effects dominate.

Single-nucleon properties
In interacting many-body systems, single-particle states are in general not uniquely defined. Owing to translation invariance, however, in nuclear matter they can be identified using momentum eigenvalues 4 , and their properties can be derived within the framework of Landau's theory of normal Fermi liquids (Baym & Pethick 1991). As pointed out by Benhar & Lovato (2017) for the case of vanishing temperature, the single-nucleon energy of Eq. (45) can be interpreted as the energy of a quasiparticle carrying momentum k and spin-isospin quantum numbers specified by the index λ.
The chemical potentials µ λ are obtained from Eq. (46). At T = 0, they depend on the spectra e k,λ through the Fermi distributions, while for T = 0 they reduce to the Fermi energies e k F λ , as dictated by the generalisation of Hugenholtz-Van Hove theorem to the case of a multicomponent system. As an example, in the left panel of Fig. 8 the neutron and proton chemical potentials of nuclear matter at proton fraction Y p = 0.1-denoted µ n and µ p , respectively-are displayed as a function of baryon density and temperature.
The description of quasiparticle dynamics is largely based on the effective mass m λ , defined by the equation The effective mass, dictating the nucleon dispersion relation in matter, plays a critical role in determining the rates of a number of processes relevant to neutron star properties, such as, e.g., neutrino emission and absorption. The right panel of Fig. 8 illustrates the density and temperature dependence of the proton and neutron effective masses in nuclear matter at proton fraction Y p = 0.1, expressed in units of the nucleon mass m N . The results show that the temperature dependence, often neglected in astrophysical applications, is, in fact, significant. At B = 0 (3 0 ) the values of m p and m n at T = 50 Mev turn out to differ from the corresponding zero-temperature values by ∼ 20% and 10% (15% and 8%), respectively.

NEUTRON STAR MATTER
The calculation of neutron star properties requires the determination of the EOS of charge neutral matter consisting of neutrons, protons and leptons in equilibrium with respect to the neutron β-decay and lepton capture processes where denotes the flavour of the lepton participating in the reactions. For T = 0, the results of calculations of the EOS, described in the next section, provide the input needed to solve the Tolman-Oppenheimer-Volkoff (TOV) equations (Tolman, R. C. 1939;Oppenheimer, J. R. and Volkoff, G. M. 1939), determining mass and radius of cold neutron stars.

Equation of state of β-stable matter
Let us consider matter comprising electrons and muons, hereafter referred to as npeµ matter. Under the assumption of transparency to neutrinos and antineutrinos, implying that these particles can freely escape and have vanishing chemical potentials, the equilibrium condition reduces to where µ denotes the chemical potential of leptons of flavour . The above condition must be fulfilled together with the additional constraint of charge neutrality, requiring that with Y = / B , being the lepton density. For any given values of B and T , Eqs. (55) and (56) univocally determine the proton and lepton fractions, needed to obtain the EOS of β-stable matter. Figure 9 shows the density and temperature dependence of the proton fraction of charge neutral npeµ matter in β equilibrium. The most prominent thermal effect turns out to be a departure from the monotonically increasing behavior of Y p as a function of density, leading to the appearance of pronounced minimum for T 30 MeV.
In addition to the baryon contributions discussed in the previous sections, the free energy and pressure of β-stable npeµ matter comprise contributions arising from the presence of leptons, which can be safely treated as non interacting relativistic fermions. The solid and dashed lines of Fig. 10 illustrate the density dependence of the total pressure of npeµ matter where the indices B and L label the baryon and lepton components, for T = 0 and 50 MeV, respectively. A comparison between the zero-temperature total and lepton pressure, displayed by the dot-dash line, shows that the baryon contribution is largely dominant at all densities.

Neutron star properties
The theoretical framework described in the previous sections has been employed to evaluate a variety of properties of cold neutron stars (Sabatucci & Benhar 2020;Maselli et al. 2021;Tonetto et al. 2021). The results of these analyses show that the predictions of our approach are compatible with the determinations of maximum mass, tidal deformability and radius based on recent multimessenger observations (Antoniadis et al. 2013;Cromartie et al. 2019;Abbott et al. 2017a;Riley et al. 2019). A pioneering application of the finite-temperature EOSs to study the evolution of proto-neutron stars has been also carried out by Camelio et al. (2017).
The determination of the mass-radius relation of hot neutron stars involves non trivial additional difficulties, because hydrodynamic equilibrium and heat transfer-determining the temperature profile in the star interior-must be consistently taken into account. In order to provide an admittedly rough estimate of the impact of thermal effects on the maximum neutron star mass, in Fig. 11 we show results obtained by solving the TOV equations with EOSs evaluated at constant temperature using our approach.
In the region B / 0 > 4, the EOSs computed following the procedure discussed in Sect. 5.1 with 0 ≤ T ≤ 50 MeV have been smoothly matched to the zero-temperature EOS obtained by Akmal et al. (1998) using the AV18+UIX nuclear Hamiltonian, which is used as baseline for the determination of the CBF effective interaction at 0.25 ≤ B / 0 ≤ 4; see Section 2. The validity of this prescription is supported by the observation that the contribution of thermal pressure becomes vanishingly small at ∼ 4 0 ; see Fig. 10.
The crust region has been described using the model adopted by Akmal et al. (1998), which provides a description of matter at B < 0.1 fm −3 . The use of a zero-temperature EOS appears to be justified for the purpose of our exploratory analysis, because the value of the maximum mass is largely unaffected by the structure of matter in the low-density crust region.
The results of Fig. 11 show that the impact of thermal effects on the maximum mass M max , leading to its increase, is not large. The values of M max at T = 0 and 50 MeV turn out to be 2.37 M and 2.49 M , with M being the solar mass, respectively.
It has to be pointed out that the pattern emerging from our analysis-while being in qualitative agreement with the results of existing studies carried out within the framework of the relativistic mean field (RMF) approach or using Skyrme-type models (Prakash et al. 1997;Kaplan et al. 2014)-turns out to be at variance with the findings of studies based on G-matrix perturbation theory and phenomenological nuclear Hamiltonians, which predict a small decrease of the maximum mass at nonzero temperature (Figura et al. 2020(Figura et al. , 2021.

SUMMARY AND OUTLOOK
We have developed a consistent formalism suitable to carry out perturbative calculations of the properties of both cold and hot nuclear matter at arbitrary proton fraction. The distinctive feature of our method is the use of an effective nuclear Hamiltonian, based on a phenomenological potential providing an accurate description of nucleonnucleon scattering in the energy range relevant to dense matter in the neutron star interior.
Our results reproduce the zero-temperature EOSs of PNM and SNM obtained from highly refined many-body approaches, such as AFDMC and FHNC/SOC, by construction. Moreover, SNM calculations performed using the CBF effective interaction yield the correct equilibrium density, and an interaction energy within ∼ 15% of the empirical value (Benhar & Lovato 2017).
Recently, Lovato et al. (2022), have carried out extensive analysis of the EOS of zero-temperature PNM, aimed at comparing different nuclear Hamiltonians and computational techniques. The results of this study, showing that the AV6P+UIX Hamiltonian yields results in close agreement with those obtained from the full AV18+UIX up to B = 2 0 , strongly supports the use of this somewhat simplified model to describe nuclear dynamics at supranuclear density.
Within our computational scheme, based on finite-temperature perturbation theory, thermal effects on the energy per baryon and the single-nucleon spectrum are taken into account in a consistent fashion, without introducing any simplifying assumptions, and the thermodynamic consistency condition turns out to be fulfilled to very high accuracy.
The results reported in this work suggest that our approach has the potential to provide accurate evaluations of the variety of properties of dense nuclear matter needed for numerical simulations of astrophysical processes.
It is very important to note that the availability of a reliable description of neutron star matter in the regime in which nucleons are the relevant degrees of freedom is also essential to firmly establish the occurrence of transitions to more exotic forms of matter-involving baryons other than protons and neutrons as well as deconfined quarks-which are expected to become energetically favoured at higher densities.
The analyses of recent astrophysical data have provided information about the radius of a neutron star of mass M = 1.4 M , whose central density does not exceed ∼ 3 0 . The reported values-R = 12.45 ± 0.65 km (Miller et al. 2021), 12.33 +0.76 −0.81 km (Riley et al. 2019), and 12.18 +0.56 −0.79 km (Raaijmakers et al. 2021)-turn out to be compatible with the predictions of theoretical calculations performed using EOSs of purely nucleonic matter (Sabatucci & Benhar 2020), thus suggesting that in this case the paradigm of nuclear many-body theory can still be applied.
The equatorial radius of the neutron star J0740+6620, having mass M = 2.072 +0.067 −0.066 M and central density ∼ 4 0 , has been also evaluated to be R = 12.39 +1.30 −0.98 km (Riley et al. 2019). The small difference between the radii of stars of mass 1.4 and 2.071 M , implying that the EOS is still rather stiff at > 3 0 , appears to rule out the occurrence of a strong first order phase transition in the density range 3 0 B 4 0 . Independent empirical information about the limits of applicability of the description of dense nuclear matter in terms of nucleons is obtained from the analysis of the large body of electron-nucleus scattering data. Scaling in the variable y-the definition and physical interpretation of which is thoroughly discussed by West (1975) and Day et al. (1990)-has been observed by experiments using a broad range of targets, extending from 2 H (Ciofi degli Atti et al. 1987) to nuclei as heavy as 197 Au (Day, et al. 1987;Arrington et al. 1999). The results of these studies unambiguously demonstrate that in the kinematical region corresponding to momentum transfer q 1 GeV and large negative y the beam particles primarily couple to nucleons belonging to correlated pairs, the momentum of which can be in excess of 500 MeV. Based on the results of experimental analyses carried out at Jefferson Lab, Subedi et al. (2008), argued that the occurrence of strong nucleon-nucleon correlations is associated to large fluctuations of the nuclear density, that can locally reach values as high as ∼ 5 0 . According to the argument proposed in Section 2.1, this estimate is consistent with the observation of scaling at y ∼ −600 MeV, which in turn implies that at B ∼ 5 0 nuclear matter largely behaves as a collection of nucleons.
As a final remark, it has to be pointed out that the description of nuclear matter based on non relativistic nuclear many-body theory, providing the conceptual framework underlying our approach, unavoidably leads to violation of causality in the high-density limit. It has to be pointed out, however, that within the approach discussed in this paper the speed of sound in cold β-stable matter becomes larger than the speed of light at > 4.5 0 . Therefore, our EOS turns out to be relativistically consistent in the region relevant to the stability of a 2.2 M neutron star; see Fig. 11. The difficulties arising from violation of causality can be further alleviated taking into account relativistic boost corrections to the nucleon-nucleon potential, along the line discussed by Akmal et al. (1998). Work is presently being carried out to include these correction in our formalism.