11institutetext: Dipartimento di Fisica “E. Fermi”, Università di Pisa, Largo
B. Pontecorvo, 3 I-56127 Pisa, Italy 22institutetext: INFN, Sezione di Pisa,
Largo B. Pontecorvo, 3 I-56127 Pisa, Italy 33institutetext: Dipartimento di
Fisica, Universitá di Trento, Via Sommarive 14, 38123 Trento, Italy
44institutetext: INFN-TIFPA, Trento Institute for Fundamental Physics and
Applications, via Sommarive 14, I-38123 Trento, Italy

# Microscopic equation of state of hot nuclear matter for numerical relativity
simulations

Domenico Logoteta, email<EMAIL_ADDRESS>2 2 Albino Perego,
email<EMAIL_ADDRESS>Ignazio Bombaci, email:
<EMAIL_ADDRESS>2 2

(Received ………; accepted ……)

###### Abstract

Context. A precise understanding of the equation of state (EOS) of dense and
hot matter is key to modeling relativistic astrophysical environments,
including core-collapse supernovae (CCSNe), protoneutron star (PNSs)
evolution, and compact binary mergers.

Aims. In this paper, we extend the microscopic zero-temperature BL (Bombaci
and Logoteta) nuclear EOS to finite temperature and arbitrary nuclear
composition. We employ this new EOS to describe hot $\beta$-stable nuclear
matter and to compute various structural properties of nonrotating PNS. We
also apply the EOS to perform dynamical simulations of a spherically symmetric
CCSN.

Methods. The EOS is derived using the finite temperature extension of the
Brueckner–Bethe–Goldstone quantum many-body theory in the
Brueckner–Hartree–Fock approximation. Neutron star properties are computed by
solving the Tolman–Oppenheimer–Volkoff structure equations numerically. The
sperically symmetric CCSN simulations are performed using the AGILE-IDSA code.

Results. Our EOS models are able to reproduce typical features of both PNS and
spherically symmetric CCSN simulations. In addition, our EOS model is
consistent with present measured neutron star masses and particularly with the
masses: $M=2.01\pm 0.04\,M_{\odot}$ and $M=2.14^{+0.20}_{-0.18}\,M_{\odot}$ of
the neutron stars in PSR J0348+0432 and PSR J0740+6620 respectively. Finally,
we suggest a feasible mechanism to produce low-mass black holes ($M\sim
2M_{\odot}$) that could have far-reaching consequences for interpreting the
gravitational wave event GW190814 as a BH–BH merger.

###### Key Words.:

Dense matter – Equation of state – Stars: neutron

##  1 Introduction

One of the main inputs for the study of high-energy astrophysical systems,
like core-collapse supernovae (CCSNe; see, e.g., Burrows, 2013; Janka et al.,
2016, for recent reviews), protoneutron stars (PNSs; e.g., Prakash et al.,
1997; Roberts & Reddy, 2017), and binary compact stars mergers including at
least one neutron star (e.g., Baiotti & Rezzolla, 2017; Shibata & Hotokezaka,
2019; Radice et al., 2020), is the relationship between matter pressure ($P$),
energy density ($\varepsilon$), baryonic density ($n$), and temperature ($T$).
The relationships between these physical quantities, together with information
on the matter composition, make up the equation of state (EOS) of the system
(Prakash et al., 1997; Lattimer & Prakash, 2016; Oertel et al., 2017). In
recent years, increasing interest surrounding the above-mentioned
astrophysical phenomena has stimulated several studies that aimed to provide
as accurate a description as possible of matter under extreme conditions of
density and temperature. Because of the huge density and temperature
variations expected during the dynamical evolution of collapsing stellar cores
and merging neutron stars (Fischer et al., 2014; Hanauske & Bovard, 2018;
Perego et al., 2019), the derivation of an EOS suitable for their description
requires an enormous theoretical and numerical effort. Such an EOS needs to be
capable of describing matter in a wide range of conditions, ranging from
homogeneous nuclear matter down to a mixture of fully ionized atoms. To this
end, several finite-temperature, composition-dependent EOSs have been
developed according to different frameworks and philosophies; among the many,
we only mention EOS models based on nonrelativistic Skyrme interactions
(Lattimer & Swesty, 1991; Schneider et al., 2017), variational approaches
(Togashi et al., 2017), and relativistic mean field (RMF) models (Sugahara &
Toki, 1994; Steiner et al., 2013c; Typel et al., 2010; Fattoyev et al., 2010).

The recent detection of gravitational waves (GWs) produced in binary neutron
star mergers (BNSMs; see, e.g., Abbott et al., 2017, 2020a) has provided a
strong boost to the research on dense matter physics. Gravitational wave
astronomy has thus opened up the possibility to explore matter under extreme
conditions in regimes that would be inaccessible by experiments on earth.
Gravitational waves from the inspiral phase of BNSMs (Abbott et al., 2018),
especially when combined with robust lower limits on the maximum mass of
neutron stars (NSs), can set stringent constraints on the EOS of cold (i.e.,
$T=0$) nuclear matter up to two to three times nuclear saturation density
$n_{0}=0.16\leavevmode\nobreak\ \mathrm{fm^{-3}}$ (corresponding to $\sim
2.6\times 10^{14}\mathrm{g\leavevmode\nobreak\ cm^{-3}}$). The post-merger GW
signal and the associated electromagnetic signature could provide unique
information on the EOS of hot ($T\geq 10\,\mathrm{MeV}$) and dense ($n\geq
3n_{0}$) stellar matter (see e.g., Radice et al., 2018a, b). In particular,
post-merger GWs could carry information about the possible appearance of new
particle species (e.g., hyperons) or new phases of matter (e.g., deconfined
quark matter; see e.g., Shibata et al., 2005; Bauswein & Janka, 2012;
Hotokezaka et al., 2013; Bernuzzi et al., 2015; Rezzolla & Takami, 2016;
Bernuzzi et al., 2016; Zappa et al., 2018; Bauswein et al., 2019; Most et al.,
2020). In addition, the combined mass-radius measurements provided by Neutron
Star Interior Composition Explorer (NICER; see, e.g., Miller et al., 2019;
Riley et al., 2019) set other independent, robust constraints on the high-
density EOS.

The purpose of the present work is to extend the microscopic zero-temperature
EOS of nucleonic matter derived by Bombaci & Logoteta (2018) to finite
temperature. Here ”microscopic” refers to the fact that the EOS has been
derived within a specified quantum many-body approach starting from two-body
and three-body nuclear interactions which reproduce nucleon-nucleon (NN)
scattering data, the experimental binding energy of light ($A=3,4$) atomic
nuclei, and the empirical saturation properties of nuclear matter.

Various zero-temperature microscopic EOS models are available in the
literature (e.g., Wiringa et al., 1988; Baldo et al., 1997; Akmal et al.,
1998) and have mainly been employed to describe neutron star structure. Some
of these EOSs have also been used in BNSM simulations adding a thermal
contribution from the so-called $\Gamma$-law prescription (Bauswein et al.,
2010; Kiuchi et al., 2014).

The nuclear many-body problem at finite temperature has been considered by
several authors and various techniques have been developed to handle this very
complicated topic. Among the most popular ones are those based on Green’s
function method (Frick & Müther, 2003; Frick et al., 2005; Rios et al., 2006),
on Bloch-De Dominicis (BD) diagrammatic expansion (Bloch & de Dominicis, 1958,
1959a, 1959b), and, more recently, on many-body perturbation theory (Fiorilla
et al., 2012). Thus several microscopic EOSs have been extended to finite
temperature (e.g., Bombaci et al., 1994; Baldo & Ferreira, 1999; Zuo et al.,
2004; Burgio et al., 2011; Wellenhofer et al., 2015; Drischler et al., 2016;
Carbone et al., 2018; Carbone & Schwenk, 2019).

In the present work, as in Bombaci & Logoteta (2018), we model the neutron
star core as a uniform charge-neutral fluid of neutrons, protons, electrons,
and muons. The EOS is calculated for arbitrary lepton fractions, densities,
and temperatures as is indeed required by CCSNe and BNSM simulations in order
to describe the behavior of matter which, in the more general case, may be out
of equilibrium with respect to the weak interactions. We note that even
considering this “simplified” picture which does not take into account
”exotic” degrees of freedom (like hyperons or deconfined quarks), the
determination of the EOS from the underlying nuclear interactions is still a
very challenging theoretical problem. In fact, it is necessary to calculate
the EOS under extreme conditions of density, neutron-proton asymmetry, and
temperature in a regime where the EOS is poorly constrained by nuclear data
and experiments.

In order to describe the interactions between nucleons, consistently with
Bombaci & Logoteta (2018), we use the version $b$ of the local chiral two-
nucleon potential reported by Piarulli et al. (2016, see Table 2 of that
work), supplemented with the three-nucleon force discussed by Logoteta et al.
(2016b). Finally, using the Brueckner–Bethe–Goldstone (BBG; see, e.g., Day,
1967) many-body theory within the Brueckner–Hartree–Fock (BHF) approximation
extended to finite temperature, we calculate the Helmholtz free energy of
asymmetric nuclear matter. All the other relevant thermodynamic quantities
(pressure, energy density, etc.) can then be obtained using standard
thermodynamic relations. This finite temperature EOS will be provided in
parametric form, ready for use in numerical simulations of BNSMs and CCSNe. In
the particular case of $\beta$-stable nuclear matter at finite temperature, we
use our new EOS to compute various structural properties of nonrotating
protoneutron stars considering several different stellar conditions such as
isoentropic and isothermal neutron star sequences. In addition, we apply this
new EOS to perform some spherically symmetric (1D) CCSN simulations.

We stress that in the present work we are extending the zero-temperature EOS
of Bombaci & Logoteta (2018) to finite temperature using modern two-body and
three-body nuclear interactions based on chiral effective field theory, while
employing the existing finite-temperature BHF many-body approach.

The paper is organized as follows: in section 2 we describe the extension of
the BBG many-body theory at finite temperature. In section 3, we present the
results of our new finite-temperature, microscopic EOS for stellar matter.
Section 4 is devoted to describing several astrophysical applications of our
new EOS, namely: the construction of the EOS in $\beta$-equilibrium both in
the isoentropic and the isothermal cases, the calculation protoneutron stars
structure using these EOSs, and finally its application to spherically
symmetric simulations of CCSNe. In the last section, we summarize our results
and outline the main conclusions of the present work.

##  2 Brueckner–Bethe–Goldstone many-body theory at finite temperature

The Bloch-De Dominicis diagrammatic expansion of the grand canonical
thermodynamic potential (Bloch & de Dominicis, 1958, 1959a, 1959b) is
significant in that it corresponds, in the zero temperature limit, to the BBG
expansion of the ground-state energy (Baldo & Ferreira, 1999). Thus, BD
expansion can be viewed as the natural extension of the BBG many-body theory
at finite temperature. In particular, Baldo & Ferreira (1999) showed that the
dominant terms in the BD expansion are the ones that correspond to the (zero
temperature) BBG diagrams, where the temperature is introduced in the
occupation numbers only, represented by Fermi distribution functions (Bombaci
et al., 1993, 1994).

More specifically, the BBG many-body theory is based on a linked cluster
expansion (the hole-line expansion) of the energy per nucleon
$\widetilde{E}\equiv E/A$ of nuclear matter. The various terms of the
expansion can be represented by Goldstone diagrams grouped according to the
number of independent hole lines (i.e., lines representing empty single
particle states in the Fermi sea). The basic ingredient in this approach is
the Brueckner reaction matrix $G$, which sums, in a closed form, the infinite
series of the so-called ladder diagrams and takes into consideration the
short-range strongly repulsive part of the nucleon–nucleon (NN) interaction.

In the general case of asymmetric nuclear matter with neutron number density
$n_{n}$, proton number density $n_{p}$, total nucleon number density
$n=n_{n}+n_{p}$, and isospin asymmetry (asymmetry parameter),

$\beta=\frac{n_{n}-n_{p}}{n}\,,$ (1)

