© The Authors 2026

1. Introduction

The structure and composition of neutron stars (NSs) depend on the unknown properties of dense and strongly interacting baryonic matter. The same is true for the evolution of core-collapse supernovae (CCSNe) (Janka et al. 2007; Mezzacappa et al. 2015; Schneider et al. 2017; O’Connor & Couch 2018; Burrows et al. 2020), proto-neutron stars (PNSs) (Pons et al. 1999; Pascal et al. 2022), binary neutron star (BNS) mergers (Shibata & Taniguchi 2011; Rosswog 2015; Baiotti & Rezzolla 2017; Endrizzi et al. 2018; Ruiz et al. 2020; Prakash et al. 2021; Most et al. 2023), and the formation of black holes (BHs) in failed CCSNe (Sumiyoshi et al. 2007; Fischer et al. 2009; O’Connor & Ott 2011; Hempel et al. 2012).

The tremendous progress made by multimessenger astronomy of NSs over the past decade has contributed unprecedented and valuable knowledge about the properties of states of matter that are impossible to produce and study in terrestrial laboratories. The high-density behavior of the NS equation of state (EOS) was bracketed by measurements of pulsars with masses around or larger than 2 M_⊙ (Demorest et al. 2010; Antoniadis et al. 2013; Arzoumanian et al. 2018; Cromartie et al. 2020; Fonseca et al. 2021), and by the interpretation of NS coalescence in the GW170817 event (Abbott et al. 2017) as the collapse of a hypermassive star into a BH, which can be translated into an upper limit for the maximum mass that NSs can sustain (Margalit & Metzger 2017; Rezzolla et al. 2018; Khadkikar et al. 2021). The measurement of the combined tidal deformability of NSs with masses 1.17 ≲ M/M_⊙ ≲ 1.60 in the GW170817 event (Abbott et al. 2017, 2019) provides the first ever constraint on the behavior of neutron-rich matter over the density range 1 ≲ n/n_sat ≲ 3, where n_sat ≈ 0.16 fm⁻³ ≈ 2.7 × 10¹⁴ g/cm³ represents the nuclear saturation density. The determination of equatorial radii by NICER X-ray telescope, based on the analysis of pulse profiles of millisecond pulsars with masses 1.4 ≲ M/M_⊙ ≲ 2.1 (Riley et al. 2019, 2021; Miller et al. 2019, 2021; Vinciguerra et al. 2024; Choudhury et al. 2024; Mauviard et al. 2025), has provided additional information on the stiffness of neutron-rich EOSs in the suprasaturation domain.

In recent years, a large number of statistical inferences of NS EOS have been performed, which systematically address the role of the EOS model, the set of constraints, prior distributions, and other factors. For a discussion, see (Beznogov & Raduta 2023, 2024a,b) and references therein. Despite these efforts, the NS EOS remains largely unknown. Several factors contribute to this situation, including the insufficiently understood sensitivity of different astrophysical observations to various domains of density and isospin asymmetry (δ = (n_n − n_p)/n) of nuclear matter (NM), degeneracies related to the particle composition of dense matter, the unknown effective baryon-baryon interactions, and the still large uncertainties in astrophysical measurements. Here, n_n and n_p represent the neutron and proton particle densities, respectively, and n = n_n + n_p.

In the absence of constraints from nuclear physics experiments, the thermal response of dense matter, which governs the evolution of PNSs, CCSNe, BNS mergers, and stellar BH formation, is even more uncertain. Numerical simulations of CCSNe (Schneider et al. 2019; Yasin et al. 2020; Andersen et al. 2021), BNS mergers (Fields et al. 2023; Raithel & Paschalidis 2023), and stellar BH formation (Schneider et al. 2020) demonstrate that correlations exist between the evolution of these phenomena and the value of the nucleon effective mass, m_eff, at saturation density. Schneider et al. (2019), Yasin et al. (2020), Andersen et al. (2021) show that large values of m_eff favor high (low) values of the central density (temperature) in the cores of PNSs, as well as smaller PNS radii. By affecting the compactness of PNSs, m_eff also impacts their oscillations and, consequently, the peak frequency of gravitational waves (GWs; Andersen et al. 2021). Schneider et al. (2019), Andersen et al. (2021) also show that m_eff influences the temperature, proton fraction (Y_p = n_p/n), density, and radius of the neutrinosphere, as well as neutrino energies and luminosities. In failed CCSNe, the collapse into a BH occurs earlier for EOSs with higher values of m_eff. According to Fields et al. (2023), m_eff also affects the temperature and compactness of BNS mergers, as well as the strain and spectrum of GWs. In all these circumstances, the sensitivity to m_eff exceeds that to any other EOS parameter, including the stiffnesses of symmetric nuclear matter (SNM) and pure neutron matter (PNM).

The common feature of numerical simulations by Schneider et al. (2019, 2020), Yasin et al. (2020), Andersen et al. (2021), Fields et al. (2023) is the use of non-relativistic Skyrme interactions with a monotonic behavior of m_eff(n). The major advantage of Skyrme interactions is that analytical expressions are available for most thermodynamic and microscopic quantities (Constantinou et al. 2014, 2015), making it possible to understand the role of the effective mass in thermal pressure support, specific heat, and the evolution of various phenomena. Although these studies are informative, they do not guarantee that the correlations persist when EOSs based on more realistic interactions are employed. In particular, ab initio calculations with three-body forces predict that the nucleon effective mass as a function of density exhibits a U-shaped behavior (Baldo et al. 2014; Shang et al. 2020; Somasundaram et al. 2021). Depending on the forces and theoretical approaches used, for cold SNM the position of the minimum lies between n_sat (Somasundaram et al. 2021) and 4n_sat (Shang et al. 2020), with the value at the minimum being about 70% of the bare mass. For PNM, the minimum occurs around 2n_sat/3, while the value at the minimum is about 88% of the bare mass (Somasundaram et al. 2021). Complex density behaviors clearly make it difficult to establish any connection between the evolution of phenomena, in which wide ranges of densities are explored at different instances and spatial coordinates, and the value of m_eff(n_sat). Nevertheless, the value of m_eff at each density will influence the fate of these phenomena, implying that understanding the thermal behavior of NM governed by more realistic forces is extremely important.

This work aims to systematically study the thermal response of a set of Brussels extended Skyrme interactions generated within a Bayesian inference of the EOS of dense matter (Beznogov & Raduta 2024b). The advantage of employing Brussels extended Skyrme parametrizations is, among others, the ability to qualitatively reproduce the non-trivial density dependence of the nucleon effective mass in ab initio models (Baldo et al. 2014; Shang et al. 2020; Somasundaram et al. 2021). The use of general-purpose EOS tables based on such interactions in numerical simulations is expected to yield features distinct from those obtained when EOS tables based on simpler parametrizations are employed. In this way, we aim to contribute to a better understanding of the link between properties of NM and the evolution of CCSNe, PNSs, BNS mergers, and stellar BH formation. In the longer term, our results may provide constraints on the thermal behavior of NM.

The remainder of the paper is structured as follows. Section 2 provides a brief review of the theoretical framework. Section 3 discusses two sets of effective interactions and five particular models that exhibit extreme behaviors. Values of key NM parameters are also discussed in Sect. 3. The five models correspond to RB(BBSk1), RB(BBSk2), RB(BBSk3), RB(BBSk4), and RB(BBSk5) general-purpose EOS tables publicly available on the COMPOSE online repository (Typel et al. 2022)¹; they were introduced in Raduta & Beznogov (2025), hereafter referred to as Paper I. Selected properties of NSs based on these models are presented in Sect. 4. Sect. 5 is dedicated to the thermal response. Sect. 6 analyzes the stability of PNSs and remnants of BNS mergers in terms of the maximum gravitational and baryonic masses of isentropic stars. The conclusions are presented in Sect. 7.

In Sections 2, 3, and 5, we disregard the contributions of leptons and photons. Clusterization at densities below the saturation density and at temperatures below the critical temperature of Coulomb instabilities, discussed in Paper I, is also neglected. All results presented here correspond to homogeneous NM.

2. Formalism

Here, we treat hot, strongly interacting NM within the self-consistent Hartree-Fock approach (Negele & Vautherin 1972; Vautherin 1996), using Brussels extended Skyrme effective interactions (Chamel et al. 2009). In the absence of spin polarization and assuming zero electric charge, the energy density of bulk homogeneous matter is a sum of six terms,

$$\begin{array}{c}\hfill \mathcal{H}=k+h_0+h_3+h_{\mathrm{eff}}+h_4+h_5.\end{array}$$(1)

Here, k = ℏ²τ/2m is the kinetic energy term and 2/m = 1/m_n + 1/m_p, where m_i with i = n, p denotes the bare mass of nucleons; h₀ and h₃ are interaction terms that originate from the density-independent two-body term and the density-dependent term, respectively; h_eff, h₄, and h₅ are the momentum-dependent terms of the interaction. Each interaction term can be expressed analytically in terms of particle number densities, kinetic energy densities, and the parameters of the effective interaction (Ducoin et al. 2006; Beznogov & Raduta 2024b):

$$\begin{array}{cc}& h_0=C_0{n}^2+D_0{n}_3^2,\hfill \end{array}$$(2)

$$\begin{array}{cc}& h_3=C_3{n}^{\sigma +2}+D_3{n}^{\sigma }{n}_3^2,\hfill \end{array}$$(3)

$$\begin{array}{cc}& h_{\mathrm{eff}}=C_{\mathrm{eff}}n\tau +D_{\mathrm{eff}}{n}_3{\tau }_3,\hfill \end{array}$$(4)

$$\begin{array}{cc}& h_4=\frac{t_4}{16}[3n\tau -(2x_4+1)n_3{\tau }_3]n^{\beta }\hfill \end{array}$$(5)

$$\begin{array}{cc}& h_5=\frac{t_5}{16}[(4x_5+5)n\tau +(2x_5+1)n_3{\tau }_3]n^{\gamma }.\hfill \end{array}$$(6)

Here, n = n_n + n_p and n₃ = n_n − n_p denote the isoscalar and isovector particle number densities, and τ = τ_n + τ_p and τ₃ = τ_n − τ_p denote the isoscalar and isovector kinetic energy densities. The parameters C₀, D₀, C₃, D₃, C_eff, D_eff, x₄, t₄, x₅, t₅, σ, β, and γ are constants that define the effective interaction. Their values are determined by fits to experimental nuclear data, results of ab initio models, and astrophysical observations.

At finite temperature, T, the particle number densities

$$\begin{array}{c}\hfill n_i=\frac{g_i}{2{\pi }^2}{\displaystyle \int d}k{k}^2{f}_i(k),\end{array}$$(7)

and the kinetic energy densities

$$\begin{array}{c}\hfill {\tau }_i=\frac{g_i}{2{\pi }^2}{\displaystyle \int d}k{k}^4{f}_i(k),\end{array}$$(8)

are defined in terms of momentum distributions, f_i(k), which describe the thermal population of momentum states according to a Fermi-Dirac distribution,

$$\begin{array}{c}\hfill f_i(k)=\frac{1}{1+exp[({ħ}^2{k}^2/2m_{\mathrm{eff};i}+U_i-{\mu }_i)/T]}·\end{array}$$(9)

Here, g_i = 2 is the spin degeneracy factor, μ_i denotes the chemical potential of the i-particle, and m_eff; i is the effective mass of the i-particle, which is defined through

$$\begin{array}{c}\hfill \frac{1}{m_{\mathrm{eff};i}}=\frac{1}{m_i}+\frac{2}{{ħ}^2}[{\stackrel{\sim }{C}}_{\mathrm{eff}}(n)n\pm {\stackrel{\sim }{D}}_{\mathrm{eff}}(n)n_3],\end{array}$$(10)

where

$$\begin{array}{cc}& {\stackrel{\sim }{C}}_{\mathrm{eff}}(n)=C_{\mathrm{eff}}+[3t_4{n}^{\beta }+t_5(4x_5+5)n^{\gamma }]/16,\hfill \\ \hfill & {\stackrel{\sim }{D}}_{\mathrm{eff}}(n)=D_{\mathrm{eff}}+[-t_4(2x_4+1)n^{\beta }+t_5(2x_5+1)n^{\gamma }]/16.\hfill \end{array}$$(11)