the reaction matrix depends on the isospin third components $\tau$ and
$\tau^{\prime}$ of the two interacting nucleons ($\tau,\tau^{\prime}=n,p$).
Thus, there are different $G$-matrices describing the $nn$, $pp$, and $np$ in
medium effective interactions. These are obtained by solving the generalized
Bethe–Goldstone equations

$G_{\tau\tau^{\prime}}(\omega,T)=v+v\sum_{k_{a},k_{b}}\frac{\mid\@vec{k_{a}},\@vec{k_{b}}\rangle
Q_{\tau\tau^{\prime}}\langle\@vec{k_{a}},\@vec{k_{b}}\mid}{\omega-\epsilon_{\tau}(k_{a},T)-\epsilon_{\tau^{\prime}}(k_{b},T)+i\eta}\,G_{\tau\tau^{\prime}}(\omega,T)\;,$
(2)

where ${v}$ is the bare NN interaction (or a density-dependent two-body
effective interaction where averaged three-nucleon forces are introduced) and
the quantity $\omega$ is the so-called starting energy. In the present work we
consider spin unpolarized nuclear matter, and therefore in equation (2) and in
the following equations, we drop the spin indices to simplify the formalism
111Spin polarized nuclear matter within the BHF approach has been considered
by, e.g., Vidaña & Bombaci (2002) and Bombaci et al. (2006).. The operator

$\mid\@vec{k_{a}},\@vec{k_{b}}\rangle\,Q_{\tau\tau^{\prime}}\,\langle\@vec{k_{a}},\@vec{k_{b}}\mid\leavevmode\nobreak\
\equiv\leavevmode\nobreak\
\mid\@vec{k_{a}},\@vec{k_{b}}\rangle\,Q_{\tau\tau^{\prime}}(\@vec{k_{a}},\@vec{k_{b}};n,\beta,T)\,\langle\@vec{k_{a}},\@vec{k_{b}}\mid\,=\,\Big{[}1-f_{\tau}(\@vec{k_{a}},\tilde{\mu}_{\tau},T)\Big{]}\;\Big{[}1-f_{\tau{{}^{\prime}}}(\@vec{k_{b}},\tilde{\mu}_{\tau},T)\Big{]}$
(3)

is the finite temperature Pauli operator, which can be written in terms of the
Fermi distribution function

$f_{\tau}(k,\tilde{\mu}_{\tau},T)=\frac{1}{1+{\rm
exp}[(\epsilon_{\tau}(k,T)-\tilde{\mu}_{\tau})/T]}$ (4)

of the two nucleons in the intermediate scattering states. In the above
expression, $\tilde{\mu}_{\tau}$ is the so-called auxiliary chemical potential
and is obtained using the normalization condition

$n_{\tau}=\sum_{k}f_{\tau}(k,\tilde{\mu}_{\tau},T)\,,$ (5)

where $n_{\tau}$ is the number density of neutrons ($\tau=n$) or protons
($\tau=p$). The Bethe–Goldstone equation therefore describes the two nucleons
scattering in the presence of other nucleons, and the Brueckner $G$-matrix
represents the effective interaction between two nucleons in the nuclear
medium and properly takes into account the short-range correlations arising
from the strongly repulsive core in the bare NN interaction.

The single-particle energy $\epsilon_{\tau}(k,T)$ of a nucleon ($\tau=n,p$)
with momentum $\@vec{k}$ and mass $m_{\tau}$ is given by

$\epsilon_{\tau}(k,T)=\frac{\hbar^{2}k^{2}}{2m_{\tau}}+U_{\tau}(k,T)\ ,$ (6)

where $U_{\tau}(k,T)$ is a single-particle potential that represents the mean
field felt by a nucleon due to its interaction with the surrounding nucleons
in the medium. In the BHF approximation of the BBG theory, $U_{\tau}(k,T)$ is
calculated through the real part of the on-energy-shell $G$-matrix (Bethe et
al., 1963; Hüfner & Mahaux, 1972) and is given by:

$U_{\tau}(k,T)=\sum_{\tau^{\prime}}\sum_{k^{\prime}}\mbox{Re}\
\langle\@vec{k},\@vec{k^{\prime}}\mid
G_{\tau\tau^{\prime}}(\omega^{*},T)\mid\@vec{k},\@vec{k^{\prime}}\rangle_{a}\
f_{\tau^{\prime}}(\@vec{k^{\prime}},\tilde{\mu}_{\tau},T)\;,$ (7)

where
$\omega^{*}=\epsilon_{\tau}(k,T)+\epsilon_{\tau^{\prime}}(k^{\prime},T),$ and
the matrix elements are properly antisymmetrized. We make use of the so-called
continuous choice (Jeukenne et al., 1976; Grangè et al., 1987; Baldo et al.,
1990, 1991) for the single-particle potential $U_{\tau}(k,T)$ when solving the
Bethe–Goldstone equation. As shown by Song et al. (1998) and Baldo et al.
(2000), the contribution of the three-hole-line diagrams to the energy per
nucleon $E/A$ is minimized in this prescription and a faster convergence of
the hole-line expansion for $E/A$ is achieved with respect to the gap choice
for $U_{\tau}(k,T)$.

We finally performed the usual angular average of the Pauli operator and of
the energy denominator (Grangè et al., 1987; Baldo et al., 1991) in the
Bethe–Goldstone equation (2). We note that, as opposed to the zero-temperature
case, the average of the Pauli operator should be performed numerically at
finite temperature. In this way, the Bethe-Goldstone equation can be expanded
in partial waves. In the present work, we consider partial wave contributions
up to a total two-body angular momentum $J_{\rm max}=8$. We verified that the
inclusion of partial waves with $J_{\rm max}>8$ does not appreciably change
our results.

In this scheme, Eqs. (2)–(7) have to be solved self-consistently using an
iterative numerical procedure. Once a self-consistent solution is achieved,
the energy per nucleon of asymmetric nuclear matter is given by:

$\displaystyle\widetilde{E}(n,\beta,T)=\frac{E}{A}(n,\beta,T)$
$\displaystyle=$
$\displaystyle\frac{1}{A}\sum_{\tau}\sum_{k}\left[\frac{\hbar^{2}k^{2}}{2m_{\tau}}+\frac{1}{2}U_{\tau}(k,T)\right]f_{\tau}(k,\tilde{\mu}_{\tau},T)\;.$
(8)

The total entropy $S$ has been calculated in the approximation of a mixture of
two ideal Fermi gases of quasi-particles in the mean fields
$U_{\tau}(k;n,\beta,T)$ (Bombaci et al., 1993):

$\displaystyle S(n,\beta,T)$ $\displaystyle=$ $\displaystyle-
k_{B}\sum_{\tau}\sum_{k}\bigg{[}f_{\tau}(k,\tilde{\mu}_{\tau},T)\ {\rm
ln}\,\Big{(}f_{\tau}(k,\tilde{\mu}_{\tau},T)\Big{)}+\Big{(}1-f_{\tau}(k,\tilde{\mu}_{\tau},T)\Big{)}\
{\rm ln}\,\Big{(}1-f_{\tau}(k,\tilde{\mu}_{\tau},T)\Big{)}\bigg{]}\,.$ (9)

Subsequently, the Helmholtz fee energy $F$ can be calculated using the
standard thermodynamic relation:

$\displaystyle F=E-TS\,.$ (10)

In this scheme, the bare two-body nuclear interaction ${v}$ is the only
physical input we need for the numerical solution of the Bethe-Goldstone
equation. As already mentioned, in the present work, as two-body nuclear
interaction, we employ the model $b$ of the local chiral N3LO$\Delta$
interaction derived by Piarulli et al. (2016).

It is however well established that within the most advanced nonrelativistic
quantum many-body approaches it is not possible to reproduce the empirical
saturation point of symmetric nuclear matter (SNM) when using two-body nuclear
interactions only. The saturation points obtained using different NN
potentials lie within a narrow band called the Coester band (Coester et al.,
1970; Day, 1981), with either an overly high saturation density or an overly
low binding energy ($B=-E/A$) compared to the empirical values. In particular,
SNM turns out to be overbound with an overly high saturation density when
using modern high-precision NN potentials, fitting NN scattering data up to an
energy of $350$ MeV, with a $\chi^{2}$ per datum close to $1$ (Li et al.,
2006). For the nuclear matter case, as in the case of few-nucleon systems
(Kalantar-Nayestanak et al., 2012; Hammer et al., 2013; Binder et al., 2016),
three-nucleon-forces (TNFs) are considered to be the missing physical effect
of the whole picture. The inclusion of TNFs is therefore required in order to
reproduce a realistic saturation point (Friedman & Pandharipande, 1981; Baldo
et al., 1997; Akmal et al., 1998; Li et al., 2008; Zuo et al., 2014; Kievsky
et al., 2018). In addition, TNFs are crucial in the case of dense
$\beta$-stable nuclear matter to obtain a stiff EOS (Baldo et al., 1997; Akmal
et al., 1998; Li & Schulze, 2008; Chamel et al., 2011) compatible with the
measured masses 222Stellar masses are given in units of solar mass
$M_{\odot}=1.989\times 10^{33}$ g $M=1.97\pm 0.04\,M_{\odot}$ (Demorest et
al., 2010), $M=2.01\pm 0.04\,M_{\odot}$ (Antoniadis et al., 2013), and
$M=2.14^{+0.20}_{-0.18}\,M_{\odot}$ (Cromartie et al., 2019) of the neutron
stars in PSR J1614-2230, PSR J0348+0432, and PSR J0740+6620 respectively.

In the present work we use the same TNFs adopted by Bombaci & Logoteta (2018)
namely the N2LO$\Delta$. The parameters of such force were determined in
Logoteta et al. (2016b) in order to reproduce good saturation properties of
SNM. The specific parametrization that we use in the present work is the
N2LO$\Delta 1$ (Logoteta et al., 2016b; Bombaci & Logoteta, 2018). As in
Logoteta et al. (2016b) and in Bombaci & Logoteta (2018), we derive an
effective density-dependent two-body force $v_{NN}^{\rm eff}(n)$ starting from
the original three-body force by averaging over one of the three nucleons
(Holt et al., 2010). The Bethe–Goldstone equation (2) is then solved adding
this effective density-dependent two-body force to the bare NN interaction.
This strategy allows us to avoid the need to solve three-body Faddeev
equations in the nuclear medium (Bethe–Faddeev equations) (Bethe, 1965;
Rajaraman & Bethe, 1967) which are currently too complicated to deal with.

We note that adopting the two- and three-nucleon interactions described above,
the resulting EOS is consistent with the one derived in Bombaci & Logoteta
(2018) at zero temperature.

Finally, we want to emphasize that the compatibility between the BBG many-body
method and other two independent many-body methods, namely the Fermi
hypernetted chain/single-operator chain, and the auxiliary-field diffusion
Monte Carlo, was recently checked in benchmark calculations (Piarulli et al.,
2020) of the energy per particle of pure neutron matter.

##  3 Equation of state of hot nuclear matter

Figure 1: Free energy per particle ($F/A$) of SNM as a function of the nucleon
number density ($n$) for the interaction model considered in this work. The
square symbols represent the results of our microscopic BHF calculations,
whereas the continuous lines represent the energy per particle obtained using
the parametrization given by eq. (13). Temperatures (T) in the figure are in
MeV.

Figure 2: Same as in Figure 1, but for PNM.

For baryon densities greater than about $0.05$ ${\rm fm^{-3}}$, the stellar
matter can be described as an homogeneous fluid of neutrons and protons
(nuclear matter) with leptons (electrons and muons) to guarantee electric
charge neutrality. In this density region, the finite temperature BHF approach
outlined in the previous section is appropriate to describe the system. For
densities smaller than about $0.05$ ${\rm fm^{-3}}$, nucleons start to form
clusters and, as a consequence, the stellar matter becomes inhomogeneous333We
note however that this consideration strongly depends on the temperature of
the system: above some critical temperature ($\sim 18$ MeV), at all densities
clusters disappear and matter becomes homogeneous.. In these conditions, the
BHF approach does not provide a reliable description of the nuclear component
and another theoretical framework is therefore required to compute the EOS. We
return to this point below.

In the density range $0.05\leavevmode\nobreak\ \mathrm{fm^{-3}}\lesssim
n\lesssim 1.2\leavevmode\nobreak\ \mathrm{fm^{-3}}$, the free energy per
nucleon $\widetilde{F}=F/A$ of asymmetric nuclear matter can be calculated by
solving Eqs. (2)–(9) numerically for various values of the asymmetry parameter
($0\leq\beta\leq 1$) and temperatures ($0\lesssim T\lesssim
100\leavevmode\nobreak\ \mathrm{MeV}$). In order to evaluate the free-energy
contribution of asymmetric nuclear matter, we employ the so-called parabolic
approximation in the asymmetry parameter $\beta$ (Bombaci & Lombardo, 1991;
Bombaci et al., 1993, 1994; Wellenhofer et al., 2016):

$\widetilde{F}(n,\beta,T)=\widetilde{F}(n,\beta=0,T)+F_{\rm
sym}(n,T)\,\beta^{2}\,,$ (11)

where $F_{sym}(n,T)$ is the so-called nuclear symmetry free energy which, in
the $T=0$ limit, gives the nuclear symmetry energy. In the parabolic
approximation of Eq. (11), the symmetry free energy can therefore be obtained
as the difference between $\widetilde{F}$ of pure neutron matter (PNM,
$\beta=1$) and that of SNM ($\beta=0)$:

$F_{\rm sym}(n,T)=\widetilde{F}(n,\beta=1,T)-\widetilde{F}(n,\beta=0,T)\,.$
(12)

It has been shown that such approximation provides accurate results both at
zero (Bombaci & Lombardo, 1991) and finite temperature (Bombaci et al., 1993,
1994; Wellenhofer et al., 2016). Therefore, the two basic quantities appearing
in Eqs. (11) and (12) are the free energies per nucleon of SNM and of PNM as a
function of density and temperature. Their behaviors are shown in Fig. 1 and
Fig. 2 for SNM and PNM, respectively. The filled square symbols represent the
results of our microscopic BHF calculations described in the previous section,
whereas the continuous lines are the results of the fit obtained using the
following parametrization:

$\widetilde{F}(n,T)\equiv\frac{F}{A}(n,T)=(a_{0}+a_{1}t^{2}+a_{2}t^{3}+a_{3}t^{4})\
{\ln}\leavevmode\nobreak\
\Bigg{(}\frac{n}{\mathrm{fm}^{-3}}\Bigg{)}+(a_{4}+a_{5}t^{2})\
n+(a_{6}+a_{7}t^{2})\
n^{a_{8}}+\frac{1}{n}(a_{9}t^{2}+a_{10}t^{a_{11}})+a_{12}\;,$ (13)

where we define $t\equiv T/(100\leavevmode\nobreak\ \mathrm{MeV})$. The fit
parameters $a_{i}$ for SNM and PNM are reported in Table 1. The two-
dimensional fit in the variables ($n,T$) produces a very good $\chi^{2}/{\rm
datum}$ $\sim 1$ for both SNM and PNM.

Table 1: Value of the coefficients $a_{i}$ (with $i=0,1,...,12$) in the parametrization of the free energy per particle $F/A$ (see Eq. 13) for SNM (first row) and PNM (second row). The coefficients $a_{0}$, $a_{1}$, $a_{2}$, $a_{3}$ and $a_{12}$ are given in MeV, the coefficients $a_{4}$ and $a_{5}$ are given in MeV fm3, $a_{6}$ and $a_{7}$ in MeV fm${}^{3a_{8}}$, $a_{9}$ and $a_{10}$ in MeV fm-3, $a_{8}$ and $a_{11}$ are dimensionless quantities. The $\chi^{2}/{\rm datum}$ resulting from the fit is 1.05 for SNM and 1.1 for PNM. | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$ | $a_{6}$ | $a_{7}$ | $a_{8}$ | $a_{9}$ | $a_{10}$ | $a_{11}$ | $a_{12}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
SNM | -2.33 | 206.30 | -87.49 | 20.54 | -111.36 | -172.48 | 310.13 | 67.20 | 1.84 | 13.35 | -11.42 | 1.49 | -12.37
PNM | 3.205 | 158.61 | -54.64 | 12.11 | -28.17 | -123.14 | 415.96 | 47.55 | 1.91 | 19.43 | -18.37 | 1.66 | 17.00

Using these parametrizations of the free energy per nucleon we can easily
obtain analytic expressions for all the relevant thermodynamic quantities
defining the EOS of hot asymmetric nuclear matter. To this purpose we
introduce the free-energy density as
$\mathfrak{f}_{N}=F/V=n\leavevmode\nobreak\ (F/A)$, and then the neutron
($\mu_{n}$) and proton ($\mu_{p}$) chemical potentials can be obtained using:

$\mu_{n}\,(n,\beta,T)=\frac{\partial\mathfrak{f}_{N}}{\partial
n_{n}}\bigg{|}_{T,V,n_{p}}\,,\quad\quad\mu_{p}\,(n,\beta,T)=\frac{\partial\mathfrak{f}_{N}}{\partial
n_{p}}\bigg{|}_{T,V,n_{n}}\,,$ (14)

and the entropy density $s_{N}=S/V=n\leavevmode\nobreak\ (S/A)$ (in units of
Boltzmann constant $k_{B}$) can be obtained using:

$s_{N}(n,\beta,T)=-\frac{\partial\mathfrak{f}_{N}}{\partial
T}\bigg{|}_{V,n_{n},n_{p}}\,.$ (15)

The nuclear pressure is then given by:

$P_{N}(n,\beta,T)=\mu_{n}n_{n}+\mu_{p}n_{p}-\mathfrak{f}_{N}\,,$ (16)

and finally the energy density by:

$\varepsilon_{N}=\mathfrak{f}_{N}+Ts_{N}\,.$ (17)

Figure 3: Entropy per nucleon ($S/A$) as a function of the nucleon density $n$
at different temperatures, $T=20,50,70$ MeV (from bottom to top) for SNM (left
panel) and PNM (right panel). The continuous lines represent $S/A$ calculated
using Eq. (15) in which we employ the fit formula (Eq. (13)) for the free
energy. The circular symbols represent the entropy per nucleon, calculated
using Eq. (9) in the approximation of a mixture of two ideal Fermi gases of
quasi-particles in the BHF mean fields.

Figure 4: Energy per particle ($E/A$) as a function of the nucleon density $n$
at different temperatures $T=20,50,70$ MeV (from bottom to top) for SNM (left
panel) and PNM (right panel). The continuous lines represent
$E/A=\varepsilon_{N}/n$ calculated using Eq. (17) whereas the circular symbols
represent the results of the self-consistent BHF calculations, i.e., Eq. (8).

Before proceeding, we wish to comment on the procedure that we have outlined
above. It is well known that in the BHF approximation the nucleon chemical
potentials, calculated according to the thermodynamic definition in Eq. (14),
differ from the auxiliary chemical potentials in Eq. (5) used in the many-body
approach. As a consequence, the Hugenholtz-Van Hove theorem (Hugenholtz & Van
Hove, 1958) is not fulfilled. In the case of SNM at zero temperature, the
Hugenholtz-Van Hove theorem states that the energy per nucleon at the
saturation density must be equal to the Fermi energy of the quasi-particles. A
violation of this theorem is therefore a signal of the thermodynamic
inconsistency of the BHF approximation. This problem can be solved by going
beyond the BHF approximation of the single-particle energy by including some
extra diagrams of the hole-line expansion of the mass operator (see e.g.,
Hüfner & Mahaux, 1972) as discussed in Baldo et al. (1990); Zuo et al. (1999,
2006). In particular, it turns out that the inclusion of the so-called
rearrangement contribution $M_{2}$, which is a second-order diagram in the
$G$-matrix, and the renormalization contributions of third and fourth order in
the $G$-matrix, strongly improve the thermodynamic consistency of the
perturbative many-body approach. Inclusion of such diagrams is beyond the
scope of the present work. Here we employ a different scheme: we use the
expression of the free energy per particle obtained from our microscopic
calculations and provided by Eqs. (11)-(13), and afterwards we derive all the
other thermodynamic quantities from it; see Eqs. (14) - (17). Within this
scheme we obtain a thermodynamically consistent EOS. This strategy has been
successfully adopted in various works (Burgio & Schulze, 2010; Burgio et al.,
2011) which make use of the finite temperature extension of the BHF
approximation.

In order to compare the microscopic BHF calculations with respect to the
thermodynamically consistent EOS derived from the fit of the free energy, we
show in Fig. 3 the entropy per nucleon ($S/A$) as a function of the nucleon
density $n$ at different temperatures ($T=$ 20, 50, 70 MeV) for SNM (left
panel) and PNM (right panel). The continuous lines represent $S/A$ calculated
using Eq. (15) in which we employ the fit formula shown in Eq. (13) for the
free energy. The circular symbols represent the entropy per nucleon,
calculated using Eq. (9), in the approximation of a mixture of two ideal Fermi
gases of quasi-particles in the BHF mean fields $U_{\tau}(k;n,\beta,T)$. The
agreement between the two results is relatively good but deteriorates at the
highest temperatures and densities. In Fig. 4 we report the same type of
analysis as shown in Fig. 3 but for the energy per particle $E/A$. The
continuous lines represent $E/A=\varepsilon_{N}/n$ calculated using Eq. (17)
whereas the circular symbols represent the results of the self-consistent BHF
calculations of Eq. (8). The agreement between the two results is particularly
good at all the considered temperatures and densities.

We finally note that in the present paper we do not employ any approximation
in the calculation of the BHF single-particle neutron and proton potentials,
which are evaluated at finite temperature by solving Eqs. (2), (6), and (7)
consistently. A so-called frozen correlation approximation for $U(k,T)$,
namely the use of the values calculated at zero temperature, has often been
applied instead in other works (e.g., Burgio et al., 2011).

Additionally, we note that nuclear chiral interactions are usually employed up
to two to three times the nuclear saturation density $n_{0}$. As shown for
instance in Drischler et al. (2019), Piarulli et al. (2020), the uncertainty
estimates on many-body nuclear observables increase for increasing nuclear
density. However, to describe astrophysical systems like NSs or BNS mergers we
need to extend the density range of our many-body calculations above the
trusted region of chiral effective field theory. This scheme is more
consistent than the use of piecewise polytropic extensions of the EOS for
example, which are sometimes used in numerical general relativity simulations.

We now consider the low-density region, that is, baryon densities $n\leq 0.05\
{\rm fm}^{-3}$ where nucleon clusters appear. For $T\gtrsim
0.5\leavevmode\nobreak\ {\rm MeV}$, matter is in nuclear statistical
equilibrium (NSE), i.e., strong and electromagnetic interactions are in
equilibrium. Under these conditions, the possible nuclear species are also
fully determined by the baryon density, temperature, and asymmetry parameter.
The latter parameter is often expressed in terms of the total 444In this
context, ”total” refers to both free and bound protons in nuclei. proton
fraction, $Y_{p}$. In the present work in order to describe all the
thermodynamical quantities for $n\leq 0.05\ {\rm fm}^{-3}$, we consider the HS
EOS hempel10 in its publicly available tabular forms. This EOS was derived in
the framework of the RMF approach and models the appearance of nuclei through
an excluded volume approach (see Hempel & Schaffner-Bielich, 2010, for more
details). In particular, we consider two specific EOS tables derived from two
different RMF nuclear interaction parametrizations: the TM1 (Sugahara & Toki,
1994) and SFHo (Steiner et al., 2013c) ones. As detailed in these latter
studies, these EOSs take into account matter clusterization in the low-density
regions and therefore complement the EOS provided by our microscopic BHF
approach. Specifically, we performed a linear interpolation in all the
relevant thermodynamical quantities between the microscopic BHF and the RMF
calculations in the region $0.05\ {\rm fm}^{-3}\leq n\leq 0.08\ {\rm
fm}^{-3}$. We checked that the match between the EOSs provided by the two
different approaches is smooth enough to avoid the creation of fictitious
phase transitions. Despite its simplicity, such an interpolation preserves the
thermodynamical consistency of the two underlying EOSs. In addition, the
spherically symmetric CCSN simulations discussed in the following section
serve as additional verification of the construction of the complete EOS. Such
simulations are indeed very sensitive in a wide density range: from very low
density ($\sim 10^{-13}$ fm-3 ) up to $\sim 0.32$ fm-3. The closeness of the
match between the BHF EOS and the EOSs based on the RMF approach can therefore
be additionally tested by checking that the typical features of 1D CCSNe
simulations are correctly reproduced. We finally note that we have also
compared our results with others obtained using different EOSs for the low-
density region, such as the one reported by Togashi et al. (2017) which were
derived using a variational approach, or one of the Lattimer-Swesty (Lattimer
& Swesty, 1991) EOS family (specifically the one with compressibility modulus
$K_{\infty}=220$ MeV). We did not find appreciable differences at least for
our purposes.