The single-particle potential of the i-particle is U_i = ∂ℋ/∂n_i = U_0i + U_3i + U_effi + U_4i + U_5i.

The locality of Skyrme interactions ensures that nucleon effective masses are independent of temperature. In the zero-temperature limit, m_eff; i corresponds to the Landau effective mass, which is defined in terms of density of single-particle states at the Fermi surface,

$$\begin{array}{c}\hfill \frac{1}{m_{\mathrm{eff};i}}={\frac{1}{{ħ}^2{k}_i}\frac{d{e}_i}{d{k}_i}|}_{k=k_{F;i}},\end{array}$$(12)

where k_F; i is the Fermi momentum of the i-particle. The expressions of the various terms that enter U_i are as follows:

$$\begin{array}{cc}& U_{0_i}=2C_{0}n\pm 2D_0{n}_3,\hfill \end{array}$$(13)

$$\begin{array}{cc}& U_{3_i}=(\sigma +2)C_3{n}^{\sigma +1}+\sigma D_3{n}^{\sigma -1}{n}_3^2\pm 2D_3{n}^{\sigma }{n}_3,\hfill \end{array}$$(14)

$$\begin{array}{cc}& U_{{\mathrm{eff}}_i}=C_{\mathrm{eff}}\tau \pm D_{\mathrm{eff}}{\tau }_3,\hfill \end{array}$$(15)

$$\begin{array}{cc}& U_{4_i}=\frac{t_4}{8}{n}^{\beta -1}[n\tau (2+x_4+\frac{3\beta }{2})-(1+2x_4)(\frac{\beta }{2}{n}_3{\tau }_3+n{\tau }_i)]\hfill \end{array}$$(16)

$$\begin{array}{cc}& U_{5_i}=\frac{t_5}{8}{n}^{\gamma -1}[n\tau (2+x_5+\frac{5\gamma }{2}+2x_5\gamma )+(1+2x_5)\hfill \\ \hfill & \times (\frac{\gamma }{2}{n}_3{\tau }_3+n{\tau }_i)].\hfill \end{array}$$(17)

In Eqs. (10), (13), (14), and (15), the ± sign distinguishes neutrons (+) from protons (−).

The EOS stiffness is determined by pressure,

$$\begin{array}{c}\hfill P=\frac{2}{3}k+h_0+(\sigma +1)h_3+\frac{5}{3}{h}_{\mathrm{eff}}+(\frac{5}{3}+\beta )h_4+(\frac{5}{3}+\gamma )h_5.\end{array}$$(18)

Given that σ, β, and γ are positive (see Table III in Beznogov & Raduta 2024b), it is clear that the pressure depends on the EOS stronger than the energy.

Thermal contributions to state variables are conveniently gauged by taking the difference between the values that the quantity takes at finite and zero temperatures, X_th(n_n, n_p, T) = X(n_n, n_p, T)−X(n_n, n_p, T = 0). In the case of energy density and pressure, we find (Constantinou et al. 2014)

$$\begin{array}{c}\hfill e_{\mathrm{th}}={\displaystyle \sum _{i=n,p}}\frac{{ħ}^2}{2m_{\mathrm{eff};i}}[{\tau }_i(T)-{\tau }_i(T=0)]\end{array}$$(19)

and

$$\begin{array}{c}\hfill P_{\mathrm{th}}={\displaystyle \sum _{i=n,p}}\frac{{ħ}^2}{3m_{\mathrm{eff};i}}(1-\frac{3}{2}\frac{n}{m_{\mathrm{eff};i}}\frac{\partial m_{\mathrm{eff};i}}{\partial n})[{\tau }_i(T)-{\tau }_i(T=0)],\end{array}$$(20)

respectively. Equations (19) and (20) show that e_th and P_th explicitly depend on the effective masses of the nucleons. Equation (20) indicates that P_th additionally depends on the density dependence of the effective masses of the nucleons. However, through τ_i, both the thermal energy density and the thermal pressure depend on every term of the effective interaction.

In the low-temperature limit, this dependence can be derived explicitly using the Sommerfeld expansion. At the lowest order in temperature,

$$\begin{array}{c}\hfill {\tau }_i(T)-{\tau }_i(T=0)=\frac{T^2}{{ħ}^4}{\left(\frac{\pi }{3}\right)}^{2/3}{m}_{\mathrm{eff};i}^2{n}_i^{1/3}.\end{array}$$(21)

Substituting this expansion into Eqs. (19) and (20) gives

$$\begin{array}{cc}& e_{\mathrm{th}}\approx \frac{T^2}{2{ħ}^2}{\left(\frac{\pi }{3}\right)}^{2/3}{\displaystyle \sum _{i=n,p}}{m}_{\mathrm{eff};i}{n}_i^{1/3},\hfill \end{array}$$(22)

$$\begin{array}{cc}& P_{\mathrm{th}}\approx \frac{T^2}{3{ħ}^2}{\left(\frac{\pi }{3}\right)}^{2/3}{\displaystyle \sum _{i=n,p}}[m_{\mathrm{eff};i}-\frac{3}{2}\frac{\partial m_{\mathrm{eff};i}}{\partial n}n]n_i^{1/3}.\hfill \end{array}$$(23)

The expressions above are accurate for T ≲ 10 − 20 MeV, depending on the effective interaction and the baryon number density. In the low temperature limit, e_th depends only on m_eff, while P_th depends on m_eff and its derivative with respect to n.

3. Nuclear matter

In this work, we examined finite-temperature effects for two sets of effective interactions, together with the five effective interactions for which general purpose EOS tables were constructed in Paper I.