##  4 Astrophysical applications

In astrophysical applications involving nuclear matter, the astrophysical
plasma is assumed to be charge neutral and composed of hadrons (nuclei,
nucleons), massive leptons (electron, muons, and their antiparticles), and
photons in thermodynamical equilibrium. Hadrons are usually considered in NSE.
Under these assumptions, the most general full EOS could be expressed as a
function of baryon density, temperature, and a couple of particle abundances,
such as for example the net (i.e., subtracted by the antiparticle
contribution) abundance of electrons, $Y_{e}\equiv
n_{e}/n=\left(n_{e^{-}}-n_{e^{+}}\right)/n$ and of muons, $Y_{\mu}\equiv
n_{\mu}/n=\left(n_{{\mu}^{-}}-n_{{\mu}^{+}}\right)/n$. The abundances $Y_{e}$
and $Y_{\mu}$ , together with that of protons ($Y_{p}$), must satisfy the
local charge neutrality. Depending on the matter conditions and on the
explored timescales, a gas of trapped neutrinos in weak equilibrium can also
be present.

In order to obtain a complete EOS suitable for astrophysical applications we
have to add the contributions of massive leptons ($L$), photons ($\gamma$),
and possibly neutrinos ($\nu$, in three different species) to the hadronic
entropy, energy density, and pressure, namely:
$S=S_{N}+S_{L}+S_{\nu}+S_{\gamma}$,
$\varepsilon=\varepsilon_{N}+\varepsilon_{L}+\varepsilon_{\nu}+\varepsilon_{\gamma}$,
and $P=P_{N}+P_{L}+P_{\nu}+P_{\gamma}$. The energy density
($\varepsilon_{i}$), pressure ($P_{i}$), and entropy ($S_{i}$) for the massive
leptons and neutrinos ($i=L,\ \nu$) have been computed using the expressions
for relativistic ideal Fermi gases at finite temperature with $m_{e}c^{2}$ =
0.511 MeV, $m_{\mu}c^{2}=105.658$ MeV and $m_{\nu}=0$ (Timmes & Arnett, 1999).
Additional requirements, including weak equilibrium conditions, are usually
needed to fix the relative amounts of electrons, muons, and neutrinos (see
below for a few concrete examples). The resulting number densities fix the
particle chemical potentials. Photon contributions are instead calculated
according to the Bose statistic and result in well-known analytic expressions.
In the remainder of this section, we show results related to the application
of our new EOS to various astrophysical systems. We first discuss the
calculation of the EOS for stellar $\beta$-stable nuclear matter; we then
determine the corresponding structure for nonrotating protoneutron stars, and
finally analyze the results of simulations of CCSNe. We also note that both
the $T=0$ and the finite temperature versions of the BL EOS have already been
used in numerical relativity simulations of BNSMs (Endrizzi et al., 2018;
Bernuzzi et al., 2020).

###  4.1 Equation of state for $\beta$-stable nuclear matter

As a first application of our new EOS, we derive a finite-temperature,
$\beta$-stable EOS suitable for describing hot neutron star matter. The latter
is required to model PNSs produced in successful CCSN explosions and the
possibly metastable remnant of a BNSM in the case where a black hole (BH) does
not form promptly at merger. In both cases, the large temperatures resulting
from the collapse or merger dynamics produce a copious amount of neutrinos,
which initially form a trapped neutrino gas and diffuse out over the diffusion
timescales ($\sim$ seconds). The fate of these hot neutron stars, as well as
their properties, depends on the progenitor system. Numerical simulations show
that in both scenarios the bulk of matter is characterized by a quasi-uniform,
low-entropy-per-baryon profile $\widetilde{S}=S/A\approx 1-3$ (in units of the
Boltzmann constant $k_{B}$), which decreases only over the cooling timescale
because of neutrino emission (Burrows & Lattimer, 1986; Prakash et al., 1997;
Pons et al., 1999; Fischer et al., 2010; Hüdepohl et al., 2010; Roberts et
al., 2012; Kastaun & Galeazzi, 2015; Perego et al., 2019). In the PNS case,
the total lepton fraction is initially rather high
($Y_{l_{e}}\equiv(n_{e}+n_{\nu_{e}})/n\sim 0.3-0.35$, where $n_{\nu_{i}}$ is
the net electron neutrino density) and the cooling process is associated with
a significant deleptonization towards cold neutrino-less weak equilibrium. In
the BNSM case, matter is characterized by an initially low $Y_{l_{e}}\lesssim
0.1$, and experiences a leptonization process because of decompression and
temperature increase.

To account for the different NS profiles and cooling phases, we model hot
nuclear matter by adopting an isoentropic EOS and fixing the value of the
entropy per baryon $S/A$ in the range $0\leq S/A\leq 3$. In the following
section, we also consider the case of a isothermal EOS although this case is
not particularly realistic. We consider both the case of neutrino-less and
neutrino trapped matter. In the latter case we fix the total lepton fractions
to $Y_{l_{e}}=0.33$ and$Y_{l_{\mu}}=0$.

In order to calculate the composition of $\beta$-stable nuclear matter we
solve the equations for chemical equilibrium at a given nuclear density $n$
and temperature $T$:

$\mu_{n}-\mu_{p}=\mu_{e}-\mu_{\nu_{e}}\;,\ \ \ \ \ \ \
\mu_{\mu}=\mu_{e}-\mu_{\nu_{e}}+\mu_{\nu_{\mu}}\;.$ (18)

In addition, charge neutrality requires

$n_{p}=n_{e}+n_{\mu}\,.$ (19)

In the case of neutrino-free matter we set: $\mu_{\nu_{e}}=\mu_{\nu_{\mu}}=0$.
Once the solutions of the above equations have been determined, all the other
thermodynamical quantities can be calculated according to standard relations
as discussed in the previous section.

In Fig. 5 we report the temperature profiles as a function of the nuclear
density for isoentropic, $\beta$-stable EOSs with constant entropy per baryon
$S/A=1$ (continuous lines) and $S/A=2$ (dashed lines), in the case of
neutrino-free and neutrino-trapped matter. As matter is compressed, the
temperature rises. As expected, larger temperatures are reached when
increasing the value of the entropy per baryon. For a fixed nuclear density,
the effect of neutrino-trapping reduces the value of the temperature compared
with the neutrino-free case. This effect is well understood noting that when
neutrinos are trapped, a new degree of freedom is introduced in the system
which therefore becomes more ”disordered” causing the total entropy to
increase. For a fixed value of $S/A$, and a given nuclear density $n$, a lower
temperature is therefore required.

Figure 5: Temperature profiles as a function of nuclear density $n$ for
isoentropic ($S/A=$ const) $\beta$-stable EOSs in the case of neutrino-free
(red lines) and neutrino trapped (blue lines) matter.

The composition of isoentropic ($S/A=2$) $\beta$-stable nuclear matter, that
is, the particle fractions $Y_{i}=n_{i}/n$ (with $i=n,p,e^{-},\mu^{-}$)
without (with) neutrino trapping, is shown in Fig. 6. In the same figure we
report the $T=0$ particle composition of the BL EOS for comparison (Bombaci &
Logoteta, 2018). The main effect of thermal contribution is to increase the
number of protons and to move the muon onset to smaller density. This behavior
is similar to that found by many other phenomenological and microscopic
approaches (Burgio et al., 2011). We note that in the case of matter including
trapped neutrinos (right panel in Fig. 6), the electron fraction rises and the
composition of matter becomes more symmetric in the neutron and proton
fractions than in the untrapped case. As a consequence, the EOS of neutrino-
trapped matter is expected to become softer (Prakash et al., 1997, see Bombaci
(1996) for a detailed discussion on the role of neutrino trapping and the
value of the neutron star maximum mass).

Figure 6:  Composition of dense matter for different physical regimes. In the
left (right) panel we show the composition of $\beta$-stable neutrino-free
(trapped) nuclear matter, i.e., the particle fractions $Y_{i}=n_{i}/n$ (with
$i=n,p,e^{-},\mu^{-}$), in the case of an isoentropic EOS with $S/A=2$. In the
right panel we fix $Y_{l_{e}}=0.33$. The dashed lines in the left panel
represent the composition of the $T=0$ BL EOS (Bombaci & Logoteta, 2018) which
is reported for comparison. See text for further details.

Figure 7 shows the EOSs (namely the relation between pressure and energy
density ($\varepsilon$)) for $\beta$-stable nuclear matter assuming several
thermodynamical stellar conditions in logarithmic (left panel) and linear
(right panel) scale. The two different scales are adopted in order to
highlight the low- and high-density parts of each EOS. We considered the
following cases: $T=0.1$ MeV, isoentropic with $S/A=1$ and $S/A=2$ without
neutrino trapping, and isoentropic with $S/A=1$ and $S/A=2$ with neutrino
trapping. We also report the $T=0$ EOS of Bombaci & Logoteta (2018) for
comparison. We note that for nuclear densities below $0.08$ $\rm fm^{-3}$, the
BL EOS is joined with the Sly4 EOS (Haensel & Pichon, 1994; Douchin & Haensel,
2001). This combination has been adopted in Endrizzi et al. (2018) in order to
simulate a BNSM of equal-mass neutron stars. In this work, as explained in the
previous section, for $n\leq 0.05\ {\rm fm}^{-3}$, we employ the TM1 model
(Sugahara & Toki, 1994) (and very similar results were obtained in the case of
the SFHo EOS). We therefore note that for each EOS reported in Fig. 7 we have
consistently employed the same thermodynamical conditions assumed for the core
as for low-matter density ($n\leq 0.05\ {\rm fm}^{-3}$). Additionally, we note
that in the case of trapped neutrinos we have enforced the weak equilibrium
conditions discussed above in the same density interval. An additional comment
about this point is in order. Neutrinos rapidly escape from the star in low-
density matter. This is usually schematized introducing the concept of a
neutrino-sphere inside which neutrinos are trapped. However, as discussed in
several works (e.g., Strobel et al., 1999; Fisher et al., 2009; Endrizzi et
al., 2020), the exact location of the ”radius” of the neutrino-sphere is
energy dependent. In the references quoted above this ”radius” was found in a
range of densities, namely $10^{-5}-10^{-3}\ {\rm fm}^{-3}$ , for relevant
neutrino energies. In the present work, fixing $S/A=1$ in $\beta$-stable
matter, we find that at a nuclear density of around $4.7\times 10^{-5}\ {\rm
fm}^{-3}$ neutrinos disappear. Increasing $S/A,$ this threshold drops
slightly; for instance setting $S/A=2$ we find that the threshold becomes
$1.4\times 10^{-5}\ {\rm fm}^{-3}$. We therefore locate the radius of the
neutrino-sphere according to the density at which neutrinos disappear. We
consider this criterion to deal with the low-density neutrino component as a
natural and reasonable one for our purposes.