The first set (S1) corresponds to the ensemble of effective interactions from run 1 in Beznogov & Raduta (2024b). These interactions were generated within a Markov chain Monte Carlo (MCMC) procedure that constrained both nuclear and NS matter. For NM we constrained: (i) the behavior of saturated SNM via the values of the saturation density (n_sat) and, the energy per particle (E_sat), and compression modulus (K_sat) at n_sat, (ii) the symmetry energy (J_sym) at n_sat, (iii) the energy per neutron (E/A) in PNM at n = 0.08, 0.12 and 0.16 fm⁻³, for which we used χEFT data from Somasundaram et al. (2021; iv) for n ≤ n_l, $0\le m_{\mathrm{eff};n}^{\mathrm{SNM}}/m_n\le 1$ and $0\le m_{\mathrm{eff};n}^{\mathrm{PNM}}/m_n\le 1$, (v) for n ≤ n_l, $v_{F;n}^{\mathrm{SNM}}/c\le 1$ and $v_{F;n}^{\mathrm{PNM}}/c\le 1$ (Duan & Urban 2023), where v_F; i = ℏk_F; i/m_eff; i denotes the Fermi velocity of the i-particle. We required NS EOSs to be: (vi) causal up to the density that corresponds to the central density of the maximum mass configuration (n_c; TOV), (vii) thermodynamically stable, (viii) capable of producing maximum gravitational masses (M_G; TOV) greater than 2 M_⊙. We chose the value of n_l = 0.8 fm⁻³ for the upper density limit for applying conditions (iv) and (v). This value represented a compromise between the extension of the validity domain of the interaction and computational efficiency.

The second set (S2) corresponds to the effective interactions in the set above that: (a) satisfy conditions (iv) and (v) for both neutrons and protons at arbitrary proton fractions 0.0 ≤ Y_p ≤ 0.5 and for n ≤ n_c; TOV, (b) satisfy condition (vi) for arbitrary proton fractions 0.0 ≤ Y_p ≤ 0.5 and for n ≤ n_c; TOV. The number of interactions in S2 represents ≈0.7% of those in S1. Although no explicit constraint was imposed on the thermodynamic stability of NM with arbitrary isospin asymmetry, all interactions in S2 automatically fulfilled this condition.

The interactions BBSk2, BBSk3, BBSk4, and BBSk5 belong to S1 and exhibit extreme behaviors in their effective masses and thermal pressure as functions of density. BBSk2 and BBSk3 lie near the 0.98 (upper) and 0.02 (lower) quantiles of m_eff(n), while BBSk4 and BBSk5 provide large negative and positive thermal pressures. BBSk1 belongs to run 2 of Beznogov & Raduta (2024b), which, beyond the constraints of run 1, also imposes constraints on the neutron effective mass in PNM and on the nucleon effective mass in SNM up to n = 0.16 fm⁻³. For BBSk1, $m_{\mathrm{eff}}^{\mathrm{SNM}}(n)$ corresponds to the “median” in that run. The values of NM parameters can be found in Table 2 of Paper I.

The properties of NM EOSs in the two sets are summarized in Table 1 in terms of medians and 90% confidence intervals (CI; lower and upper quantiles). The additional constraints applied in S2 remove models that exhibit the softest increase in energy per particle (E/A) as a function of n in SNM, as well as those with the steepest rise of E/A(n) in PNM. Consequently, the symmetry energy (E_sym) is substantially reduced for densities in excess of 3n_sat. While most S1 models show a steep increase with density (see Fig. 1 in Beznogov & Raduta 2024b), the median in S2 increases only slightly with density (not shown). This implies that almost half of the S2 models feature a decreasing E_sym(n) at high densities. The exclusion of models with soft SNM EOSs is a consequence of the upper limit imposed on m_eff at densities above n_l (see below). The exclusion of models with stiff PNM EOSs arises from the condition on PNM, which is stiffer than NS matter, to be causal up to densities exceeding n_l. This results in slightly reduced maximum gravitational masses for NSs (see Sect. 4). Table 1 shows that the additional constraints increase the effective neutron mass in the SNM at 0.16 fm⁻³ by 18%.

Figure 1 provides insight into the density dependence of the neutron effective mass in SNM and PNM. In SNM, modulo the neutron-proton mass split, $m_{\mathrm{eff};n}^{\mathrm{SNM}}(n)=m_{\mathrm{eff};p}^{\mathrm{SNM}}(n)$. In Figs. 1–11, the hatched regions denote the 90% CI domains corresponding to S1 (yellow) and S2 (gray). Imposing the additional S2 constraints suppresses interactions with low $m_{\mathrm{eff};n}^{\mathrm{SNM}}$ values at n/n_sat ≲ 4 and high values at larger densities. In other words, while the original set of interactions (S1) favors a rather pronounced U-shaped behavior of $m_{\mathrm{eff};n}^{\mathrm{SNM}}(n)$, the more restricted set (S2) shows a weaker density dependence. Nevertheless, the median and both quantiles in S2 show U-shaped profiles. The additional constraints slightly reduce the dispersion bands of both $m_{\mathrm{eff};n}^{\mathrm{SNM}}(n)$ and $m_{\mathrm{eff};n}^{\mathrm{PNM}}(n)$. For n/n_sat ≲ 4, $m_{\mathrm{eff};n}^{\mathrm{PNM}}$ decreases with n in both sets. At higher densities, $m_{\mathrm{eff};n}^{\mathrm{PNM}}$ in S1 continues to decrease, whereas it becomes density independent in S2. The extra constraints of S2 suppress interactions with low $m_{\mathrm{eff};n}^{\mathrm{PNM}}$ values for n/n_sat ≲ 2 and those with a steep decrease at n/n_sat ≳ 3. In Figs. 1–11, we plot the quantiles and bands up to the density where at least one model in a set violates the S2 conditions.

Fig. 1. Neutron effective mass (meff; n) in SNM (top) and PNM (bottom) in units of bare neutron mass (mn) as a function of density. Medians and upper and lower quantiles at 90% CI of the two sets of models in Sect. 3 are depicted with solid and short dashed curves, respectively. The other curves correspond to the BBSk1 – BBSk5 forces (see Paper I).

Figure 2 shows the density dependence of neutron (m_eff; n) and proton (m_eff; p) effective masses in neutron-rich matter. For both sets, the medians and upper quantiles of m_eff; n and m_eff; p, as well as the lower quantiles of m_eff; p, show U-shaped behavior, similar to $m_{\mathrm{eff};n}^{\mathrm{SNM}}(n)$. The lower quantile of m_eff; n in S1 exhibits two extrema, suggesting that many models yielding low values of this quantity show similar behavior. The additional S2 constraints eliminate all these models. The dispersions of m_eff; n and m_eff; p are only slightly reduced after applying the additional constraints. The m_eff; p curves exhibit a stronger density dependence than the m_eff; n curves.

Fig. 2. Neutron (meff; n) and proton (meff; p) effective masses, expressed in units of bare neutron (mn) and proton (mp) masses, respectively, as functions of density in NM with Yp = 0.2. Line styles and colors are the same as in Fig. 1.

For both sets of interactions, the dispersions of all considered m_eff remain significant even at n/n_sat < 1. Additional information on the diversity of m_eff(n) behaviors accommodated by the MCMC procedure is provided by the BBSk1-BBSk5 interactions. For BBSk2, $m_{\mathrm{eff};n}^{\mathrm{SNM}}$ and both m_eff; n and m_eff; p in neutron-rich matter are U-shaped, whereas $m_{\mathrm{eff};n}^{\mathrm{PNM}}$ shows almost no density dependence. For BBSk3, $m_{\mathrm{eff};n}^{\mathrm{PNM}}$ has an inflection point at n/n_sat ≈ 2 and decreases steeply for n/n_sat > 4. Its $m_{\mathrm{eff};n}^{\mathrm{SNM}}(n)$ and m_eff; p(n) in neutron-rich matter are U-shaped, while m_eff; n in NM with Y_p = 0.2 shows two extrema and only a minor variation over 1 ≲ n/n_sat ≲ 6. For BBSk4, all the effective masses considered show a strong U-shaped density dependence. For BBSk5, all the effective masses considered decrease dramatically with increasing density. Neither BBSk4 nor BBSk5 is compatible with the S2 constraints. The m_eff(n) for BBSk1 is similar in shape to that in BBSk3, but the values differ.

Extra insight into the behavior of the effective masses and thermal response (see Sect. 5) can be obtained following the density dependence of

$$\begin{array}{c}\hfill Q_i=1-\frac{3}{2}\frac{n}{m_{\mathrm{eff};i}}\frac{\partial m_{\mathrm{eff};i}}{\partial n}·\end{array}$$(24)

Figure 3 shows Q_n for the cases considered in Fig. 1. Because most of the interactions in S1 exhibit a U-shaped behavior of $m_{\mathrm{eff};n}^{\mathrm{SNM}}(n)$, for half of the models in this set Q_n^SNM < 0 for n/n_sat ≳ 4.4. The filtering imposed by S2 constraints removes interactions that provide a steep increase of $m_{\mathrm{eff};n}^{\mathrm{SNM}}(n)$ at high density, resulting in a more moderate decrease of Q_n^SNM with density. The median of the S2 models reaches zero only at n ≈ 0.9 fm⁻³. The dominant decrease of $m_{\mathrm{eff};n}^{\mathrm{PNM}}(n)$ with density for the S1 models ensures that even the lower quantile of Q_n^PNM remains positive. The presence of models with exotic behaviors in S1 is indicated by BBSk5 and BBSk4: for BBSk5 (BBSk4), Q_n^SNM and Q_n^PNM increase (decrease) with n. For BBSk4, Q_n^SNM and Q_n^PNM become negative at n ≳ 0.5 fm⁻³ and n ≳ 0.7 fm⁻³, respectively.

Fig. 3. Same as Fig. 1 but for Qn (see Eq. (24)).

4. Neutron star matter

Table 1 lists the values of the key NS parameters corresponding to the two sets of interactions. As discussed in Sect. 3, the requirement for PNM to be causal up to n_c; TOV softens the EOSs, which results in slightly lower values of M_G; TOV and M_B; TOV, as well as higher central density (n_c; TOV), energy density (ρ_c; TOV), and pressure (P_c; TOV) for the most massive configurations. The speed of sound at n_c; TOV also increases. The medians of the radii and tidal deformabilities for 1.4 M_⊙ and 2 M_⊙ NSs are only marginally affected by the additional constraints. This suggests that the softening of PNM EOSs is compensated by stiffening in more neutron-rich cores resulting from the lower values of the symmetry energy.

Figure 4 illustrates the composition of NS matter in terms of the proton fraction, Y_p, as a function of density. The interactions in S1 allow for extremely different compositions. For most models, Y_p increases with n. Models such as BBSk2 and BBSk5, which show the opposite trend, also exist. The strong reduction of the symmetry energy for n/n_sat ≳ 3 in S2 results in interactions with decreasing Y_p(n) being slightly more numerous than those with increasing Y_p(n).

Fig. 4. Proton fraction (Yp) as a function of density in NS matter.

5. Thermal behavior

We now consider the finite-temperature behavior of NM. Following Constantinou et al. (2014, 2015), Raduta et al. (2021), Raduta (2022), we investigated the density dependence of e_th, p_th, thermal (Γ_th) and adiabatic (Γ_S) indexes over wide ranges of temperature and proton fraction. We used the same sets of interactions and specific forces as in the previous section.

Figures 5 and 6 show the density dependence of e_th (Eq. (19)) in PNM and SNM at T = 20 MeV, and in NM with Y_p = 0.2 at T = 5 and 20 MeV, respectively. For all considered (Y_p, T) combinations and n/n_sat ≲ 2, e_th(n) increases. At higher densities, we observe multiple behaviors. For NM with Y_p = 0.2 and for SNM, most interactions predict that e_th(n) increases. BBSk5 is an exception to this trend, showing a steady decrease. In PNM, a rise-and-fall behavior is common for most interactions in S1, as seen in the median and in the predictions for BBSk3 and BBSk5. Models that predict an increase of e_th(n) also exist, as shown by the curves for BBSk2, BBSk4, and the upper quantile. BBSk1, which corresponds to the median behavior of a run that imposes additional constraints on m_eff, also exhibits a rise-and-fall behavior. With the extra conditions in S2, nearly all models show an increasing e_th(n). In all cases, the ordering of e_th at fixed n replicates that of m_eff. The strongest density dependence occurs for BBSk4, which also exhibits the strongest density dependence of m_eff. The relative widths of the uncertainty bands remain constant over the considered temperature range. For S2 at 90% CI, its values at n = 0.8 fm⁻³ are 0.31 (0.30), 0.26 (0.25), and 0.31 (0.30) for Y_p = 0, 0.2 and 0.5 and T = 5 MeV (T = 20 MeV).