Figure 7: Pressure vs. energy density ($\varepsilon$) relations for various
$\beta$-stable EOSs derived in different thermodynamical stellar conditions.
In the left (right) panel we use a logarithmic (linear) scale. The BL $T=0$
(black continuous line) taken from (Bombaci & Logoteta, 2018) is shown for
reference. We note that for the BL EOS for nuclear densities below $0.08$ $\rm
fm^{-3}$, the Sly4 EOS (Haensel & Pichon, 1994; Douchin & Haensel, 2001) has
been adopted (Endrizzi et al., 2018). The BLh with $T=0.1\ {\rm MeV}$ is
represented by the black dashed line. Red and blue curves refer to the
following cases: isoentropic EOS with $S/A=1$ and $S/A=2$ without neutrino
trapping (red dashed lines) and isoentropic EOS with $S/A=1$ and $S/A=2$ with
neutrino trapping (blue continuous and dashed lines). For the low-density part
of the BLh EOSs, we employ the TM1 model (Sugahara & Toki, 1994). Temperatures
(T) in the figure are in MeV. (See text for further details).

###  4.2 Neutron star structure

Figure 8: Mass–radius (left panel) and mass–central density ($n_{c}$) (right
panel) for the EOSs reported in Fig. 7. Temperatures (T) in the figure are in
MeV.

The structural properties of nonrotating neutron stars can be obtained by
integrating the equation for hydrostatic equilibrium in general relativity
(Tolman, 1934; Oppenheimer & Volkoff, 1939). We numerically integrated these
equations using the microscopic EOSs of $\beta$-stable nuclear matter
described in the previous sections. For each static neutron star
configuration, we additionally compute the total baryon number $N_{B}$, and we
define the baryonic mass (or “rest mass”) of the neutron star as
$M_{B}=m_{u}N_{B}$, where $m_{u}$ is the baryonic mass unit that we take equal
to $m_{u}=m(^{12}C)/12=1.6605\times 10^{-24}{\rm g}$. The total binding energy
of the star is therefore given by $E_{\rm bind}=(M_{B}-M_{G})c^{2}$ which
represents the total energy liberated during the birth of the neutron star.
The results are shown in Fig. 8, where we plot the gravitational mass $M_{G}$
(in unit of solar mass $M_{\odot}$) as a function of the surface stellar
radius $R$ (panel (a)), and the gravitational mass as a function of the
central density $n_{c}$ (panel (b)) for the considered EOS models. We note
that our zero-temperature EOS model (BL) is compatible with the largest
current measured neutron star mass, which is
$M_{G}=2.14^{+0.2}_{-0.18}\,M_{\odot}$ (Cromartie et al., 2019) for the
neutron stars in PSR J0740+6620. Various structural properties of the maximum
mass configuration for the considered EOS models are listed in Table 2. We
note that thermal effects lead to a monotonic behavior (Bombaci et al., 1995)
of the gravitational maximum mass, both in the case of isoentropic and
isothermal EOSs. Neutrino trapping instead makes the EOS softer and thus leads
to a decrease in the values of both the gravitational and baryonic maximum
mass (Bombaci, 1996; Prakash et al., 1997). This can be seen when comparing
the results in Table 2 for the isoentropic EOS in the neutrino-free and
neutrino-trapped (last three rows in Table 2) cases. This effect was already
apparent in the discussion of the composition of $\beta$-stable matter in the
neutrino-free and neutrino-trapped cases (see Fig. 6). The electron chemical
potential $\mu_{e}$ in neutrino-trapped matter is indeed increased by the
value of the neutrino chemical potential $\mu_{\nu_{e}}$ (see Eq. 18). In this
way, a larger fraction of electrons and thus protons (due to charge
neutrality) are present at a given nuclear density. Consequently, the EOS
becomes softer, being more proton rich. However, we note that the presence of
other degrees of freedom in neutrino-trapped matter, like for instance
hyperons, may change the previous conclusion (Bombaci, 1996; Prakash et al.,
1997). A final remark is that a consistent treatment of neutrinos in the low-
density region produces an increase in the minimum neutron-star mass. Looking
at Fig. 7 it is clear that for $1\ {\rm MeV}\ {\rm
fm}^{-3}\lesssim\varepsilon\lesssim\ 10^{3}\ {\rm MeV}\ {\rm fm}^{-3}$ (which
corresponds to $10^{-5}\ {\rm fm}^{-3}\lesssim n\lesssim 0.1\ {\rm fm}^{-3}$)
the EOSs with trapped neutrinos (blue lines) are stiffer than the ones
obtained in the neutrino-less case. This is in agreement with the findings of
other works (Burgio & Schulze, 2010).

Table 2: Maximum mass configuration properties for various thermodynamical
stellar conditions considered in this work.

EOS model |  $M_{G}$ ($M_{\odot}$) |  $R$ (km) |  $n_{c}$ (fm-3) |  $M_{B}$ ($M_{\odot}$)
---|---|---|---|---
BL $T=0$ | 2.080 | 10.26 | 1.156 | 2.456
BLh $T=0.1$ | 2.070 | 10.18 | 1.175 | 2.457
BLh $T=5$ | 2.072 | 10.00 | 1.162 | 2.457
BLh $T=10$ | 2.073 | 10.05 | 1.158 | 2.457
BLh $T=30$ | 2.090 | 10.31 | 1.141 | 2.449
BLh $T=50$ | 2.112 | 10.91 | 1.084 | 2.406
BLh $S/A=1$ | 2.076 | 10.37 | 1.141 | 2.446
BLh $S/A=2$ | 2.082 | 10.81 | 1.098 | 2.390
BLh $S/A=3$ | 2.088 | 11.97 | 0.992 | 2.296
BLh $S/A=1$ $\nu$ | 2.021 | 10.39 | 1.190 | 2.323
BLh $S/A=2$ $\nu$ | 2.034 | 10.92 | 1.129 | 2.289
BLh $S/A=3$ $\nu$ | 2.051 | 12.07 | 1.018 | 2.223
555Stellar gravitational maximum mass $M_{G}$, corresponding radius $R$,
central baryon number density $n_{c}$, and baryonic maximum mass $M_{B}$.
Stellar masses are given in units of solar mass $M_{\odot}=1.989\times
10^{33}$ g.

To test the zero temperature limit of our finite temperature EOS model, we
calculated an isothermal NS sequence with $T=0.1$ MeV. As one can see,
comparing the results in the first two rows of Table 2, we recover the
properties of the maximum mass configuration predicted by the $T=0$ BL EOS.
The differences in the radius (see Fig. 8) for stars with masses of less than
about $1\,M_{\odot}$, between the BL EOS and the BLh at $T=0.1$ MeV, are due
to the different EOS models used at low densities
($n_{B}<0.05\leavevmode\nobreak\ {\rm fm^{-3}}$). We note in addition that in
Bombaci & Logoteta (2018), in order to construct a parametrization of the
energy per nucleon for nuclear $E/A$ densities $n\geq 0.08$ fm-3, only the
potential energy part of microscopic calculations was included in the fit. At
zero temperature, the kinetic contribution to $E/A$ is indeed analytic (Eqs.
(6) and (7) in Bombaci & Logoteta (2018)). In the present work, in which we
have extended the microscopic results of Bombaci & Logoteta (2018) to finite
temperature, we have instead performed a fit of the full free energy per
nucleon $F/A$ because in this case the kinetic term is also nonanalytic. It is
therefore important that the two $\beta$-stable EOSs BL and BLh ($T=0.1$ MeV)
provide almost indistinguishable results (see right panel of Fig. 7) in the
range $n\geq 0.08$ fm-3.

Figure 9:  Gravitational mass $M_{G}$ as a function of the baryonic mass
$M_{B}$ for our EOS model. The red dashed line represents the protoneutron
star (PNS) sequence calculated using the isoentropic EOS with $S/A=2$ and
trapped neutrinos with $Y_{l_{e}}=0.33$. The black continuous line represents
the cold deleptonized NS sequence calculated using the
$T=0.1\leavevmode\nobreak\ {\rm MeV}$ EOS for neutrino-free matter. The end
point on each curve represents the corresponding maximum mass configuration.
The values of these quantities for the present BLh EOS model are reported in
Table 2. The asterisk on the black continuous line represents the effective
maximum mass $M_{\rm G\,max}^{*}$ for the cold deleptonized NS sequence (see
text for more details).

###  4.3 Early evolution of a newborn neutron star

It has been shown (Bombaci, 1996) that the concept of neutron star maximum
mass has to be re-examined when the formation and early evolution dynamics of
a newborn PNS are taken into account. The early evolution of a PNS is driven
by thermal and neutrino trapping effects on the EOS. The main features of this
process can be schematically investigated considering the following two
snapshots of the evolution process:

  * (i)

the hot PNS at a time $t\sim 3\,\mathrm{s}$ after core bounce ($t=0$),
described by the isoentropic EOS with trapped neutrinos;

  * (ii)

the cold and deleptonized NS at $t\sim 30\,\mathrm{s}$ (neutrino diffusion
time), described by the cold and neutrino-free EOS.

As most of the matter accretion on the forming NS happens in the very early
stages after birth ($t<3\,\mathrm{s}$) (Chevalier, 1989), the neutron star
baryonic mass $M_{\rm B}$ stays almost constant during the evolution between
these two configurations. The evolution of a PNS can therefore be analyzed in
the $M_{\rm G}$–$M_{\rm B}$ plane. To this purpose we show in Fig. 9 the
gravitational mass as a function of the baryonic mass for the following
stellar sequences: (i) PNSs (red dashed line), that is, isoentropic EOS with
$S/A=2$ and trapped neutrinos with $Y_{l_{e}}=0.33$; and (ii) cold
deleptonized NSs (continuous black line) with $T=0.1\leavevmode\nobreak\ {\rm
MeV}$ and neutrino-free matter. The end point on each curve represents the
corresponding maximum mass configuration, which we denote as $M_{\rm
G\,max}^{(i)}$ and $M_{\rm B\,max}^{(i)}$ for the PNS sequence, and with
$M_{\rm G\,max}^{(f)}$ and $M_{\rm B\,max}^{(f)}$ for the final cold
deleptonized NS sequence. The values of these quantities for the present BLh
EOS model are reported in Table 2. From our results in Fig. 9 one can see that
a PNS born with $M_{\rm B}^{(i)}\leq M_{\rm B\,max}^{(i)}$ will evolve to a
cold deleptonized NS with a gravitational mass $M_{\rm G}^{(f)}$ and the same
baryonic mass $M_{\rm B}^{(i)}$ of the initial PNS configuration. The binding
energy of the star will increase by $\Delta B=(M_{\rm G}^{(i)}-M_{\rm
G}^{(f)})\,c^{2}$. This energy will be released mainly through neutrino
emission. In the specific case where $M_{\rm B}^{(i)}=M_{\rm B\,max}^{(i)}$
and for the case considered in Fig. 9, we find $M_{\rm
G}^{(f)}=1.958\,M_{\odot}$ and $\Delta B=1.358\times 10^{53}\mathrm{erg}$. We
now consider a PNS with $M_{\rm B\,max}^{(i)}<M_{\rm B}^{(i)}\leq M_{\rm
B\,max}^{(f)}$ . In this case the PNS cannot be supported by the matter
pressure against gravitational collapse because its baryonic mass $M_{\rm
B}^{(i)}$ is greater than the maximum possible baryonic mass $M_{\rm
B\,max}^{(i)}$ for the initial configuration. Thus, the collapsing stellar
core will collapse to a BH after reaching supranuclear densities. In the
classical analysis, à la Oppenheimer–Volkoff, where the dynamical evolution of
the PNS is not taken into account, stars in this baryonic mass range are
considered to have a stable equilibrium configuration in the cold deleptonized
stellar sequence. Thus, the gravitational mass $M_{\rm G\,max}^{*}\equiv
M_{\rm G}^{(f)}(M_{\rm B\,max}^{(i)})$ of the star, corresponding to the
evolution of the maximum mass PNS configuration, plays the role of an
effective maximum mass for the cold deleptonized NS sequence (Bombaci, 1996).
We stress that this is a feasible mechanism to produce low-mass black holes
($M_{\rm BH}\sim 2\,M_{\odot}$) and could have far-reaching consequences for
interpreting the gravitational wave event GW190814 (Abbott et al., 2020b) as a
BH–BH merger. Finally, when $M_{\rm B}^{(i)}>M_{\rm B\,max}^{(f)}$ the stellar
core will collapse to a BH.