Fig. 5. Behavior of eth(n) in PNM and SNM at T = 20 MeV.

Fig. 6. Same as Fig. 6 but for NM with Yp = 0.2 for T = 5 and 20 MeV.

Figure 7 shows the density dependence of p_th (Eq. (20)) in PNM and SNM at T = 20 MeV, while Fig. 8 shows p_th(n) in NM with Y_p = 0.2 at T = 5 and 20 MeV. At high densities, the U-shaped dependence of m_eff(n) (Figs. 1 and 2) leads many interactions to predict a strong decrease of p_th, reaching negative values. Negative p_th occurs more frequently in SNM than in PNM and the behavior of p_th(n) in NM with Y_p = 0.2 resembles that of p_th(n) in SNM more closely than that of p_th(n) in PNM. This indicates that the behavior of p_th is not governed by the most abundant species, but by the species with the strongest density dependence. The most extreme behaviors are associated with BBSk4 and BBSk5 models. At high densities and Y_p ≥ 0.2, even the more restricted set S2 includes a significant number of models with p_th < 0. The relative width of the uncertainty band increases with the proton fraction, but remains constant with temperature. For S2 at 90% CI, the values at n = 0.8 fm⁻³ are 0.72 (0.71), 1.64 (1.62), and 2.56 (2.56) for Y_p = 0, 0.2, and 0.5 and T = 5 MeV (T = 20 MeV).

Fig. 7. Behavior of pth(n) in PNM and SNM at T = 20 MeV.

Fig. 8. Same as Fig. 7 but for NM with Yp = 0.2 for T = 5 and 20 MeV.

The thermal index,

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{th}}=1+\frac{P_{\mathrm{th}}}{e_{\mathrm{th}}},\end{array}$$(25)

is commonly used to assess deviations from ideal gas behavior and to supplement cold EOSs with thermal contributions (Bauswein et al. 2010; Hotokezaka et al. 2013; Endrizzi et al. 2018; Camelio et al. 2019; Weih et al. 2020).

Figure 9 shows the behavior of Γ_th(n) in NM with Y_p = 0.2 at T = 5 and 20 MeV. At low densities, where interactions are negligible, Γ_th → 5/3, which is the ideal gas limit. For n ≲ 0.2 fm⁻³, the dispersion between models is small and the density dependence is weak. Nevertheless, the curves for BBSk4 and BBSk5 clearly show that BBSk4 (BBSk5) provides a Γ_th that decreases (increases) with n. These trends persist at higher densities. For n/n_sat ≳ 2, our models show a significant dependence on the density and interaction. At n = 0.7 fm⁻³ and T = 5 MeV, the most extreme values are −0.2 for BBSk4 and 2.3 for BBSk5. We obtain similar values for T = 20 MeV indicating that the temperature dependence is negligible. Equations (19) and (20) show that, for SNM and PNM, Γ_th does not depend on T. In fact, in these two cases, Γ_th = 1 + 2/3Q_n. The median in S1 (S2) shows a pronounced (moderate) decrease with density, due to the strongly (moderately) U-shaped behavior of m_eff; n and m_eff; p for most models in this set. The negative values of Γ_th arise from negative p_th values. In suprasaturated stellar matter, the dominant pressure contribution comes from nucleons and keeps Γ_th negative even after taking electrons into account. As expected from the behaviors of e_th and p_th, Γ_th for NM with Y_p = 0.2 resembles Γ_th in SNM (not shown), while Γ_th in PNM (not shown) is dissimilar.

Fig. 9. Behavior of Γth(n) in NM with Yp = 0.2 for T = 5 and 20 MeV.

The adiabatic index,

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_S={\frac{\partial lnP}{\partial lnn}|}_S,\end{array}$$(26)

measures the stiffness of the EOS in isentropic processes. Its behavior in PNM and SNM with T = 20 MeV, and in NM with Y_p = 0.2 at T = 5 and 20 MeV, is shown in Figs. 10 and 11, respectively. For temperatures below an interaction- and Y_p-dependent threshold, subsaturated SNM and NM with moderate isospin asymmetry exhibit a liquid-gas phase instability (Ducoin et al. 2006). This instability appears as a back-bending in all Y_i(X_i)_{Yj,  j ≠ i} curves, where X_i and Y_i denote an extensive variable and its conjugated intensive variable, respectively (Ducoin et al. 2006). The divergence of Γ_S in NM with Y_p = 0.2 and T = 5 MeV results from the back-bending behaviors of P(n)_T. The pronounced peak of Γ_S in SNM with T = 20 MeV reflects the proximity to the critical temperature of the liquid-gas phase transition, whose typical values lie in the range 16 MeV ≲ T_C ≲ 20 MeV. However, these instabilities are irrelevant in stellar matter, which is stabilized by clusterization (see Paper I). For n/n_sat ≥ 2, most models show a smooth and weak variation of Γ_S(n). For BBSk2 and BBSk4, Γ_S shows a sudden drop at densities near the validity limit of the models. With the extra conditions in S2, most of the models predict a slow decrease in Γ_S(n), and the dispersion among the models is relatively low.

Fig. 10. Behavior of ΓS(n) in PNM and SNM at T = 20 MeV.

6. Stability of proto-neutron stars and remnants of binary neutron star mergers

The maximum gravitational mass of isentropic compact stars is relevant for BH formation in failed CCSNe. Numerical simulations by Schneider et al. (2020) show that the onset of collapse occurs when the gravitational mass of a PNS exceeds the maximum gravitational mass corresponding to the entropy value that characterizes most of the PNS volume.

In the absence of accretion or mass loss, the baryonic mass is conserved, making it a useful quantity for analyzing the stability of PNSs and BNS merger remnants. If the baryonic mass of a PNS (or remnant) exceeds the maximum baryonic mass of the cold beta-equilibrated NS, the PNS (or remnant) is unstable to collapse into a BH, regardless of any temporary stabilizing mechanisms.