Consider now a NS in a binary stellar system in which the companion star is a
normal star. If a common envelope is formed during the evolution of the binary
system, the NS can accrete matter from its companion star with a certain
accretion rate $\dot{M}_{\rm B}$. After a sufficiently long time the NS is
able to increase its baryonic mass above the value $M_{\rm B\,max}^{(i)}$ and
begins to populate the portion of the final NS sequence with $M_{\rm
B\,max}^{(i)}<M_{\rm B}\leq M_{\rm B\,max}^{(f)}$ , eventually reaching and
then overcoming the Oppenheimer–Volkoff maximum mass configuration $M_{\rm
B\,max}^{(f)}$ . Thus, as pointed out for the first time in (Bombaci, 1996),
in the baryonic mass range $[M_{\rm B\,max}^{(i)}\;,M_{\rm B\,max}^{(f)}]$ one
can have both NSs and BHs. For the case depicted in Fig. 9, this range is
[2.289 , 2.457] $M_{\odot}$ .

###  4.4 Application to spherically symmetric core-collapse supernovae

Figure 10: Left: Temporal evolution of the shock radius (solid) and of the PNS
radius (dashed) for each of the four spherically symmetric CCSN models. Models
differ by the nuclear EOSs, including two versions of the BLh EOS and two RMF
EOSs, namely HS(TM1) and SFHo. The two BLh EOSs use the HS(TM1) and the SFHo
EOS for the low-density regime ($n<n_{\rm th}$). We note the logarithmic scale
on the time axis. Right: Same as on the left panel, but for the central
density (solid) and temperature (dashed). The cyan (blue) shaded area
represents matter below (above) the threshold density $n_{\rm th}$ (used to
join the BLh and the HS(TM1) EOS) for the BLh+TM1 model.

Figure 11: Left: Same as in Figure 10, but for the electron neutrino (solid)
and antineutrino (dashed) luminosities measured at the edge of the
computational domain. Right: Same as in the left panel, but for the electron
neutrino (solid) and antineutrino (dashed) mean energies.

To test our finite temperature EOS in an astrophysical dynamical scenario, we
set up simulations of a CCSN in spherical symmetry. For recent reviews of the
physics and the status of CCSN modeling we refer to Burrows (2013), Janka et
al. (2016). In the following, we summarize only the most significant aspects
necessary to understand results from our simulations.

A CCSN marks the end of the life of a massive star, with a zero-age main
sequence mass above $\sim 8M_{\odot}$, up to the pair instability limit.
During the so-called silicon-burning phase, in the centre of the star silicon
is converted into iron group elements in NSE conditions. The forming iron core
is mainly sustained by the pressure of degenerate electrons. However, once the
Chandrasekhar mass limit has been reached, the gravitational collapse becomes
unavoidable and the core collapses under its own gravity. At the onset of the
collapse, matter in the center has already developed a significant deviation
from isospin symmetry ($Y_{e}\lesssim 0.48$) and typical densities and
temperatures are $\rho_{c}\sim$ $(3-8)$ $\times$ $10^{9}{\rm
g\leavevmode\nobreak\ cm^{-3}}$ and $T_{c}\gtrsim 0.5\leavevmode\nobreak\ {\rm
MeV}$ , respectively, with an intrinsic variability due to the progenitor
initial mass, metallicity, and possibly binarity. The core collapses on the
dynamical timescale of $\tau\sim\left(G\rho_{c}\right)^{-1/2}\sim
100\leavevmode\nobreak\ {\rm ms}$ and the collapse is further accelerated by
the core deleptonization, through electron captures on nuclei and protons, and
the subsequent escape of the produced neutrinos. While neutrinos initially
escape freely, only when $\rho_{c}\sim 10^{12}{\rm g\leavevmode\nobreak\
cm^{-3}}$ does the typical neutrino mean free path decrease below the size of
the inner core and neutrinos start to be trapped. The collapse proceeds until
$\rho_{c}\sim\rho_{0}$. The subsequent EOS stiffening halts the collapse and
produces an outgoing shock wave. The moment when the maximum density is
reached is called core-bounce. At core bounce, the electron fraction at the
center of the core is $Y_{e}\lesssim 0.3$, while the fraction of trapped
electron neutrino is $Y_{\nu_{e}}\lesssim 0.1$. Temperatures increase up to a
few tens of MeV and the shock photodissociates the iron core into free
neutrons and protons. Neutrinos of all flavors are copiously produced and
emitted with typical neutrino luminosities of the order of a few times
$10^{52}{\rm erg\leavevmode\nobreak\ s^{-1}}$. After a few tens of
milliseconds, the shock loses its energy because of photodissociation and
neutrino emission, and then stalls and converts into an accretion shock,
inside which a PNS forms. The neutrino diffusion from the optically thick PNS,
and the persistent deleptonization and compression of the accreting matter
guarantee high neutrino luminosities over a timescale of one second. The re-
absorption of electron flavor (anti)neutrinos mainly emitted at the surface of
the PNS produces a region of net energy deposition just behind the shock (the
so-called gain region). This energy deposition is thought to be the cause of
the shock revival that eventually explodes the star. Detailed models in
spherical symmetries including all the relevant physics fail to revive the
shock through neutrino heating. Multi-dimensional simulations proved that
convection and hydrodynamical instabilities play a crucial role in helping the
explosion to occur.

Our hydrodynamic simulations are performed with the Lagrangian general-
relativistic, adaptive-grid code Agile (Liebendörfer et al., 2001). We employ
the isotropic diffusion source approximation (IDSA) for the transport of
electron (anti)neutrinos (Liebendörfer et al., 2009). The neutrino energy is
discretized using $20$ geometrically spaced energy bins in the range
$\left[3,300\right]$ MeV. The neutrino reactions included in the IDSA are
electron captures on nuclei, electron and positron captures on free nucleons,
and scattering off free nucleons and nuclei (Bruenn, 1985; Mezzacappa &
Bruenn, 1993). These represent the minimal set of the most relevant weak
processes during the collapse and in the early post-bounce phase. The emission
of $\mu$ and $\tau$ (anti)neutrinos in the post-bounce phase is handled by a
gray leakage scheme (Cabezón et al., 2018). This scheme has been gauged on
detailed Boltzmann transport models, and includes pair annihilation,
nucleon–nucleon bremsstrahlung, and scattering off nucleons. With this setup,
we simulate the collapse, bounce, and post-bounce phases (roughly the first
700ms after bounce) for a 20$M_{\odot}$ ZAMS, solar metallicity progenitor
star (Woosley & Heger, 2007). We have chosen this progenitor because it has
been extensively explored in several recent comparison papers in spherical
symmetry and in multi-dimensions (O’ Connor et al., 2018; Pan et al., 2019;
Schneider et al., 2019). For this specific progenitor, at the onset of the
collapse, $\rho_{c}\approx 3.3\times 10^{9}{\rm g\leavevmode\nobreak\
cm^{-3}}$, $T_{c}\approx 0.67{\rm MeV,}$ and $Y_{e,c}\approx 0.44$. In all our
models we use 103 radial zones which includes the innermost $\sim
2.1M_{\odot}$. This corresponds to an initial radius of $R\approx 5\times
10^{8}{\rm cm}$. We finally stress that the code that we are using is
identical to the publicly available version of AGILE-
IDSA666https://astro.physik.unibas.ch/people/matthias-
liebendoerfer/download.html.

To close the set of hydrodynamics equations and to compute the necessary
neutrino emissivities, absorptivities, and opacities, we couple a finite-
temperature, composition-dependent nuclear EOS to the CCSN model. We consider
four different EOSs: two versions of the complete BLh EOS and two RMF EOSs,
namely the SFHo and HD(TM1) EOS. The two versions of the BLh EOS differ in the
EOS used for the low-density part ($n_{B}<0.05\leavevmode\nobreak\ {\rm
fm^{-3}}$) joined to the high-density BLh EOS according to the procedure
described in Section 3. All these EOSs include electron, positrons, and
photons as outlined at the beginning of this section. In particular, we used
precisely the same SFHo and HS(TM1) EOSs. Thus, our choice allows us to
directly address the impact of the transition to the low-density EOS, as well
as the difference between our EOS and other EOSs commonly used in the
literature. Trapped neutrinos are consistently modeled by the IDSA, while
muons are neglected. We stress that the use of a CCSN model to test the
transition is challenging, because all matter inside the domain is initially
at low densities and most of it has dynamically transited to high densities by
the end of the simulation.

All the different models run without problems through the different CCSN
phases (collapse, bounce, after-bounce). Core bounce happens
$\sim$378.6(384.6) ms after the beginning of the simulation for the SFHo and
TM1 (for the SFHo+BLh and TM1+BLh) models. In Figures 10 and 11 we present the
outcome of our set of simulations. In particular, in the left panel of Figure
10 the temporal evolution of the shock radius and the PNS radius are depicted.
The latter is defined as the radius where the rest mass density is equal to
$10^{11}{\rm g\leavevmode\nobreak\ cm^{-3}}$. In the right panel of the same
figure, we present the evolution of the central density and temperature. In
Figure 11 the neutrino luminosities and mean energies are reported from the
latest phase of the collapse up to the end of the simulation.

The CCSN models employing the BLh high-density EOS show all the features
expected from a spherically symmetric CCSN. The PNS radius, shock radius,
neutrino luminosities, and mean energies follow the expected trends, with all
the models being close to one another throughout the simulated time. When the
central density reaches values above nuclear saturation density, the core
bounces, and an outgoing shock wave forms. The shock expands up to $\sim$135km
and then recedes and converts into an accretion shock. A hot mantle is built
below the shock, and the forming PNS contracts inside it due to the emission
of neutrinos.

Matter in the center compresses, steadily increasing the central density. Due
to compression, the core temperature also increases. As the neutrino diffusion
time from the PNS is of the order of seconds (much larger than the simulated
time), the lepton fraction, the electron fraction, and the entropy in the
center remain almost constant during the simulation ($Y_{l}\approx 0.32-0.34$
, $Y_{e}\approx 0.28-0.30$ and $S/A\approx 1.25-1.35$). Matter in the center
is highly degenerate and the increase in temperature is relatively small
($T_{c}\lesssim 20-30\leavevmode\nobreak\ {\rm MeV}$, even at the end of the
simulation) compared to the corresponding Fermi temperature. Among the
different EOSs, HS(TM1) is the stiffest and its central densities are
significantly lower than for the other EOSs. The largest densities are reached
by the BLh-based EOSs, with very minor differences between its two different
implementations. The BLh models also present the largest central temperature.
For this quantity, the relative difference between BLh+SFHo and BLh+TM1 varies
slightly with time, with the latter being very close to the former close to
bounce and towards the end of the simulation. Due to the lower densities and
temperatures reached by the TM1 model, its shock and PNS radii are
systematically larger than for other EOS models, while the neutrino
luminosities and mean energies are systematically smaller. Apart from this
expected difference, we note how the SFHo+BLh and TM1+BLh stay close to each
other during the entire simulation, including during the $\nu_{e}$ burst peak
and the corresponding rapid shock expansion phase, between 1 and 10 ms after
bounce. Finally, in the lower part of the right panel of Figure 10, by the end
of the simulation, for models using the BLh-based EOS, most of the matter in
the computational domain is described by the high-density part of the EOS.

##  5 Summary