The relation between the maximum gravitational mass of isentropic stars, $M_{\mathrm{G};\max }^{(S/A)}$, and the maximum gravitational mass of cold-catalyzed configurations, M_G; TOV, was examined by Raduta et al. (2020) and Wei et al. (2021), using a large number of phenomenological and microscopic models. Their results show that for models based on the covariant density functional theory of NM and for microscopic variational models, $M_{\mathrm{G};\max }^{(S/A)}>M_{\mathrm{G};\mathrm{TOV}}$ (for all values of S/A), while the relation is inverted for models based on the microscopic Brueckner-Hartree-Fock theory. Skyrme-like interactions with m_eff equal to the bare mass and K_sat within the currently accepted range also yield $M_{\mathrm{G};\max }^{(S/A)}>M_{\mathrm{G};\mathrm{TOV}}$ for all values of S/A. The behavior of the maximum baryonic mass of hot stars is ambiguous. According to Raduta et al. (2020), models based on the covariant theory of NM can predict either an increase or a decrease of $M_{\mathrm{B};\max }^{(S/A)}$ with the entropy per baryon.

Figure 12 illustrates $M_{\mathrm{G};\max }^{(S/A)}$ and $M_{\mathrm{B};\max }^{(S/A)}$ as functions of the entropy per baryon (S/A), as predicted by the BBSk1, BBSk2, and BBSk5 effective interactions. For simplicity, we also assumed that hot compact objects have constant radial profiles of Y_p, and we arbitrarily varied these values between 0.06 ≤ Y_p ≤ 0.3, producing the bands shown in the figure. For BBSk1, both $M_{\mathrm{G};\max }^{(S/A)}$ and $M_{\mathrm{B};\max }^{(S/A)}$ depend strongly on Y_p, whereas for BBSk2 and BBSk5 this dependence is rather weak. This difference arises from the high symmetry energy in BBSk1 and the low values in BBSk2 and BBSk5 (see Sect. 3). The three interactions exhibit different behaviors for $M_{\mathrm{G};\max }^{(S/A)}$ and $M_{\mathrm{B};\max }^{(S/A)}$ as functions of S/A. BBSk1 and BBSk2 predict that up to a certain value, $M_{\mathrm{G};\max }^{(S/A)}$ remains constant, while $M_{\mathrm{B};\max }^{(S/A)}$ decreases before both increase; for BBSk5, both quantities increase with S/A but at different slopes. For BBSk1, $M_{\mathrm{G};\max }^{(S/A)}\gtrless M_{\mathrm{G};\mathrm{TOV}}$ ($M_{\mathrm{B};\max }^{(S/A)}\gtrless M_{\mathrm{B};\mathrm{TOV}}$) depending on Y_p (Y_p and S/A). For BBSk2, for all Y_p and S/A, $M_{\mathrm{G};\max }^{(S/A)}\gtrsim M_{\mathrm{G};\mathrm{TOV}}$ and $M_{\mathrm{B};\max }^{(S/A)}\lesssim M_{\mathrm{B};\mathrm{TOV}}$, the latter being highly unusual. For BBSk5, for all Y_p and S/A, $M_{\mathrm{G};\max }^{(S/A)}>M_{\mathrm{G};\mathrm{TOV}}$ and $M_{\mathrm{B};\max }^{(S/A)}>M_{\mathrm{B};\mathrm{TOV}}$.

Fig. 11. The same as Fig. 11 but for NM with Yp = 0.2 for T = 5 and 20 MeV.

Fig. 12. Maximum gravitational (bottom) and baryonic (top) masses versus entropy per baryon (S/A) for isentropic PNSs with constant profiles of 0.06 ≤ Yp ≤ 0.3, as predicted by BBSk1 (black), BBSk2 (cyan) and BBSk5 (magenta) forces. Horizontal lines indicate the corresponding values, MG; TOV and MB; TOV, for cold catalyzed NSs.

Fig. 13. Central temperature of the maximum mass configuration ($T_{\mathrm{c};\max }^{(S/A)}$) versus entropy per baryon (S/A) for the same cases as in Fig. 12.

Using the criterion introduced at the beginning of this section, we conclude that: (i) in a failed CCSN, a PNS with S/A ≳ 2.6 built on BBSk1 will collapse earlier than a PNS built on BBSk5. (ii) a PNS with S/A ≲ 2.2 built on BBSk5 will collapse earlier than a PNS built on BBSk2.

We next consider the criterion based on the maximum baryonic masses. Here, the difference between the EOSs is even more pronounced. For BBSk5, the situation is standard: thermal support stabilizes the star against collapse. The hotter the star, the higher the maximum supported baryonic mass. For BBSk2, the situation is opposite: a cold, beta-equilibrated star supports higher maximum baryonic masses than hot stars. This means that, without other stabilizing mechanisms, the PNS (or remnant) exceeding $M_{\mathrm{B};\max }^{(S/A)}$ will promptly collapse into a BH, even if its mass is below the maximum baryonic mass of a cold NS. This outcome is a direct consequence of p_th being negative (in some T − n − Y_p domain) for BBSk2, while p_th remains positive for BBSk5 (see Sect. 5). A similar effect is also observed in some Brueckner-Hartree-Fock models (Lu et al. 2019).

In agreement with Schneider et al. (2019, 2020), Yasin et al. (2020), Andersen et al. (2021), Fields et al. (2023), our results confirm that the evolution of astrophysical phenomena that involve hot dense matter probes the nucleon effective mass and its density dependence. However, the effective interactions used here have more sophisticated m_eff(n) dependencies, making the link with m_eff(n_sat) non-trivial. Additionally, for models with high symmetry energy at high densities (e.g., BBSk1), the stability of hot compact objects also depends on the Y_p profile.

Figure 13 investigates the central temperature of the maximum mass configuration ($T_{\mathrm{c};\max }^{(S/A)}$) as a function of S/A for the cases considered in Fig. 12. For S/A ≲ 3.5, all forces predict that $T_{\mathrm{c};\max }^{(S/A)}$ increases with S/A; the spread of data also increases with S/A. For all effective interactions, the highest (lowest) value of $T_{\mathrm{c};\max }^{(S/A)}$ corresponds to the lowest (highest) value of Y_p. The maximum spread for given S/A corresponds to BBSk1 and is due to its high symmetry energy at densities above 3n_sat. Additionally, high values of m_eff at n_c; max are correlated with low values of $T_{\mathrm{c};\max }^{(S/A)}$.

7. Conclusions

A large number of ab initio constrained EOS models generated within a Bayesian analysis (Beznogov & Raduta 2024b) was used to investigate the thermal response of dense NM over domains of density, temperature, and proton fraction relevant for the evolution of astrophysical phenomena involving hot compact objects. These models rely on the Brussels parameterization of the Skyrme effective interaction, which offers two major advantages. First, it is flexible enough to allow for widely different behaviors in the suprasaturation regime, including a density dependence of the nucleonic effective mass in qualitative agreement with predictions from microscopic calculations with three-body forces (Baldo et al. 2014; Shang et al. 2020; Somasundaram et al. 2021). Second, the availability of analytical expressions for most thermodynamic and microscopic quantities makes it possible to assess the role that various NM parameters play at finite temperature.

Insufficient knowledge of the dense matter EOS, widely discussed in the literature, translates into thermal responses that differ significantly from one model to another. In general, effective interactions that yield a U-shaped behavior of m_eff(n) produce EOSs with substantially lower thermal pressure than those with a monotonic decrease of m_eff(n). For densities above several times n_sat, the former models feature negative thermal pressures. Negative thermal pressures (p_th < 0) were previously obtained within χEFT calculations (Keller et al. 2023).

The stiffness of the EOS at finite temperatures impacts the stability of hot stars, with consequences for the fates of PNSs and BNS mergers. Three particular interactions in our set, which manifest different behaviors of m_eff(n), lead to different dependencies of the maximum gravitational and baryonic masses as a function of S/A. In particular, a model with m_eff decreasing with n predicts that both $M_{\mathrm{G};\max }^{(S/A)}$ and $M_{\mathrm{B};\max }^{(S/A)}$ exceed their counterparts in cold, beta-equilibrated NSs. Two other models, with U-shaped m_eff(n) behaviors, show that, depending on S/A, Y_p, and the effective interaction, $M_{\mathrm{G};\max }^{(S/A)}\gtrless M_{\mathrm{G};\mathrm{TOV}}$ and $M_{\mathrm{B};\max }^{(S/A)}\gtrless M_{\mathrm{B};\mathrm{TOV}}$, with $M_{\mathrm{B};\max }^{(S/A)}<M_{\mathrm{B};\mathrm{TOV}}$ being the most notable result. The stability with respect to the collapse of models with high E_sym at high densities also depends on the Y_p-profile.

The use of EOS tables based on Brussels-Skyrme interactions (such as those introduced in Paper I and available on COMPOSE) in numerical simulations will make the correlation between m_eff and the evolution of astrophysical phenomena more difficult to establish than in studies such as (Schneider et al. 2019, 2020; Yasin et al. 2020; Andersen et al. 2021; Fields et al. 2023). Nevertheless, this is a necessary step toward a better understanding of the properties of dense matter, including the possible presence of non-nucleonic degrees of freedom.

The present work is a follow-up of the work of Raduta et al. (2024), which considers the thermal response of a set of EOSs that rely on the covariant density functional theory of nuclear matter and is generated within a Bayesian inference of the EOS of dense matter (Beznogov & Raduta 2023). Together with that study, the present analysis contributes to a deeper understanding of the role that the EOS plays in the evolution of hot compact objects.

Acknowledgments

We acknowledge support from a grant from the Ministry of Education and Research, CNCS/CCCDI–UEFISCDI, Project No. PN-IV-P1-PCE-2023-0324; partial support from Project No. PN 23 21 01 02 is also acknowledged. The two authors have contributed equally to this work. Note added in proofs. In Paper I, the thermodynamic identities in the paragraph following Eq. (38) should read: F + PV − ∑_{i = n, p, {A, Z},el}μ_iN_i = 0 and F_i + P_iV − μ_iN_i ≠ 0, respectively.

References

All Tables

All Figures

	Fig. 1. Neutron effective mass (meff; n) in SNM (top) and PNM (bottom) in units of bare neutron mass (mn) as a function of density. Medians and upper and lower quantiles at 90% CI of the two sets of models in Sect. 3 are depicted with solid and short dashed curves, respectively. The other curves correspond to the BBSk1 – BBSk5 forces (see Paper I).
In the text

	Fig. 2. Neutron (meff; n) and proton (meff; p) effective masses, expressed in units of bare neutron (mn) and proton (mp) masses, respectively, as functions of density in NM with Yp = 0.2. Line styles and colors are the same as in Fig. 1.
In the text

	Fig. 3. Same as Fig. 1 but for Qn (see Eq. (24)).
In the text

	Fig. 4. Proton fraction (Yp) as a function of density in NS matter.
In the text

	Fig. 5. Behavior of eth(n) in PNM and SNM at T = 20 MeV.
In the text

	Fig. 6. Same as Fig. 6 but for NM with Yp = 0.2 for T = 5 and 20 MeV.
In the text

	Fig. 7. Behavior of pth(n) in PNM and SNM at T = 20 MeV.
In the text

	Fig. 8. Same as Fig. 7 but for NM with Yp = 0.2 for T = 5 and 20 MeV.
In the text

	Fig. 9. Behavior of Γth(n) in NM with Yp = 0.2 for T = 5 and 20 MeV.
In the text

	Fig. 10. Behavior of ΓS(n) in PNM and SNM at T = 20 MeV.
In the text

	Fig. 11. The same as Fig. 11 but for NM with Yp = 0.2 for T = 5 and 20 MeV.
In the text

	Fig. 12. Maximum gravitational (bottom) and baryonic (top) masses versus entropy per baryon (S/A) for isentropic PNSs with constant profiles of 0.06 ≤ Yp ≤ 0.3, as predicted by BBSk1 (black), BBSk2 (cyan) and BBSk5 (magenta) forces. Horizontal lines indicate the corresponding values, MG; TOV and MB; TOV, for cold catalyzed NSs.
In the text

	Fig. 13. Central temperature of the maximum mass configuration ($T_{\mathrm{c};\max }^{(S/A)}$) versus entropy per baryon (S/A) for the same cases as in Fig. 12.
In the text