In this work, we extend the zero-temperature microscopic EOS of dense
asymmetric and $\beta$-stable nuclear matter of Bombaci & Logoteta (2018) to
finite temperature. For consistency, we used the same many-body method as well
as the same two- and three-body nuclear forces of Bombaci & Logoteta (2018).
For the low-ensity part of our EOS (namely below $0.05\ {\rm fm^{-3}}$), we
connected our microscopic EOS with a EOS derived in the RMF approach (Hempel &
Schaffner-Bielich, 2010). In order to test different choices, we used two
publicly available EOS tables obtained from two different parametrizations of
the nuclear RMF Lagrangian: the TM1 (Sugahara & Toki, 1994) and SFHo (Steiner
et al., 2013c) ones.

We considered three applications of this new EOS. The first one concerned the
calculation of the isoentropic $\beta$-stable EOS considering both the case of
$\nu$-free and $\nu$-trapped matter. We note that the effects of neutrino
trapping has been analyzed in the whole range of density of interest for hot
neutron stars, therefore including the low-density EOS part. In our second
application, we used the finite temperature $\beta$-stable EOSs to describe
the structure of nonrotating protoneutron stars. We discuss the thermal
effects on the maximum and minimum neutron star masses obtained according to
the various physical scenarios analyzed in this paper. In our last
application, we employed our new EOS to perform a dynamical simulation of
spherically symmetric CCSNe. This last application turned out to be
particularly important because it is very sensitive to the density range where
we match microscopic calculations with the RMF results. Our findings are in
agreement with those present in the literature on this kind of simulation,
supporting the plausibility of the EOS construction proposed in this paper.

We note in addition that the zero-temperature neutron star maximum mass
predicted by our new EOS is consistent with present measured neutron star
masses and particularly with the mass $M=2.14^{+0.20}_{-0.18}\,M_{\odot}$ of
the neutron stars in PSR J0740+6620.

The consistent thermodynamical approach adopted in this work is relevant for
multidimensional CCSNe and BNSM simulations (Shen et al., 1998; Hempel &
Schaffner-Bielich, 2010; Steiner et al., 2013c; Oertel et al., 2017; Togashi
et al., 2017). In the near future we indeed plan to apply this new EOS to
BNSMs as well as to 3D CCSN simulations.

We also plan to develop a more consistent treatment of the low-density part of
our EOS namely below $0.05{\rm fm^{-3}}$. A possibility may consist in finding
a skyrme parametrization able to fit microscopic calculations around
saturation density and to join the high-density part of the EOS based on fully
microscopic calculations with the latter.

We finally note that due to the very high central densities expected in the
core of NSs, various “exotic” constituents like hyperons (Glendenning, 1985;
Vidaña et al., 2011; Chatterjee & Vidaña, 2016; Haidenbauer et al., 2017) or a
deconfined quark phase (Glendenning, 1996; Bombaci et al., 2008; Logoteta et
al., 2012; Bombaci & Logoteta, 2013; Bombaci et al., 2016) cannot be excluded.
Thermal effects may favor the presence of these additional degrees of freedom
in neutron star matter. Work in this direction is ongoing. Specifically, we
are currently trying to extend our present calculations to the case of finite-
temperature hyperonic matter, adopting nucleon-hyperon (Haidenbauer et al.,
2013) and hyperon-hyperon (Haidenbauer et al., 2016) interactions derived in
the framework of ChPT, and to include hyperonic three-body forces like those
in the zero-temperature case discussed by Logoteta et al. (2019).

## Acknowledgement

This work has been partially supported by “PHAROS”, COST Action MP1304. AP
acknowledges the usage of computer resources under a CINECA-INFN (allocation
INF20_teongrav).

## References

  * Abbott et al. (2017) Abbott, B. P. et al. (LIGO Scientific Collaboration and Virgo Collaboration) 2017, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2018) Abbott, B. P. et al. (LIGO Scientific Collaboration and Virgo Collaboration) 2018, Phys. Rev. Lett., 121, 161101
  * Abbott et al. (2020a) Abbott, B. P. et al. (LIGO Scientific Collaboration and Virgo Collaboration) 2020, ApJ Lett., 892, L3
  * Abbott et al. (2020b) Abbott, B. P. et al. (LIGO Scientific Collaboration and Virgo Collaboration) 2020, ApJ Lett., 896, L44
  * Akmal et al. (1998) Akmal, A., Pandharipande V. R., & Ravenhall D. G., 1998, Phys. Rev. C, 58, 1804
  * Antoniadis et al. (2013) Antoniadis, J. et al. 2013 Science, 340, 1233232
  * Baiotti & Rezzolla (2017) Baiotti, L., & Rezzolla, L. 2017, Rep. Progr. in Phys., 80, 9
  * Baldo et al. (1997) Baldo M., Bombaci I., & Burgio G. F. 1997, A&A , 328, 274
  * Baldo et al. (1991) Baldo M., Bombaci I., Ferreira L. S., Giansiracusa G., & Lombardo U. 1991, Phys. Rev. C, 43, 2605
  * Baldo et al. (1990) Baldo, M., Bombaci I., Giansiracusa G., Lombardo U., Mahaux C., & Sartor R. 1990 Phys. Rev. C, 41, 1748
  * Baldo & Ferreira (1999) Baldo, M., & Ferreira, L.S., Phys. Rev. C, 1999, 59, 682
  * Baldo et al. (2000) Baldo M., Giansiracusa G., Lombardo U., & Song H. Q. 2000, Phys. Lett. B, 473, 1
  * Bauswein & Janka (2012) Bauswein, A., & Janka, H.-T. 2012, Phys. Rev. Lett., 108, 011101
  * Bauswein et al. (2010) Bauswein, A., Janka, H.-T. & Oechslin, R. 2010, Phys. Rev. D, 82, 084043
  * Bauswein et al. (2019) Bauswein, A., Bastian, N.-U. F., Blaschke, D. B., et al. 2019, Phys. Rev. Lett., 122, 061102
  * Bernuzzi et al. (2015) Bernuzzi, S., Dietrich, T, & Nagar, A. 2015, Phys. Rev. Lett., 115, 091101
  * Bernuzzi et al. (2016) Bernuzzi, S., Radice, D., Ott, C. D., et al. 2016, Phys. Rev. D, 94, 024023
  * Bernuzzi et al. (2020) Bernuzzi, S., Breschi, M., Daszuta, B., Endrizzi, A., Logoteta, D., Nedora, V., Perego, A., Radice, D., Schianchi, F., Zappa, F., Bombaci, I., & Ortiz, N. 2020, MNRAS, 497, 1488
  * Bethe (1965) Bethe, H. A. 1965, Phys. Rev., 138, 804B
  * Bethe et al. (1963) Bethe, H. A., Brandow, B. H., & Petschek, A. G. 1963, Phys. Rev., 129, 225
  * Binder et al. (2016) Binder, S. et al. (LENPIC Collaboration) 2016, Phys. Rev. C, 93, 044002
  * Bloch & de Dominicis (1958) Bloch, C., & De Dominicis, C. 1958, Nucl. Phys. A, 7, 459
  * Bloch & de Dominicis (1959a) Bloch, C., & De Dominicis, C. 1959, Nucl. Phys. A, 10, 181
  * Bloch & de Dominicis (1959b) Bloch, C., & De Dominicis, C. 1959, Nucl. Phys. A, 10, 509
  * Bombaci (1996) Bombaci, I. 1996, A&A, 305, 817
  * Bombaci et al. (1995) Bombaci, I., Prakash, M., Prakash, M., Ellis, P.J., Lattimer, J.M., & Brown, G.E. 1995, Nucl. Phys. A 583, 623
  * Bombaci et al. (2008) Bombaci, I., Panda, P. K., Providência, C., & Vidaña, I. 2008, Phys. Rev. D, 77, 083002
  * Bombaci & Logoteta (2013) Bombaci, I., & Logoteta, D. 2013, MNRAS Lett. 433, L79
  * Bombaci et al. (2016) Bombaci, I., Logoteta, D., Vidaña, I., & Providência, C. 2016, Eur. Phys. J. A, 52, 58
  * Bombaci & Logoteta (2018) Bombaci, I., & Logoteta, D. 2018, A&A, 609, A128
  * Bombaci & Lombardo (1991) Bombaci, I., & Lombardo, U. 1991, Phys. Rev. C, 44, 1892
  * Bombaci et al. (1993) Bombaci, I., Kuo, T. S., & Lombardo, U. 1993, Phys. Lett., B 311, 9
  * Bombaci et al. (1994) Bombaci, I., Kuo, T. S., & Lombardo, U. 1994, Phys. Rep., 242, 165
  * Bombaci et al. (2006) Bombaci, I., Polls, A., Ramos, A., Rios, A., & Vidaña, I. 2006, Phys. Lett. B, 632, 638
  * Bruenn (1985) Bruenn, S. W. 1985, ApJS, 58, 771
  * Burgio & Schulze (2010) Burgio, G. F., & Schulze, H.-J. 2010, A&A 518, A17
  * Burgio et al. (2011) Burgio, G. F., Schulze, H.-J., & Li, A. 2011, Phys. Rev. C, 83, 025804
  * Burrows (2013) Burrows, A. 2013, Reviews of Modern Physics, 85, 245
  * Burrows & Lattimer (1986) Burrows, A. & Lattimer, J. M. 1986, ApJ, 307, 178
  * Cabezón et al. (2018) Cabezón, R. M., Pan, K.-C., Liebendörfer, M., et al. 2018, A&A, 619, A118
  * Carbone et al. (2018) Carbone, A., Polls, A., & Rios, A. 2018, Phys. Rev. C, 98, 025804
  * Carbone & Schwenk (2019) Carbone, A., & Schwenk, A. 2019, Phys. Rev. C, 100, 025805
  * Chamel et al. (2011) Chamel, N., Fantina, A. F., Paearson, J. M., & Goriely, S. 2011, Phys. Rev. C, 84, 062802(R)
  * Chatterjee & Vidaña (2016) Chatterjee, D. & Vidaña, I. 2016, Eur. Phys. J. A, 52, 29
  * Chevalier (1989) Chevalier, R.A. 1989, ApJ, 346, 847
  * Coester et al. (1970) Coester F., Cohen S., Day B., & Vincent C. M. 1970, Phys. Rev. C, 1, 769
  * Cromartie et al. (2019) Cromartie, H. T. et al. 2019, Nature Astronomy 10.1038
  * Day (1967) Day, B. D. 1967, Rev. Mod. Phys., 39, 719
  * Day (1981) Day B. D. 1981, Phys. Rev. Lett., 47, 226
  * Demorest et al. (2010) Demorest P., Pennucci T., Ransom S., Roberts M., Hessels J. 2010, Nature, 467, 1081
  * Douchin & Haensel (2001) Douchin, F., & Haensel, P. 2001, A&A 380, 151
  * Drischler et al. (2016) Drischler, C., Hebeler, K., & Schwenk, A. 2016, Phys. Rev. C, 93, 054314
  * Drischler et al. (2019) Drischler, C., Hebeler, K., & Schwenk, A. 2019, Phys. Rev. Lett., 122, 042501
  * Endrizzi et al. (2018) Endrizzi, A., Logoteta, D., Giacomazzo, B., Bombaci, I., Kastaun, W., & Ciolfi, R. 2018, Phys. Rev. D, 98, 043015
  * Endrizzi et al. (2020) Endrizzi, A., Perego, A., Fabbri, F. M., Branca, L., Radice, D., Bernuzzi, S., Giacomazzo, B., Pederiva, F. & Lovato, A. 2020, Eur. Phys. J. A, 56, 15
  * Fattoyev et al. (2010) Fattoyev, F. J., Horowitz, C. J., Piekarewicz, J., & Shen, G. 2010, Phys. Rev. C, 82, 055803
  * Fiorilla et al. (2012) Fiorilla, S., Kaiser, N. & Weise, W. 2012, Nucl. Phys. A, 880, 65
  * Fisher et al. (2009) Fischer, T., Whitehouse, S. C., Mezzacappa, A., Thielemann, F.-K., & Liebendörfer, M. 2009, A&A, 499, 1
  * Fischer et al. (2010) Fischer, T., Whitehouse, S. C., Mezzacappa, A., et al. 2010, A&A, 517, A80
  * Fischer et al. (2014) Fischer, T., Hempel, M., Sagert, I., et al. 2014, European Physical Journal A, 50, 46
  * Friedman & Pandharipande (1981) Friedman, B., & Pandharipande, V. R. 1981, Nucl. Phys. A, 361, 502
  * Frick & Müther (2003) Frick, T., & Müther, H. 2003, Phys. Rev. C, 68, 034310
  * Frick et al. (2005) Frick, T., Muther, H., Rios, A., Polls, A., & Ramos, A. Phys. Rev. C, 71, 014313
  * Glendenning (1985) Glendenning, N.K. 1985, ApJ., 293, 470
  * Glendenning (1996) Glendenning, N.K., 1996, Compact Stars: Nuclear Physics, Particle Physics, and General Relativity, Springer Verlag
  * Grangè et al. (1987) Grangé, P., Cugnon, J., & Lejeune, A. 1987, Nucl. Phys. A, 473, 365
  * Haensel & Pichon (1994) Haensel, P., & Pichon, B. 1994, A&A 283, 313
  * Haidenbauer et al. (2017) Haidenbauer, J., Meißner, U.-G., Kaiser, N., & Weise, W. 2017, Eur. Phys. J. A, 53, 121
  * Haidenbauer et al. (2016) Haidenbauer, J., Meißner, U.-G. & Petschauer, S. 2016, Nucl. Phys. A, 954, 273
  * Haidenbauer et al. (2013) Haidenbauer, J., Petschauer, S., Kaiser, N., Meißner, U.-G., Nogga, A., & Weise, W. 2013, Nucl. Phys. A, 915, 24
  * Hammer et al. (2013) Hammer, H. W., Nogga, A., & Schenk A. 2013, Rev. Mod. Phys., 85, 197
  * Hanauske & Bovard (2018) Hanauske, M., & Bovard, L. 2018, Journal of Astrophysics and Astronomy, 39, 45
  * Hempel & Schaffner-Bielich (2010) Hempel, M., & Schaffner-Bielich, J. 2010, Nucl. Phys. A, 837, 210
  * Holt et al. (2010) Holt, J. W., Kaiser, N., & Weise, W. 2010, Phys. Rev. C, 81, 024002
  * Hotokezaka et al. (2013) Hotokezaka, K., Kiuchi, K., Kyutoku, K., et al. 2013, Phys. Rev. D, 88, 044026
  * Hüdepohl et al. (2010) Hüdepohl, L., Müller, B., Janka, H.-T., et al. 2010, Phys. Rev. Lett., 104, 251101
  * Hüfner & Mahaux (1972) Hüfner, J., & Mahaux, C. 1972, Ann. Phys. (N.Y.), 73, 525
  * Hugenholtz & Van Hove (1958) Hugenholtz, N. M., & Van Hove, L. 1958, Physica (Amsterdam), 24, 363
  * Janka et al. (2016) Janka, H.-T., Melson, T., & Summa, A. 2016, Annual Review of Nuclear and Particle Science, 66, 341
  * Jeukenne et al. (1976) Jeukenne, J. P., Lejeunne, A., & Mahaux, C. 1976, Phys. Rep., 25, 83
  * Kalantar-Nayestanak et al. (2012) Kalantar-Nayestanaki, N., Epelbaum, E., Messchendorp, J. G., & Nogga, A. 2012, Rep. Prog. Phys., 75, 016301
  * Kastaun & Galeazzi (2015) Kastaun, W., & Galeazzi, F. 2015, Phys. Rev. D, 91, 064027
  * Kievsky et al. (2018) Kievsky, A., Viviani, M., Logoteta, D., Bombaci, I., & Girlanda, L. 2018 Phys. Rev. Lett., 121, 072701
  * Kiuchi et al. (2014) Kiuchi, K., Kyutoku, K., Sekiguchi, Y., Shibata, M. & Wada, T. 2014 Phys. Rev. D, 90, 041502
  * Lattimer & Swesty (1991) Lattimer, J., & Swesty D. F. 1991, Nucl. Phys. A, 535, 331
  * Lattimer & Prakash (2016) Lattimer, J., & Prakash, M. 2016, Physics Reports, 621, 127
  * Li et al. (2008) Li, Z. H., Lombardo, U., Schulze, H.-J., & Zuo, W. 2008, Phys. Rev. C, 77, 034316
  * Li et al. (2006) Li, Z. H., Lombardo, U., Schulze, H.-J., Zuo, W., Chen, L. W., & Ma, H. R. 2006, Phys. Rev. C, 74, 047304
  * Li & Schulze (2008) Li, Z. H. & Schulze, H.-J. 2008, Phys. Rev. C , 78, 028801
  * Liebendörfer et al. (2001) Liebendörfer, M., Mezzacappa, A., & Thielemann, F.-K. 2001, Phys. Rev. D, 63, 104003
  * Liebendörfer et al. (2009) Liebendörfer, M., Whitehouse, S. C., & Fischer, T. 2009, ApJ, 698, 1174
  * Logoteta et al. (2012) Logoteta, D., Bombaci, I., Providência, C., & Vidaña, I. 2012, Phys. Rev. D, 85, 023003
  * Logoteta et al. (2016b) Logoteta, D., Bombaci, I., & Kievsky, A. 2016b, Phys Rev. C, 94, 064001
  * Logoteta et al. (2019) Logoteta, D., Vidaña, I. & Bombaci, I. 2019, European Physical Journal A, 55, 207
  * Mezzacappa & Bruenn (1993) Mezzacappa, A., & Bruenn, S. W. 1993, ApJ, 405, 637
  * Miller et al. (2019) Miller, M. C. et al., 2019, Astrophys. J. Lett., 887, L24
  * Most et al. (2020) Most, E. R., Jens Papenfort, L., Dexheimer, V., et al. 2020, European Physical Journal A, 56, 59
  * Oertel et al. (2017) Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, Rev. Mod. Phys., 89, 015007
  * O’ Connor et al. (2018) E. O’ Connor et al., 2018, J. Phys. G: Nucl. Part. Phys., 45, 104001
  * Oppenheimer & Volkoff (1939) Oppenheimer, J., & Volkoff, G. 1939, Phys. Rev, 55, 374
  * Pan et al. (2019) Pan, K.-C., Mattes, C., O’Connor, E. P., et al. 2019, Journal of Physics G Nuclear Physics, 46, 014001
  * Perego et al. (2019) Perego, A., Bernuzzi, S. & Radice, D. Eur. Phys. J. A 2019, 55, 124
  * Piarulli et al. (2016) Piarulli M., Girlanda, L., Schiavilla, R., Kievsky, A., Lovato, A., Marcucci, L. E., Pieper, S. C., Viviani, M. & Wiringa, R. B. 2016, Phys. Rev. C, 94, 054007
  * Piarulli et al. (2020) Piarulli M., Bombaci, I., Logoteta, D., Lovato A., & Wiringa, R. B. 2020, Phys. Rev. C, 101, 045801
  * Pons et al. (1999) Pons, J. A., Reddy, S., Prakash, M., Lattimer, J. M., & Miralles, J. A. 1999, ApJ, 513, 780
  * Prakash et al. (1997) Prakash M., Bombaci I., Prakash M., Ellis P. J., Lattimer J. M., & Knorren R. 1997, Phys. Rep., 280, 1
  * Radice et al. (2018a) Radice, D., Perego, A., Hotokezaka, K., et al. 2018, ApJ, 869, 130
  * Radice et al. (2018b) Radice, D., Perego, A., Bernuzzi, S., et al. 2018, MNRAS, 481, 3670
  * Radice et al. (2020) Radice, D., Bernuzzi, S., & Perego, A. 2020, arXiv e-prints, arXiv:2002.03863
  * Rajaraman & Bethe (1967) Rajaraman, R., & Bethe, H. A. 1967, Rev. Mod. Phys., 39, 745
  * Rezzolla & Takami (2016) Rezzolla, L., & Takami, K. 2016, Phys. Rev. D, 93, 124051
  * Riley et al. (2019) Riley, T. E. et al., 2019, Astrophys. J. Lett., 887, L21
  * Rios et al. (2006) Rios A., Polls, A., Ramos, A., & Muther, H., 2006 Phys. Rev. C 74, 054317
  * Roberts et al. (2012) Roberts, L. F., Shen, G., Cirigliano, V., et al. 2012, Phys. Rev. Lett., 108, 061103
  * Roberts & Reddy (2017) Roberts, L. F., & Reddy, S. 2017, Handbook of Supernovae, 1605
  * Schneider et al. (2019) Schneider, A. S., Roberts, L. F., Ott, C. D., et al. 2019, Phys. Rev. C, 100, 055802
  * Shen et al. (1998) Shen, H., Toki, H., Oyamatsu, K., & Sumiyoshi, K. 1998, Nucl. Phys. A, 637, 435
  * Shibata & Hotokezaka (2019) Shibata, M., & Hotokezaka, K. 2019, Annual Review of Nuclear and Particle Science, 69, 41
  * Shibata et al. (2005) Shibata, M., Taniguchi, K., & Uryu, K. 2005, Phys. Rev. D, 71, 084021
  * Schneider et al. (2017) Schneider, A. S., Roberts, L. F., & Ott, C. D. 2017, Phys. Rev. C 96, 065802
  * Song et al. (1998) Song H. Q., Baldo M., Giansiracusa G., & Lombardo U. 1998 Phys. Rev. Lett., 81, 1584
  * Steiner et al. (2013c) Steiner, A. W., Hempel, M., & Fisher, T. 2013, Astrophys. J., 774, 17
  * Strobel et al. (1999) Strobel, K., Schaab, C., & Weigel, M. K. 1999, A&A, 350, 497
  * Sugahara & Toki (1994) Sugahara, Y., & Toki, H. 1994 Nucl. Phys. A, 579, 557
  * Timmes & Arnett (1999) Timmes, F. X., & Arnett, D. 1999, ApJS, 125, 277
  * Togashi et al. (2017) Togashi, H., Nakazato, K., Takehara, Y., Yamamuro, S., Suzuki, H., & Takano, M. 2017, Nucl. Phys. A, 961, 78
  * Tolman (1934) Tolman, R. C. 1934, Proc. Nat. Acad. Sci. (USA), 20, 3
  * Typel et al. (2010) Typel, S., Ropke, G., Klahn, T., Blaschke, D., & Wolter, H. H. 2010, Phys. Rev. C, 81, 015803
  * Vidaña & Bombaci (2002) Vidaña, I., & Bombaci, I. 2002, Phys. Rev. C, 66, 045801
  * Vidaña et al. (2011) Vidaña, I. Logoteta, D., Providência, C., Polls, A., & Bombaci, I. 2011, EPL (Europhysics Letters), 94, 11002
  * Wellenhofer et al. (2016) Wellenhofer, C., Holt, J. W., & Kaiser, N. 2016, Phys. Rev. C, 93, 055802
  * Wellenhofer et al. (2015) Wellenhofer, C., Holt, J. W., & Kaiser, N. 2015, Phys. Rev. C, 92, 015801
  * Wiringa et al. (1988) Wiringa, R.B., Fiks, V., & Fabrocini, A. 1998, Phys. Rev. C, 38, 1010
  * Woosley & Heger (2007) Woosley, S. E., & Heger, A. 2007, Phys. Rep, 442, 269
  * Zappa et al. (2018) Zappa, F., Bernuzzi, S., Radice, D., et al. 2018, Phys. Rev. Lett., 120, 111101
  * Zuo et al. (2014) Zuo, W., Bombaci, I., & Lombardo, U. 2014, Eur. Phys. J. A, 50, 12
  * Zuo et al. (1999) Zuo, W., Bombaci, I., & Lombardo, U. 1999, Phys. Rev. C., 60, 024605
  * Zuo et al. (2006) Zuo, W., Li, Z. H., Lombardo, U., Lu, G. C., & Schulze H.-J. 2006, Phys. Rev. C., 73, 035208
  * Zuo et al. (2004) Zuo, W., Li, Z. H., Li, A., & Lu, G. C. 2004, Phys. Rev. C., 69, 064001