© The Authors 2025

1. Introduction

Numerical simulations 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); 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); evolution of proto-neutron stars (PNSs) (Pons et al. 1999; Pascal et al. 2022) and formation of black holes (BHs) in failed CCSNe (Sumiyoshi et al. 2007; Fischer et al. 2009; O’Connor & Ott 2011; Hempel et al. 2012) require accurate microphysics input of particle composition, thermodynamic properties and reaction rates (Oertel et al. 2017). Composition and thermodynamic information is customarily provided in tabular form of the so-called general purpose equation of state (EOS) tables over wide domains of temperature (0 ≤ T ≤ 100 MeV), density (10⁻¹⁴ fm⁻³ ≤ n_B ≤ 1.5 fm⁻³) and electron fraction (0 ≤ Y_e = n_e/n_B ≤ 0.6) and with a mesh fine enough to allow interpolation and computation of additional thermodynamic state variables by differentiation (Typel 2022). Some EOS tables also contain information on specific microscopic quantities, for example, effective masses and interaction potentials as well as neutrino opacities.

The increased interest in the physics of neutron stars (NSs) and BNS mergers, largely motivated by the progress allowed by multi-messenger astronomy over the past decade, has resulted in a significant number of general purpose EOS tables being built and made available for numerical simulations. The COMPOSE online repository¹ hosts more than 120 of such tables and is still expanding. This collection reflects the wealth of assumptions regarding the particle degrees of freedom, theoretical approaches, and effective interactions. It is also illustrative of the uncertainties that affect states of matter that are impossible to produce and study in terrestrial laboratories.

Systematic studies 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) have demonstrated that the evolution of these phenomena is linked to the behavior of dense matter, with the nucleon effective mass (m_eff) playing the most important role. During the post-bounce evolution, m_eff impacts the properties of both PNSs and the neutrinosphere. Schneider et al. (2019), Yasin et al. (2020) showed that large values of m_eff result in higher (lower) central densities (temperatures) in the PNS core and lower PNS radii. Clear positive (negative) correlations have also been identified between m_eff and the neutrinosphere’s temperature and proton fraction (density and radius) as well as with neutrino energies and luminosities (Schneider et al. 2019). Schneider et al. (2020) proved that for a failed CCSN, the collapse into a BH occurs earlier for higher values of m_eff.

Studies of CCSN have also indicated some sensitivity to the composition of subsaturated (n_B ≤ n_sat ≈ 0.16 fm⁻³ ≈ 2 ⋅ 10¹⁴ g/cm³) nuclear matter (NM). An example in this sense is offered by Pascal et al. (2020), who showed that the single nucleus approximation and nuclear statistical equilibrium (NSE) versions of the same EOS lead to different neutronization degrees of stellar matter during the in-fall stage of a CCSN and different neutrino mean free paths. Pascal et al. (2020) explain these discrepancies in terms of differences between the electron capture rates of the unique representative nucleus in the single nucleus approximation and the electron capture rates of the ensemble of nuclei in the NSE distribution, respectively. The role of light clusters for CCSNe and neutrino driven winds was considered by Sumiyoshi & Roepke (2008) and Arcones et al. (2008), respectively. According to Sumiyoshi & Roepke (2008), the abundant production of light clusters affects the effectiveness of neutrino reheating of the shock wave as well as the structure and composition of PNSs. According to Arcones et al. (2008), light nuclei have a small impact on the average energy of the emitted electron neutrinos but are significant for the average energy of antineutrinos. At temperatures around 0.5 MeV, the properties of subsaturated NM also govern the processes of formation and crystallization of NS crusts, thus determining their final composition and their thermal and electrical conductivities (Carreau et al. 2020a; Fantina et al. 2020; Potekhin & Chabrier 2021).

The main aim of this paper is to propose a new set of general purpose EOSs that, upon input in numerical astrophysical simulations, will contribute to a better understanding of the role that the nuclear EOS plays. All of these EOSs employ Brussels extended Skyrme effective interactions (Chamel et al. 2009). They are built within a Bayesian inference of the EOS of dense matter (Beznogov & Raduta 2024) and are selected such that they feature distinct thermal behaviors in the suprasaturation regime. The subsaturation density domain, where NM is inhomogeneous, is treated within the NSE approximation. Homogeneous NM is treated within the non-relativistic mean-field model (Negele & Vautherin 1972; Vautherin 1996).

The high convergence performance of our new NSE code made it possible to approach the nuclear saturation density for a wide range of electron fractions and temperatures. In the first place, this allowed us to study in detail the mechanisms of transitioning to homogeneous matter. We show that, depending on T and Y_p, the transition occurs either due to the gas of nuclei or the gas of self-interacting unbound nucleons filling the entire volume. Special attention is also devoted to the stability of stellar matter as a mixture of nuclear clusters, unbound nucleons, and a gas of electrons. We demonstrate that the density-driven instabilities featured by homogeneous NM do not manifest in clusterized NM with electrons.

The paper is organized as follows. In Sect. 2, we present the NSE formalism along with how it is implemented in our model. The properties of a number of newly proposed Brussels extended Skyrme effective interactions are also described. In Sect. 3 we discuss the various mechanisms by which the transition to homogeneous matter occurs in our model. The thermodynamic stability of a cold mixture of nuclei, nucleons, and electrons is addressed in Sect. 4, and special attention is given to the stability of beta-equilibrated matter. In Sect. 5, the results of our NSE model are investigated over wide domains of temperature, electron fraction, and density. The composition of matter is analyzed in terms of three generic components: heavy and light nuclei, and unbound nucleons. Then, the density dependence of a few thermodynamic state variables is analyzed for fixed values of T (Y_p) and variable values of Y_p (T). We used the NSE results at low temperatures to infer the composition of NS crusts. The conclusions are drawn in Sect. 6.

2. Extended nuclear statistical equilibrium

For densities lower than the normal nuclear saturation density n_sat and temperatures lower than a dozen mega electron volt, stellar matter consists of a mixture of nuclei, nucleons, leptons, and photons. Charged leptons are essential for compensating for the positive charge of the protons. We assumed that charge neutrality holds locally: $n_p={\overline{n}}_e+{\overline{n}}_{\mu }$, where ${\overline{n}}_{e/\mu }=n_{e^{-}/{\mu }^{-}}-n_{e^{+}/{\mu }^{+}}$ stands for the net number of electrons and muons and n_p, n_e⁻, n_μ⁻, n_e⁺, and n_μ⁺ respectively represent the number densities of protons, electrons, muons, positrons, and anti-muons. The only significant interaction between electrons and protons is the electrostatic one, for which the Wigner-Seitz approximation (Baym et al. 1971) is customarily employed. The electrons are assumed to be uniformly distributed throughout the volume, and the electron gas is treated as an ideal Fermi gas. Photons are thermalized and are described as an ideal Bose gas. Neutrinos are not necessarily in equilibrium with the rest of matter and are disregarded here. For simplicity, we also disregard muons. Accounting for them is straightforward.

For very low temperatures, the nuclear composition of stellar matter is decided upon a reaction network. As the temperature increases, more nuclear reactions become possible and their cross sections increase. At T ≈ 4 − 5 GK, all nuclear reactions proceed at the same rate as their inverses, which makes that the species’ abundances are not decided any longer by the reaction rates (Wiescher & Rauscher 2018). This is the so-called NSE regime, where species abundances are determined by nuclear binding energies and internal partition functions.

In the evolution of the CCSNe, PNSs, and BNS mergers, densities higher than n_sat and temperatures up to several tens of MeV are reached. The obvious need for a unified treatment of the sub and suprasaturation density domains urged the community to further develop the NSE approach to the point of approaching the transition to homogeneous matter and more appropriately dealing with configurations dominated by unbound nucleons. As such, the so-called extended NSE (eNSE) approach arose. In contrast with the standard NSE, where nuclei are treated as an ideal gas, the extended version accounts for interactions among unbound nucleons on the one hand and among unbound nucleons and nuclei on the other hand. The ensemble of self-interacting unbound nucleons is treated within the mean-field theory of NM, which is also used to describe the matter in the suprasaturation regime. Interactions among nuclei and with unbound nucleons are phenomenologically implemented via the excluded volume. As we show in the following, this schematic treatment turns out to be sufficient for a reasonable description of matter with densities in the proximity of the transition to homogeneous matter.

Several eNSE models have been proposed in the last 15 years. The models by Hempel & Schaffner-Bielich (2010), Pais & Typel (2017), Typel (2018) and Furusawa et al. (2011, 2013, 2017a,b) treat the gas of unbound self-interacting nucleons within the covariant density functional theory of NM while the models by Raduta & Gulminelli (2010, 2019), Gulminelli & Raduta (2015) employ the non-relativistic mean-field model of NM. While in all of these models nuclei are assumed to form an ideal gas, their description varies significantly from one implementation to another. For the binding energies, experimental values, evaluations by dedicated models, or liquid drop parametrizations have been employed. Temperature-induced excitations are accounted for via various level density expressions; the assumptions made in regard to the upper limit of the excitation energy range from the minimum between neutron and proton separation energies up to infinity. The generalized relativistic mean-field model by Pais & Typel (2017), Typel (2018) implements in-medium modifications of the nuclear energy functional via parametrized mass shifts derived from a quantum statistical model (Typel et al. 2010). These shifts mimic the Pauli blocking of clusters’ states by the nucleons in the medium and guarantee that nuclei dissolve in high-density media. They are so efficient in accounting for interactions among nuclei and nucleons that they render the excluded volume approximation used by other eNSE models useless. Temperature and in-medium effects on nuclear energy functional are phenomenologically accounted for also by Furusawa et al. (2011, 2013, 2017a,b). Of course, many nucleon-nucleon effective interactions have been used, which is mainly reflected in different behaviors of matter with densities higher than the saturation density and temperatures higher than 10 − 20 MeV. For a comparative analysis of various NSE models, see Buyukcizmeci (2013), Raduta et al. (2021). In the following, we describe the model proposed in this work.

2.1. Formalism

Here, the canonical ensemble is used. The controlled quantities are the total number of nucleons (A_tot), the total number of protons (Z_tot), the temperature (T) and the volume (V). Net charge neutrality is assumed to hold locally, which means that at any point the proton density n_p = Z_tot/V is identically equal to the net density of electrons ${\overline{n}}_{\mathrm{el}}={\overline{N}}_{\mathrm{el}}/V$. For convenience, in the following we refer to the electron-positron gas as electron gas and omit the bar over n_el.

The ensemble of nuclei is assumed to form an ideal gas. Unbound nucleons are assumed to have a homogeneous distribution. Their self-interactions are accounted for within the mean-field theory of NM, for which the non-relativistic model (Negele & Vautherin 1972; Vautherin 1996) is employed here. Interactions among nuclei and with unbound nucleons occur via the excluded volume. As we show later, this is essential for the transition to homogeneous matter at densities n_B of the order of n_sat/2. The nuclear component interacts with the electron gas via electrostatic interactions only.

These assumptions allow one to decompose the total free energy of the system into the sum of free energies:

$$\begin{array}{c}\hfill F={\displaystyle \sum _{A,Z}}(F_{A,Z}^{\mathrm{ideal}}+F_{A,Z}^{\mathrm{int}})+F_{\mathrm{Coul}}+F_{\mathrm{MF}}+F_{\mathrm{el}}.\end{array}$$(1)

The ideal-part of the free energy of a generic nucleus with the mass number A and charge number Z is

$$\begin{array}{cc}\hfill F_{A,Z}^{\mathrm{ideal}}(N_{A,Z},T)& =N_{A,Z}{M}_{A,Z}{c}^2\hfill \\ \hfill & -N_{A,Z}T[ln\left(\frac{V_{\mathrm{cl}}^{\mathrm{free}}}{N_{A,Z}}{\left(\frac{M_{A,Z}T}{2\pi {ħ}^2}\right)}^{3/2}\right)+1],\hfill \end{array}$$(2)

where M_A, Zc² is the rest-mass contribution and M_A, Z and N_A, Z denote the mass and number of the (A, Z) nuclei; $V_{\mathrm{cl}}^{\mathrm{free}}=V-V_{\mathrm{gas}}-V_{\mathrm{cl}}$ represents the volume available for the translational movement of nuclei after the volumes occupied by nucleons and all nuclear species is subtracted; V_gas = (N_n + N_p)v₀ represents the intrinsic volume of the unbound N_n neutrons and N_p protons, with each occupying the volume v₀ = 1/n_sat; V_cl = ∑_A, ZN_A, Zv_A, Z stands for the intrinsic volume of the all nuclei. To account for the isospin asymmetry (δ = 1 − 2Z/A) dependence of the volume, we take v_A, Z = A/n_sat(δ).

Population of excited states, made possible by finite temperatures, asks for consideration of the contribution from internal degrees of freedom in addition to the one due to translation. This contribution for the species (A, Z) is

$$\begin{array}{c}\hfill F_{A,Z}^{\mathrm{int}}=-Tln{z}_{A,Z}^{\mathrm{int}}=-T{N}_{A,Z}ln{\stackrel{\sim }{z}}_{A,Z}^{\mathrm{int}},\end{array}$$(3)

and it depends on the species abundance and its excitation energy spectrum. Here,

$$\begin{array}{cc}\hfill {\stackrel{\sim }{z}}_{A,Z}^{\mathrm{int}}& ={\displaystyle \sum _j}{g}_{j}exp(-{ϵ}_j/T)\hfill \\ \hfill & =g_{\mathrm{G}.\mathrm{S}.}+{\int }_0^{B(A,Z)}dϵ\rho (A,Z,ϵ)exp(-ϵ/T),\hfill \end{array}$$(4)

represents the internal partition function of the (A, Z) nucleus; g_j and ϵ_j stand for the spin degeneracy and energy of the j excited state; g_G.S. is the spin degeneracy of the ground state; ρ(A, Z, ϵ) is the level density; B(A, Z) = (A − Z)m_nc² + Zm_pc² − M_A, Zc² is the binding energy of the nucleus.

The Coulomb contribution to the free energy is the sum of the Coulomb energies of each nucleus and its surrounding gas of electrons,

$$\begin{array}{c}\hfill F_{\mathrm{Coul}}={\displaystyle \sum _{A,Z}}{N}_{A,Z}{E}_{A,Z}^{\mathrm{Coul}},\end{array}$$(5)

where (Lattimer et al. 1985),

$$\begin{array}{c}\hfill E_{A,Z}^{\mathrm{Coul}}=-\frac{3}{5}\frac{Z^2{e}^2}{R_{A,Z}}(\frac{3}{2}x-\frac{1}{2}{x}^3),\end{array}$$(6)

with x = [A/Z⋅n_el/n_sat(δ)]^1/3, R_A, Z = [3v_A, Z/4π]^1/3, n_el = [N_p+∑_A, ZZN_A, Z]/V. $E_{A,Z}^{\mathrm{Coul}}$ above has two components of opposite signs. The attractive component comes from the interaction between protons and electrons. The negative one is due to the mutual repulsion between electrons. Both terms depend on the electron density and spatial extension of the nucleus. We note that the electrostatic interaction among the protons in a nucleus does not appear in Eq. (6), because it is accounted for in the nuclear mass entering Eq. (2).

The total free energy of the unbound nucleons can be expressed as

$$\begin{array}{cc}\hfill F_{\mathrm{MF}}(N_n,N_p,T)& =N_n{m}_n{c}^2+N_p{m}_p{c}^2\hfill \\ \hfill & +V_{\mathrm{gas}}^{\mathrm{free}}{f}_{\mathrm{MF}}^0(\frac{N_n}{V_{\mathrm{gas}}^{\mathrm{free}}},\frac{N_p}{V_{\mathrm{gas}}^{\mathrm{free}}},T),\hfill \end{array}$$(7)

where $V_{\mathrm{gas}}^{\mathrm{free}}=(V-V_{\mathrm{cl}})$ is the volume available for these nucleons. As in Eq. (2), the rest-mass contribution is explicitly written. We note that, at variance with the nuclei, unbound nucleons do not geometrically exclude each other. This treatment is consistent with the mean-field approach used to model their behavior, where nucleons are described by plane waves. f_MF⁰ represents the mean-field value of the free energy density of a gas with the neutron and proton densities ${\stackrel{\sim }{n}}_n=N_n/V_{\mathrm{gas}}^{\mathrm{free}}$ and ${\stackrel{\sim }{n}}_p=N_p/V_{\mathrm{gas}}^{\mathrm{free}}$ and that does not interact with other species (thus, the superscript “0”). As a matter of fact, throughout this paper we mark by a “0” superscript the values that state variables take when the subsystem under consideration is treated as standalone. Examples of interactions that are already considered are the interactions among unbound nucleons that are accounted for in the mean-field model, the interactions among unbound nucleons and nuclei that are modeled via the excluded volume, the Coulomb interactions between the nuclei and the electrons, and among electrons that are explicitly dealt with in the Coulomb terms.

The free energy of the electron gas with density n_el is

$$\begin{array}{c}\hfill F_{\mathrm{el}}(N_{\mathrm{el}},T)=V{f}_{\mathrm{el}}^0(n_{\mathrm{el}},T).\end{array}$$(8)

In the following, we derive the expressions for some basic thermodynamic state variables upon differentiation of the free energy in Eq. (1). We start with the chemical potential of an (A, Z) nucleus defined as

$$\begin{array}{cc}\hfill {\mu }_{A,Z}& ={\frac{\partial F}{\partial N_{A,Z}}|}_{T,V,N_{\mathrm{el}},N_n,N_p,\{N_{A\prime ,Z\prime }\}}\hfill \\ \hfill & ={\mu }_{A,Z}^{\mathrm{ideal}}+{\mu }_{A,Z}^{\mathrm{int}}+{\mu }_{A,Z}^{\mathrm{Coul}}+{\mu }_{A,Z}^{\mathrm{MF}}+{\mu }_{A,Z}^{\mathrm{el}},\hfill \end{array}$$(9)

and explicitate the contribution coming from each of the terms in Eq. (1).

We note that in addition to an explicit dependence on N_A, Z, $F_{A,Z}^{\mathrm{ideal}}$ depends on N_A, Z also via $V_{\mathrm{cl}}^{\mathrm{free}}$. The first of these dependencies gives rise to

$$\begin{array}{c}\hfill {\mu }_{A,Z}^{\mathrm{ideal};1}=M_{A,Z}{c}^2-Tln\left[\frac{V_{\mathrm{cl}}^{\mathrm{free}}}{N_{A,Z}}{\left(\frac{M_{A,Z}T}{2\pi {ħ}^2}\right)}^{3/2}\right]+T\frac{v_{A,Z}{N}_{A,Z}}{V_{\mathrm{cl}}^{\mathrm{free}}},\end{array}$$(10)

while the second leads to

$$\begin{array}{c}\hfill {\mu }_{A,Z}^{\mathrm{ideal};2}=T\frac{v_{A,Z}}{V_{\mathrm{cl}}^{\mathrm{free}}}{\displaystyle \sum _{A\prime ,Z\prime \ne A,Z}}{N}_{A\prime ,Z\prime }.\end{array}$$(11)

From these two contributions, one obtains

$$\begin{array}{cc}\hfill {\mu }_{A,Z}^{\mathrm{ideal}}=M_{A,Z}{c}^2& -Tln\left[\frac{V_{\mathrm{cl}}^{\mathrm{free}}}{N_{A,Z}}{\left(\frac{M_{A,Z}T}{2\pi {ħ}^2}\right)}^{3/2}\right]+T\frac{v_{A,Z}}{V_{\mathrm{cl}}^{\mathrm{free}}}{\displaystyle \sum _{A\prime ,Z\prime }}{N}_{A\prime ,Z\prime },\hfill \end{array}$$(12)

where the summation in the last term applies to all species, without exception.

The contributions from internal excitation and Coulomb interaction are trivial and write as

$$\begin{array}{cc}& {\mu }_{A,Z}^{\mathrm{int}}=-Tln{\stackrel{\sim }{z}}_{A,Z}^{\mathrm{int}},\hfill \end{array}$$(13)

$$\begin{array}{cc}& {\mu }_{A,Z}^{\mathrm{Coul}}=E_{A,Z}^{\mathrm{Coul}},\hfill \end{array}$$(14)

respectively.

Now, N_A, Z enters the mean-field contribution to the total free energy via the volume of the (A, Z) nuclei that is forbidden to the unbound nucleons. As such,

$$\begin{array}{cc}\hfill {\mu }_{A,Z}^{\mathrm{MF}}& =-v_{A,Z}{f}_{\mathrm{MF}}^0({\stackrel{\sim }{n}}_n,{\stackrel{\sim }{n}}_p)+(V-V_{\mathrm{cl}})[\frac{\partial f_{\mathrm{MF}}^0}{\partial {\stackrel{\sim }{n}}_n}\frac{\partial {\stackrel{\sim }{n}}_n}{\partial N_{A,Z}}+\frac{\partial f_{\mathrm{MF}}^0}{\partial {\stackrel{\sim }{n}}_p}\frac{\partial {\stackrel{\sim }{n}}_p}{\partial N_{A,Z}}]\hfill \\ \hfill & =-v_{A,Z}{f}_{\mathrm{MF}}^0({\stackrel{\sim }{n}}_n,{\stackrel{\sim }{n}}_p)+{\mu }_{\mathrm{MF};n}^0\frac{v_{A,Z}{N}_n}{(V-V_{\mathrm{cl}})}+{\mu }_{\mathrm{MF};p}^0\frac{v_{A,Z}{N}_p}{(V-V_{\mathrm{cl}})}\hfill \\ \hfill & =v_{A,Z}{p}_{\mathrm{MF}}^0({\stackrel{\sim }{n}}_n,{\stackrel{\sim }{n}}_p),\hfill \end{array}$$(15)

where μ_MF; i⁰ with i = n, p represent the mean-field chemical potentials of neutrons and protons with densities ${\stackrel{\sim }{n}}_i$.

Considering now that the electron free energy depends only on control quantities, that is, the total number of protons, volume, and temperature, one gets ${\mu }_{A,Z}^{\mathrm{el}}=0$. Putting together the expressions in Eqs. (12), (13), (14) and (15), it turns out that the chemical potential of an (A, Z) nucleus is given by

$$\begin{array}{cc}\hfill {\mu }_{A,Z}& =M_{A,Z}{c}^2-Tln\left[\frac{V_{\mathrm{cl}}^{\mathrm{free}}}{N_{A,Z}}{\left(\frac{M_{A,Z}T}{2\pi {ħ}^2}\right)}^{3/2}\right]+T\frac{v_{A,Z}}{V_{\mathrm{cl}}^{\mathrm{free}}}{\displaystyle \sum _{A\prime ,Z\prime }}{N}_{A\prime ,Z\prime }\hfill \\ \hfill & -Tln{\stackrel{\sim }{z}}_{A,Z}^{\mathrm{int}}+E_{A,Z}^{\mathrm{Coul}}+v_{A,Z}{p}_{\mathrm{MF}}^0({\stackrel{\sim }{n}}_n,{\stackrel{\sim }{n}}_p).\hfill \end{array}$$(16)

The chemical potential of unbound neutrons and protons can be obtained from the thermodynamic definition,

$$\begin{array}{c}\hfill {\mu }_{n/p}={\frac{\partial F}{\partial N_{n/p}}|}_{T,V,N_{\mathrm{el}},N_{p/n},\{N_{A,Z}\}},\end{array}$$(17)

and it leads to

$$\begin{array}{c}\hfill {\mu }_{n/p}=m_{n/p}{c}^2+\frac{T{v}_0}{V_{\mathrm{cl}}^{\mathrm{free}}}{\displaystyle \sum _{A,Z}}{N}_{A,Z}+{\mu }_{\mathrm{MF};n/p}^0.\end{array}$$(18)

The above equation shows that due to the interaction with the gas of nuclei, the chemical potentials of the unbound nucleons are shifted with respect to the values they would have if other particles did not exist in the considered volume. On the same line, the chemical potential of electrons,

$$\begin{array}{c}\hfill {\mu }_{\mathrm{el}}={\frac{\partial F}{\partial N_{\mathrm{el}}}|}_{T,V,N_n,N_p,\{N_{A,Z}\}},\end{array}$$(19)

gets modified by the Coulomb interaction with the nuclei and among themselves, which brings about the expression

$$\begin{array}{c}\hfill {\mu }_{\mathrm{el}}={\mu }_{\mathrm{Coul}}+{\mu }_{\mathrm{el}}^0,\end{array}$$(20)

where

$$\begin{array}{c}\hfill {\mu }_{\mathrm{Coul}}=-\frac{3}{10}{\displaystyle \sum _{A,Z}}{N}_{A,Z}\frac{Z^2{e}^2}{R_{A,Z}}(1-x^2)\frac{x}{N_{\mathrm{el}}}·\end{array}$$(21)

Here, μ_el⁰ represents the chemical potential the electrons would have if they did not suffer any interaction.

Next, we turn to the pressure defined as

$$\begin{array}{c}\hfill P=-{\frac{\partial F}{\partial V}|}_{T,N_{\mathrm{el}},N_n,N_p,\{N_{A,Z}\}},\end{array}$$(22)

which can be expressed as a sum of contributions that stem from the various terms in the free energy (see Eq. (1)):

$$\begin{array}{c}\hfill P={\displaystyle \sum _{A,Z}}{P}_{A,Z}^{\mathrm{ideal}}+{\displaystyle \sum _{A,Z}}{P}_{A,Z}^{\mathrm{int}}+P_{\mathrm{Coul}}+P_{\mathrm{MF}}+P_{\mathrm{el}},\end{array}$$(23)

with

$$\begin{array}{cc}& P_{A,Z}^{\mathrm{ideal}}=T{N}_{A,Z}\frac{1}{V_{\mathrm{cl}}^{\mathrm{free}}},\hfill \end{array}$$(24)

$$\begin{array}{cc}& P_{A,Z}^{\mathrm{int}}=0,\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill P_{\mathrm{Coul}}& =\frac{3}{5}{\displaystyle \sum _{A,Z}}{N}_{A,Z}\frac{Z^2{e}^2}{R_{A,Z}}\frac{3}{2}(1-x^2){\left(\frac{A}{Z{n}_{\mathrm{sat}}(\delta )}\right)}^{1/3}\frac{1}{3}{n}_{\mathrm{el}}^{-2/3}\frac{\partial n_{\mathrm{el}}}{\partial V}\hfill \\ \hfill & =-\frac{3}{10}{\displaystyle \sum _{A,Z}}{N}_{A,Z}\frac{Z^2{e}^2}{R_{A,Z}}(1-x^2)\frac{x}{V}\hfill \\ \hfill & =n_{\mathrm{el}}{\mu }_{\mathrm{Coul}},\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill P_{\mathrm{MF}}& =-f_{\mathrm{MF}}^0({\stackrel{\sim }{n}}_n,{\stackrel{\sim }{n}}_p)+(V-V_{\mathrm{cl}})[{\mu }_{\mathrm{MF};n}^0\frac{{\stackrel{\sim }{n}}_n}{V-V_{\mathrm{cl}}}+{\mu }_{\mathrm{MF};p}^0\frac{{\stackrel{\sim }{n}}_p}{V-V_{\mathrm{cl}}}]\hfill \\ \hfill & =p_{\mathrm{MF}}^0({\stackrel{\sim }{n}}_n,{\stackrel{\sim }{n}}_p),\hfill \end{array}$$(27)

and

$$\begin{array}{c}\hfill P_{\mathrm{el}}=P_{\mathrm{el}}^0.\end{array}$$(28)

Similarly, the entropy of the mixture of unbound nucleons, nuclei, and electrons,

$$\begin{array}{c}\hfill S=-{\frac{\partial F}{\partial T}|}_{V,N_{\mathrm{el}},N_n,N_p,\{N_{A,Z}\}},\end{array}$$(29)

can be decomposed into

$$\begin{array}{c}\hfill S={\displaystyle \sum _{A,Z}}{S}_{A,Z}^{\mathrm{ideal}}+{\displaystyle \sum _{A,Z}}{S}_{A,Z}^{\mathrm{int}}+S_{\mathrm{Coul}}+S_{\mathrm{MF}}+S_{\mathrm{el}},\end{array}$$(30)

where

$$\begin{array}{cc}& S_{A,Z}^{\mathrm{ideal}}=N_{A,Z}[ln\left(\frac{V_{\mathrm{cl}}^{\mathrm{free}}}{N_{A,Z}}{\left(\frac{M_{A,Z}T}{2\pi {ħ}^2}\right)}^{3/2}\right)+\frac{5}{2}],\hfill \end{array}$$(31)

$$\begin{array}{cc}& S_{A,Z}^{\mathrm{int}}=N_{A,Z}[ln{\stackrel{\sim }{z}}_{A,Z}^{\mathrm{int}}+T\frac{\partial ln{\stackrel{\sim }{z}}_{A,Z}^{\mathrm{int}}}{\partial T}],\hfill \end{array}$$(32)

$$\begin{array}{cc}& S_{\mathrm{Coul}}=0,\hfill \end{array}$$(33)

$$\begin{array}{cc}& S_{\mathrm{MF}}=S_{\mathrm{MF}}^0({\stackrel{\sim }{n}}_n,{\stackrel{\sim }{n}}_p),\hfill \end{array}$$(34)

$$\begin{array}{cc}& S_{\mathrm{el}}=S_{\mathrm{en}}^0(N_{\mathrm{en}}^0).\hfill \end{array}$$(35)

The total energy can be computed from the thermodynamic identity,

$$\begin{array}{c}\hfill E=F+TS,\end{array}$$(36)

and it provides,

$$\begin{array}{cc}\hfill E& ={\displaystyle \sum _{A,Z}}{N}_{A,Z}[M_{A,Z}{c}^2+\frac{3}{2}T+⟨{ϵ}_{A,Z}⟩+E_{A,Z}^{\mathrm{Coul}}]\hfill \\ \hfill & +(V-V_{\mathrm{cl}})[m_n{c}^2{\stackrel{\sim }{n}}_n+m_p{c}^2{\stackrel{\sim }{n}}_p+e_{\mathrm{MF}}^0({\stackrel{\sim }{n}}_n,{\stackrel{\sim }{n}}_p)]\hfill \\ \hfill & +V{e}_{\mathrm{el}}^0(n_{\mathrm{el}}),\hfill \end{array}$$(37)

where

$$\begin{array}{cc}\hfill ⟨{ϵ}_{A,Z}⟩& =\frac{{\displaystyle {\int }_0^{B(A,Z)}}dϵϵ\rho (ϵ)exp(-ϵ/T)}{{\displaystyle {\int }_0^{B(A,Z)}}dϵ\rho (ϵ)exp(-ϵ/T)}\hfill \\ \hfill & =T^2\frac{\partial ln{\stackrel{\sim }{z}}_{A,Z}^{\mathrm{int}}}{\partial T}\hfill \end{array}$$(38)

represents the average excitation energy of the nucleus (A, Z).

It is interesting to point out that the thermodynamic identity F − PV + ∑_{i = n, p, {A, Z},el}μ_iN_i = 0 holds for the mixture as a whole but not for individual components, that is, F_i − P_iV − μ_iN_i ≠ 0 where i corresponds to each of the following: the gas of unbound self-interacting nucleons, the ensemble of nuclei and the electron gas. This situation is due to the strong coupling between the three subsystems and the fact that it is impossible to attribute the Coulomb terms to any particular one.

Strong and electromagnetic equilibrium between the homogeneous and clusterized nuclear components requires that

$$\begin{array}{c}\hfill {\mu }_{A,Z}=(A-Z){\mu }_n+Z{\mu }_p.\end{array}$$(39)

Substitution of expressions in Eqs. (16) and (18) leads to

$$\begin{array}{cc}& -Tln\left[\frac{V_{\mathrm{cl}}^{\mathrm{free}}}{N_{A,Z}}{\left(\frac{M_{A,Z}T}{2\pi {ħ}^2}\right)}^{3/2}{\stackrel{\sim }{z}}_{A,Z}^{\mathrm{int}}\right]=(A-Z){\mu }_{\mathrm{MF},n}^0+Z{\mu }_{\mathrm{MF},p}^0\hfill \\ \hfill & +(A-Z)m_n{c}^2+Z{m}_p{c}^2-E_{A,Z}^{\mathrm{Coul}}-A{v}_0{p}_{\mathrm{MF}}^0-M_{A,Z}{c}^2,\hfill \end{array}$$(40)

from which it is straightforward to extract the number of (A, Z) nuclei in the volume, V,

$$\begin{array}{cc}& N_{A,Z}=V_{\mathrm{cl}}^{\mathrm{free}}{z}_{A,Z}^{\mathrm{int}}{\left(\frac{M_{A,Z}T}{2\pi {ħ}^2}\right)}^{3/2}\hfill \\ \hfill & \times exp\left(\frac{(A-Z){\mu }_{\mathrm{MF},n}^0+Z{\mu }_{\mathrm{MF},p}^0-E_{A,Z}^{\mathrm{Coul}}-A{v}_0{p}_{\mathrm{MF}}^0+B(A,Z)}{T}\right).\hfill \end{array}$$(41)

Equation (41) reveals that an essential piece of information for determining the composition of the inhomogeneous phase is represented by the mean-field values of the chemical potentials of the unbound self-interacting nucleons. For fixed temperature T, density n_B and proton fraction Y_p, these are determined by solving the equations for mass and charge number conservation,

$$\begin{array}{cc}& n_{B}V=A_{\mathrm{tot}}=N_n+N_p+{\displaystyle \sum _{A,Z}}{N}_{A,Z}A,\hfill \end{array}$$(42)

$$\begin{array}{cc}& n_B{Y}_{p}V=Z_{\mathrm{tot}}=N_p+{\displaystyle \sum _{A,Z}}{N}_{A,Z}Z.\hfill \end{array}$$(43)

Despite its simplicity, the numerical resolution of the above system is highly non-trivial. The reason is that, via $V_{\mathrm{cl}}^{\mathrm{free}}$, N_A, Z depends on itself and on the intrinsic volume of unbound nucleons, which, in turn, depends on the unknown variables μ_MF, n/p⁰. We found it convenient to introduce a third unknown variable, V_cl, which makes the computation of $V_{\mathrm{cl}}^{\mathrm{free}}$ in Eq. (41) straightforward, and constrain its value by a third equation,

$$\begin{array}{c}\hfill V{\prime }_{\mathrm{cl}}={\displaystyle \sum _{A,Z}}{N}_{A,Z}{v}_{A,Z}.\end{array}$$(44)

Once μ_MF, n/p⁰ are known, both the composition of matter and its energetics can be determined using the equations in this section. However, a word of caution is in order here. The definition of the chemical potential of an (A, Z) nucleus in Eq. (9), which is a key ingredient of our derivation, assumes that variations of the total free energy with respect to the number of particles belonging to one species can be taken without altering the rest of the cluster gas composition. This assumption is consistent with the ideal gas working frame, where no interaction exists among the nucleus under consideration and other nuclei that form the gas of nuclear clusters. Similarly, Eqs. (22) and (29) assume that derivatives of the total free energy with respect to volume and temperature can be taken without modifying the composition of the mixture. Yet, Eq. (41) shows that {N_A, Z} are strongly coupled and the coupling becomes more important at high densities; it also shows that {N_A, Z} has an explicit dependence on temperature. Implicit dependence also exists through all ingredients other than M_A, Z and B(A, Z). We note that all of these dependencies arise exclusively due to the excluded volume. This corroborates previous microscopic evidence (Typel 2016) that, in spite of its apparent simplicity, the excluded volume approximation is very efficient in effectively accounting for interactions.

All of these equations are similar to those of Hempel & Schaffner-Bielich (2010). Nevertheless, the derivation is different.

2.2. Nuclei

In this work, we allow all nuclei with A ≥ 2, for which experimental data and/or theoretical evaluations of the mass are available, to be populated in the inhomogeneous phase. For experimental values of nuclear masses, we use the AME2020 (Wang et al. 2021) table, which is the most recent available to date. For theoretical estimations, we use the 10 parameter mass table by Duflo & Zuker (1995)². The heaviest nuclei present in (Wang et al. 2021) and (Duflo & Zuker 1995) have A = 295 and 297, respectively. The advantage of supplementing the pool of nuclei in the experimental mass table with nuclei for which mass estimates have been performed with dedicated theoretical models consists in accounting for nuclei far beyond the drip lines, which are expected to exist in astrophysical environments like those produced at different stages in the evolution of supernovae or binary mergers. With peculiar reaction rates, exotic nuclei impact the dynamical evolution as well as the nucleosynthesis of heavy nuclei at astrophysical sites, where they serve as initial conditions for nuclear reaction network calculations.

The isospin dependence of the nuclear volumes is implemented via the isospin dependence of the saturation density of NM, v_A, Z = A/n_sat(δ), as already explained in the previous section. For nuclei with an isospin asymmetry larger than the maximum value (δ_max) for which NM manifests saturation, v_A, Z = A/n_sat(δ_max).

The ground states of even (odd) mass nuclei are considered to have a spin degeneracy factor g_G.S. = 1 (2). For nuclei with A ≥ 16, the temperature-triggered population of excited states is accounted for via the energy-dependent back-shifted Fermi gas level density; in particular, we adopt the parametrization proposed by Bucurescu & von Egidy (2005), von Egidy & Bucurescu (2005). The upper limit of the excitation energy spectrum in the expression for the internal partition function of an (A, Z)-species is taken here to be B(A, Z), see Eq. (4).

The presented modeling of nuclear clusters in a hot and dense medium suffers from two obvious drawbacks. First, in-medium modifications of the nuclear surface energy, which increase with the density of the surrounding medium and are essential for the dissolution of nuclei in dense nucleonic environments, are treated phenomenologically. A more accurate implementation of such effects requires dedicated microscopic studies, which are not available yet. Second, by disregarding temperature effects on the binding energy, our approach allows nuclei to survive beyond the limiting temperature of the Coulomb instability (Jaqaman et al. 1984; Song & Su 1991). This feature is common for many eNSE implementations (Hempel & Schaffner-Bielich 2010; Furusawa et al. 2011, 2013, 2017b; Raduta & Gulminelli 2019) and explains why a significant fraction of nuclei with various masses exists at T > 15 MeV. For a comparison, see Fig. 17 in (Raduta et al. 2021).

To correct for the spurious persistence of nuclei at temperatures higher than 8–12 MeV, in our model inhomogeneous matter with subsaturation densities and T ≥ 2T_C/3, where T_C stands for the critical temperature of the liquid-gas phase transition of subsaturated homogeneous matter, is superseded by homogeneous matter. The value we choose for the limiting temperature of the Coulomb instability is based on microscopic calculations of Jaqaman et al. (1984), Song & Su (1991) and, in contrast with the results of Jaqaman et al. (1984), Song & Su (1991), disregards extra finite size and isospin dependencies.

2.3. Nucleons and effective interactions

Unbound nucleons are treated within the non-relativistic mean-field model (Negele & Vautherin 1972; Vautherin 1996) with zero-range interactions, which successfully describes the properties of ground-state and excited nuclei and NM at both zero and finite temperatures.

In the absence of spin polarization, the energy density of homogeneous matter with no Coulomb interaction can be expressed as the sum of four terms,

$$\begin{array}{c}\hfill \mathcal{H}=k+h_0+h_3+h_{\mathrm{eff}},\end{array}$$(45)

where k = ℏ²τ/2m is the kinetic energy term, h₀ is a density-independent two-body term, h₃ is a density-dependent term, and h_eff is a momentum-dependent term. Each of the interaction terms can be expressed analytically in terms of densities of particles and kinetic energies, and parameters of the effective interaction.

The effective interactions used here belong to the family of Brussels extended Skyrme interactions (Goriely et al. 2013), situation in which the density-independent two-body term and the density-dependent term have expressions identical to those obtained in the case of standard Skyrme interactions (Ducoin et al. 2006),

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

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

The situation of the momentum-dependent term is different. While keeping the structure it has for standard interactions,

$$\begin{array}{c}\hfill h_{\mathrm{eff}}={\stackrel{\sim }{C}}_{\mathrm{eff}}(n)n\tau +{\stackrel{\sim }{D}}_{\mathrm{eff}}(n)n_3{\tau }_3,\end{array}$$(48)

this term features a more complex density dependence due to the density dependence of ${\stackrel{\sim }{C}}_{\text{eff}}(n)$ and ${\stackrel{\sim }{D}}_{\text{eff}}(n)$. The latter relate to the C_eff and D_eff parameters of the momentum-dependent terms of the standard interactions (Ducoin et al. 2006) via (Beznogov & Raduta 2024)

$$\begin{array}{cc}\hfill {\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}$$(49)

In the equations above, n = n_n + n_p and n₃ = n_n − n_p stand for the isoscalar and isovector particle number densities; τ = τ_n + τ_p and τ₃ = τ_n − τ_p denote the isoscalar and isovector densities of kinetic energy; 2/m = 1/m_n + 1/m_p, where m_i with i = n, p denotes the bare mass of nucleons.

Our preference for Brussels extended interactions is due to the extra density-dependent terms of the momentum interaction, which makes it possible to heal a series of drawbacks commonly associated with the standard Skyrme interactions. In the first place, the increased flexibility of the effective Hamiltonian can be exploited to make the effective mass of the nucleons

$$\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}$$(50)

have a density dependence similar to the one put forward by ab initio models with three body forces (Shang et al. 2020; Somasundaram et al. 2021). In the limit of zero temperature, m_eff;i corresponds to the Landau effective mass 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}{k_i}\frac{d{e}_i}{d{k}_i}|}_{k=k_{F;i}},\end{array}$$(51)

where k_F; i stands for the Fermi momentum of the i-particle. More precisely, instead of continuously decreasing with density, as is the case of the standard Skyrme interactions, extended interactions allow m_eff;i to increase with density for densities exceeding a certain value. The bonus of banning m_eff;i from excessively decreasing with density is that the Fermi velocity of nucleons

$$\begin{array}{c}\hfill v_{F;i}=\frac{k_{F;i}}{{m_{\mathrm{eff}}}_{;i}},\end{array}$$(52)

in dense matter does not exceed the speed of light, which would make the interaction unphysical (Duan & Urban 2023).

The five effective interactions for which we build here general purpose EOS tables were extracted from the bunches of several hundred thousand interactions in runs 2 (BBSk1) and 1 (BBSk2 to BBSk5) in (Beznogov & Raduta 2024), built within a Bayesian inference of the EOS of dense matter. As every other effective interaction in (Beznogov & Raduta 2024), they comply by construction with standard constraints on: i) the best known parameters of NM, ii) ab initio calculations of energy per particle in pure neutron matter (PNM), (E/A)_PNM, with densities lower than n_sat, as computed by Somasundaram et al. (2021) iii) a 2 M_⊙ lower limit on the maximum gravitational mass of NS, iv) thermodynamic stability of NS EOS, v) causality of NS EOS up to the density corresponding to the central density of the maximum mass configuration as well as a number of extra constraints like: vi) for n ≤ 0.8 fm⁻³, the Landau effective mass m_eff of neutrons in PNM and of nucleons in symmetric NM (SNM) take values in the domain [0, m_N], where m_N is the bare mass, v) for n ≤ 0.8 fm⁻³, the neutron Fermi velocity v_F; n in both SNM and PNM does not exceed the speed of light. For BBSk1, which is our benchmark interaction, the effective mass of neutrons in PNM and the effective mass of nucleons in SMN also comply with ab initio calculations over 0.08 ≤ n ≤ 0.16 fm⁻³. Similarly to the values of (E/A)_PNM, we use m_eff predictions by Somasundaram et al. (2021). As PNM is stiffer than NS matter, we make sure that, for the forces selected in this work, PNM is causal at least up to 0.8 fm⁻³. This set of interactions was selected in such a way as to manifest the utmost variability in the behaviors of effective masses and thermal pressure as functions of density. BBSk1 corresponds to the “mean” model of run 2; BBSk2 and BBSk3 lie close to the 0.98 upper and 0.02 lower quantiles of m_eff(n) in run 1; BBSk4 and BBSk5 have extreme behaviors of thermal pressure.

The parameters of these interactions are provided in Table 1 while some of their nuclear properties are listed in Table 2 in terms of NM parameters. Also given in Table 2 are the values of Landau effective mass of nucleons in SNM and Landau effective mass of neutrons in PNM at n = 0.16 fm⁻³ as well as the values of the critical temperature of the liquid-gas phase transition in SNM.

3. Inhomogeneous to homogeneous matter transition

General purpose EOS tables on COMPOSE are tabulated as a function of temperature (T), baryonic charge fraction, here equal to the proton fraction (Y_p), and number density (n_B). This format commands to transition from inhomogeneous to homogeneous matter for every value of the grid temperature lower than the value of the Coulomb instability temperature of the nuclei (Jaqaman et al. 1984; Song & Su 1991) and Y_p.

Implementation of nucleus-nucleus and nucleus-nucleon interactions via the excluded volume suggests that this transition arises because either (I) the uniform gas of unbound nucleons occupies the entire volume or (II) the gas of clusters occupies the entire volume. A third mechanism consisting of a sudden replacement of the inhomogeneous phase by the homogeneous one (III) may be foreseen whenever the homogeneous phase becomes energetically favored, which corresponds to f_homo < f_inhomo where f stands for the free energy density.

Figure 1 investigates the behavior of three quantities, representative for each of the three transition mechanisms, as functions of density. These are the fraction of mass in the unbound nucleon gas, the relative volume occupied by nuclei, and the relative difference between the baryonic free energy densities of inhomogeneous and homogeneous matter. For better readability, the contributions of electrons and photons, which depend only on n_el and T, are removed from the free energy density.

Fig. 1. Selected properties of the mixture of unbound nucleons and nuclei while approaching the transition to homogeneous matter for T = 0.23 MeV (left column) and 5.75 MeV (right column) and various proton fractions, as indicated in the legend. Considered are the relative difference between the baryonic free energy densities of inhomogeneous and homogeneous matter (bottom row), intrinsic volume of nuclei with A ≥ 2, Vcl, relative to the full volume, V, (middle row) and relative density of unbound nucleons (top row) as functions of density. For better readability, the quantities in the bottom panel are multiplied with arbitrary factors shown in the legend. The considered interaction is BBSk1.

The case of low temperatures is addressed in the left panels, where T = 0.23 MeV and 0.03 ≤ Y_p ≤ 0.3. It comes out that for Y_p values exceeding a certain threshold, nuclear clusters bind a considerable fraction of matter. The increase in density results in larger fractions of volume being occupied by clusters. Nonetheless, the transition to homogeneous matter occurs before the clusters fill the entire volume, due to f_homo < f_inhomo. For lower values of Y_p, nuclear clusters bind a smaller fraction of matter. When the density exceeds a certain Y_p-dependent threshold, the gas of unbound nucleons suddenly becomes more dense, which forces a sudden diminish of the volume occupied by the clusters and of the mass fraction bound in the clusters. The sudden change in the density of the gas of unbound nucleons can be explained by the jump from the low-density stable branch of μ_p/n(n_p/n)|_μn/p curve(s) to either the intermediate-density unstable branch or the high-density stable branch. We remind at this point that for T ≤ T_C, where T_C = 15 − 20 MeV is the critical temperature of the liquid-gas phase transition, uncharged homogeneous matter exhibits thermodynamic instabilities with respect to density fluctuations, which are manifested as decreasing regions of the μ_p/n(n_p/n) curves at constant values of μ_n/p (Ducoin et al. 2006). Consequently, multiple solutions for (n_n,n_p) are possible for given (μ_n, μ_p)-sets that correspond to the phase instability window. The residual Coulomb interaction between the relatively few and light clusters in the clusterized phase is important enough to prevent f_inhomo from exceeding f_homo and thus being superseded by homogeneous matter. Indeed, for Y_p = 0.03 and 0.06 and high densities, f_inhomo/f_homo − 1 is almost constant and slightly lower than 0. It is clear that, under the conditions in the left panels, only mechanisms I and III are at work.

The right panels investigate what happens at higher temperatures, in particular T = 5.75 MeV. As before, several values of Y_p are considered. For the lowest considered value of proton fraction, Y_p = 0.18, the situation is similar to what we observed for the two lowest values of Y_p in the case of low T, that is, the transition occurs because the unbound nucleons kick out the clusters (I). For the second low value of Y_p, Y_p = 0.22, the transition occurs because homogeneous matter becomes energetically favored (III). For Y_p = 0.3, the transition occurs because nuclear clusters occupy the entire volume (II). For the highest considered Y_p the transition occurs again because homogeneous matter is energetically more favored (III).

A synoptic outline of the transition mechanisms as a function of Y_p and T is offered in Fig. 2 (top panel) for 0.01 ≤ Y_p ≤ 0.6 and T ≤ 10 MeV. The criteria upon which the transitions via mechanisms II and III are considered to take place are V_cl/V ≥ 0.9 and f_inhomo/f_homo − 1 ≥ −10⁻⁶, respectively; the value of −10⁻⁶ arises from the precision with which the values of the free energy density of inhomogeneous and homogeneous matter are saved into the files. The difficulty in achieving convergence of the eNSE equations when the gas of strongly interacting unbound nucleons is unstable or very dense made it preferable to consider that the transition via the mechanism I occurs at the density where V_cl(n)|_Yp reaches the maximum. Figure 2 shows that, for low values of Y_p, matter transitions via mechanism I; for 5 ≲ T ≲ 10 MeV and Y_p ≳ 0.25, matter transitions via mechanism II; in all other cases, the transition is done by supersession (III). Some arbitrariness of the transition criteria is reflected in narrow domains bearing the color of mechanism III in a region where the transition holds via mechanism I. There is a single set of (T, Y_p) (marked in red) where difficulties in solving the eNSE equations made it impossible to reach a density high enough to see one of our transition mechanisms in operation. In this situation, the transition is done after the last converged density point.

Fig. 2. Mechanisms (top panel) and the density (bottom panel) for the transition to homogeneous matter as a function of temperature and proton number. The color legend in the top panel is as follows. Mechanism I (blue): ntr is taken as nB where Vcl(n) reaches the maximum for fixed values of T and Yp; mechanism II (green): ntr corresponds to the lowest value of nB where Vcl/V ≥ 0.9; mechanism III (cyan): ntr corresponds to the lowest value of nB where finhomo/fhomo − 1 ≥ −10−6. Cases where the transition is considered to occur at the highest density where convergence of Eqs. (43) and (44) is reached are marked in red. In the bottom panel the color levels correspond to ntr in units of fm−3. The considered interaction is BBSk1.

The transition densities are indicated in the bottom panel in Fig. 2. The same domains of Y_p and T as in the top panel are considered. It comes out that $0.05\lesssim n_{\mathrm{tr}}\lesssim 0.14{\mathrm{fm}}^{-3}$ and there is a significant T- and Y_p-dependence. The transition via mechanism II gives the highest values of n_tr. In cold matter, the transition density varies between n_sat/3 and n_sat/2, depending on Y_p. A question that arises now is how close these transition densities are to those predicted by more sophisticated models, which implement in-medium and temperature modifications of nuclear cluster functionals. In Sect. 5.1 we show that for cold NS matter our results are in fair agreement with those of microscopic models.

Determination of the subsaturation density where clusterized stellar matter is superseded by homogeneous matter is not enough to build the kind of EOS tables that numerical simulations of CCSN and BNS mergers require. The reason is that the Coulomb interaction present in the clusterized phase and absent in the homogeneous one results in the neutron and proton chemical potentials in one phase and the other being offset, see Eq. (18). A thermodynamically stable matching between one phase and the other, both having the same Y_p value, can be done by identifying the low density counterpart of the homogeneous matter just after the transition and smoothly merging them. The quantity we choose for identifying this low-density counterpart is the pressure. This smooth matching consists in replacing all eNSE values in between with the results of the linear interpolation between the eNSE solution at low density and its homogeneous counterpart at high density. A similar solution was adopted by Hempel & Schaffner-Bielich (2010). The composition of matter in this transition domain is derived by linearly diminishing the abundances of all light and heavy nuclei that exist at the low density border, without altering their mass and charge; the amount of unbound nucleons results from mass and charge conservation.

This approach may remind one of the Maxwell construction. We insist that it is only a smooth matching, and there are two arguments against the interpretation of this approach as a Maxwell construction. First, a Maxwell construction assumes that a phase transition takes place but, as is demonstrated in Sect. 4, clusterized matter with electrons is stable. Second, two or more phases can be put in equilibrium if all the intensive variables are equal. This is not the case here, as only the pressure has the same value in one phase and in the other.

4. Thermodynamic stability of stellar matter

Relativistic (Müller & Serot 1995; Liu et al. 2002; Providencia et al. 2006; Avancini et al. 2006) and non-relativistic (Ducoin et al. 2006, 2007; Rios 2010) mean-field calculations as well as ab initio approaches (Vidana & Polls 2008; Rios et al. 2008; Carbone et al. 2018) demonstrate that over a large density domain dilute homogeneous NM without Coulomb interactions manifests density-driven instabilities. According to Typel et al. (2010), Raduta & Gulminelli (2009, 2010), the same is the case of dilute clusterized NM.

The thermodynamic stability of clusterized NM is addressed in Fig. 3 (left panel) in terms of neutron chemical potential μ_n as a function of neutron density n_n for constant values of proton chemical potential μ_p (Ducoin et al. 2006), magenta curves. For simplicity, only the case of T = 0.1 MeV is addressed. For all considered values of μ_p and over the entire illustrated density domain, μ_n(n_n) decreases, which indicates that the system is thermodynamically unstable. The set of green curves in Fig. 3 corresponds to clusterized stellar matter, that is, a mixture of nuclei, nucleons, electrons, and photons with a positive charge of protons fully compensated by the negative charge of electrons. In this case, the baryon chemical potential μ_B(=μ_n) is plotted as a function of the baryon density n_B for constant values of the lepton chemical potential μ_L (Ducoin et al. 2007). The curve with μ_L = 0 corresponds to beta-equilibrated matter. None of the μ_B(n_B)|_μL-curves decreases with density, which means that no instability occurs in that domain of density and isospin asymmetry. This demonstrates that the necessity for Maxwell constructions does not arise at T = 0.1 MeV. For homogeneous NM, the phase instability domain shrinks with temperature; the same holds true for clusterized NM without Coulomb (Raduta & Gulminelli 2010). Considering that cold clusterized stellar matter does not feature any instability, it is reasonable to assume that the matter is stable at higher temperatures, too.

Fig. 3. Thermodynamic (in)stability of cold clusterized matter. Green curves represent the properties of clusterized stellar matter with fixed lepton chemical potentials (expressed in MeV) in terms of the baryon chemical potential, μB(=μn), as a function of baryonic number density, nB, (left panel) and total pressure, P, (right panel). Magenta curves represent the properties of clusterized baryonic matter with fixed proton chemical potentials (expressed in MeV) in terms of neutron chemical potential, μn, as a function of neutron number density, nn, (left panel) and baryonic pressure, PB, (right panel). The considered interaction is BBSk1. Here, μn and μp do not contain rest-mass contributions.

Additional information is presented on the right panel, which shows the curves of nuclear pressure P_B vs. μ_n for constant values of μ_p (magenta) and the curves of total pressure P vs. μ_B for constant values of μ_L (green). P_B is obtained from the total pressure P stored in COMPOSE tables after subtracting the electron and photon contributions. The magenta (green) family of curves resembles the corresponding curves of homogeneous nuclear (stellar) matter, i.e., are convex (concave) whenever matter is unstable (stable) with respect to phase separation.

We note that the conclusions of this section are not incompatible with subsaturated NM with electrons being unstable with respect to finite size density fluctuations (Baym et al. 1971; Pethick et al. 1995; Ducoin et al. 2007). On the contrary, present results demonstrate that after clusterization due to periodic microscopic instabilities the matter becomes stable.

5. Results

This section is devoted to the results of our eNSE model. First, we investigate the composition of the NS crusts and its dependence on the effective interaction. A systematic comparison with the predictions of the extended Thomas-Fermi plus Strutinsky integral (ETFSI) method by Pearson et al. (2018) is also carried out. Then, we turn to the finite-temperature regime. The composition and energetics of subsaturated matter are investigated for a wide range of temperatures and proton fractions.

5.1. Neutron star crusts

Here, we investigate the predictions of our eNSE model in the limit of low temperatures. Because T = 0.1 MeV considered in this subsection falls outside the validity domain of the nuclear statistical equilibrium, some explanations are necessary.

The reason why we expect NSE to give a fair description of the composition is that, as the temperature approaches zero, NSE distributions are practically reduced to a single nucleus (Hempel & Schaffner-Bielich 2010; Gulminelli & Raduta 2015). Then, an NSE state is by construction a state of equilibrium. In the limit of vanishing temperatures, that state corresponds to the minimum energy of the system. If the controlled proton fraction is tuned to comply with μ_n = μ_p + μ_e, NSE should provide a good approximation for the NS crust.

The composition of cold β-equilibrated matter is addressed in Fig. 4 in terms of neutron and proton numbers of the nuclear species present at various depths inside the NS crust. The effective interactions considered are BSk22 and BSk25 (Goriely et al. 2013). It turns out that the crust composition is stable over density domains of various widths and when it gets modified the changes are drastic. This is the well-known outcome of shell and odd-even effects in nuclear masses. Indeed, one can recognize the Z = 28 and N = 50,  82 magic numbers along with many even-Z numbers in the inner crust. The constraints of the AME2020 (Wang et al. 2021) and DZ10 (Duflo & Zuker 1995) nuclear mass tables, which we employ and which limit the crust composition, are such that the size of the cluster does not always increase with density. In particular, over 5 ⋅ 10⁻⁷ fm⁻³ ≲ n_B ≲ 8 ⋅ 10⁻⁵ fm⁻³ and 4 ⋅ 10⁻⁵ fm⁻³ ≲ n_B ≲ 10⁻² fm⁻³ nuclei in the crust have N = 50 and N = 82 while their charge number decreases from 36 to 28 and from 42 to 28, respectively. Over 6.3 ⋅ 10⁻⁷ fm⁻³ ≲ n_B ≲ 2.5 ⋅ 10⁻⁶ fm⁻³ and 4 ⋅ 10⁻⁵ fm⁻³ ≲ n_B ≲ 8 ⋅ 10⁻⁵ fm⁻³ two species with significantly different Z- and N-numbers are in competition, and, for some densities, neighboring nuclides are present, too. The novel feature predicted by our model is the occurrence of a thick ¹⁴He shell at the bottom of the crust. This nuclide is present in the DZ10 table (Duflo & Zuker 1995); its one- and two-neutron separation energies are −1.29 MeV and −5.47 MeV, respectively. The occurrence of light neutron-rich nuclides in the innermost layers of the crust was checked by limiting the pool of NSE nuclei to those present in the AME2020 (Wang et al. 2021) table only. The results indicate that, in this case, a layer of ⁷H is present instead of ¹⁴He, see the discussion below.

Fig. 4. Composition of an NS crust as a function of baryon number density in terms of numbers of neutrons (bottom panel) and protons (top panel) of nuclei. The eNSE results obtained in this work for β-equilibrated matter with T = 0.1 MeV are compared with the results of Pearson et al. (2018). The considered effective interactions are BSk22 and BSk25 (Goriely et al. 2013). We note that BSk25 results of Pearson et al. (2018) are scaled by a factor of 1/5. The vertical dashed lines respectively mark the transition to homogeneous matter for BSk22, ntr = 5.70 ⋅ 10−2 fm−3 and 6.35 ⋅ 10−2 fm−3, for the eNSE results and the results of Pearson et al. (2018). Horizontal dotted lines correspond to magic numbers: Z = 8, 20, 28 and N = 8, 20, 28, 50, 82.

Even if non-spherical degrees of freedom are not explicitly implemented in our work, we can speculate that the onset of light neutron-rich nuclides might challenge the pasta phases, should the model account for them. A detailed study of the consequences of such layers on the mechanical and transport properties of the crust is beyond the scope of this article. For the issue of mechanical stability, see Beznogov & Raduta (2025); some comments on the impurity parameter are given below. Effective interaction effects manifest in the inner crust, but the representation in Fig. 4 does not allow one to quantify their amplitude. It is nevertheless clear that much stronger effective interaction effects would manifest if larger and more neutron-rich nuclei are allowed to exist and their binding energies are allowed to reflect current uncertainties in the isovector channel.

For the sake of comparison, the results of ETFSI calculations with the BSk22 and BSk25 interactions (Goriely et al. 2013) by Pearson et al. (2018) are plotted, too. For BSk22 (BSk25) they show that up to n_B = 3.4 ⋅ 10⁻² fm⁻³ (n_B = 1.4 ⋅ 10⁻² fm⁻³) only Zr (Sn) is present in the inner crust and its neutron enrichment increases steeply and smoothly with density; for BSk22, a Ca-layer exists over 3.5 ⋅ 10⁻² fm⁻³ ≲ n_B ≲ 5.3 ⋅ 10⁻² fm⁻³; shell effects manifest only in the charge number. The size of the nuclear cluster predicted by ETFSI is much larger than the one we obtain.

The vertical dashed lines mark the crust-core transition densities for BSk22. In spite of performing calculations on a grid and limiting the crust composition to species present in the mass tables, the value we obtain for the crust-core transition density is remarkably close to the value obtained by Pearson et al. (2018).

More insight into the inner crust composition is provided in Fig. 5. The top panel depicts the mass fractions of nuclei and unbound nucleons as functions of number density. The contributions of ¹⁴He are indicated with dot-dashed curves. Despite the stark differences with the predictions of Pearson et al. (2018) in regard to the properties of nuclei, the mass sharing between the nuclei and nucleons is fairly similar in our model and the ETFSI model of Pearson et al. (2018), at least until the emergence of ¹⁴He. This suggests that the gross behavior in the inner crust is governed by the effective interaction. The nucleation of ¹⁴He in the deep crust entails a steep depletion of the gas of unbound nucleons with possible effects on ¹S₀ superfluidity. The effective interaction dependence of the inner crust composition is further investigated in the bottom panel, where, in addition to the results corresponding to BSk22 and BSk25, we also include the results corresponding to BBSk1 and BBSk2 forces. Results corresponding to the case where the NSE pool of nuclei is limited to AME2020 table only (eNSE^*) are illustrated, too. One can see that ⁷H has the same onset density and slightly smaller mass fraction (≈14% vs. ≈20%) compared to ¹⁴He. The onset density of ¹⁴He and the relative mass sharing between ¹⁴He and the nucleon gas depend on the underlying interaction. The same holds for the width of the ¹⁴He layer and the transition density to homogeneous matter.

Fig. 5. Mass fractions of various species present in the inner crust and their dependence on the effective interaction as functions of baryon number density. The eNSE results are for T = 0.1 MeV. Top panel: eNSE predictions confronted with the zero-temperature results of Pearson et al. (2018). Bottom panel: Mass fractions of 14He, 7H and unbound nucleons, as predicted by eNSE. The considered effective interactions are BSk22 and BSk25 from (Goriely et al. 2013) and BBSk1 and BBSk2 proposed here. As in Fig. 4, vertical dashed lines mark the transition to homogeneous matter.

The crust of an isolated NS forms during the neo-NS evolutionary phase (Beznogov et al. 2020). This stage begins at the end of the proto-NS epoch, when the star’s temperature drops below ≈2 MeV and its material becomes transparent to neutrinos, and lasts for about a day. Two episodes are of particular interest for the matters discussed here. One corresponds to the instance when NSE freezes out and marks the beginning of the epoch when only a limited subset of nuclear reactions can change the composition of matter. The other corresponds to the crystallization of the crust. The freeze-out temperature of NSE depends on the density and, even if further composition changes are possible, its value plays an important role for crust composition and thermal properties. The crystallization temperature of a Coulomb plasma depends on the ion’s charge squared (Chamel & Haensel 2008). This means that the inner crust of BSk22 by Pearson et al. (2018) will have a crystallization temperature higher by a factor of one hundred compared to our inner crust with ¹⁴He layer. In the absence of dynamical simulations of crust formation, it is difficult to say which of these instances occurs first and what exactly the crust composition and transport properties are. As such, only statistical estimates can be made. A useful quantity to this aim is the impurity parameter (Potekhin & Chabrier 2021),

$$\begin{array}{c}\hfill Q_{\mathrm{imp}}={\displaystyle \sum _j}{y}_j{(Z_j-⟨Z⟩)}^2,\end{array}$$(53)

where ⟨Z⟩=∑_jy_jZ_j is the mean charge number, y_j = Y_j/∑_iY_i and Y_i stands for the yield per unit volume of species i.

Figure 6 (bottom panel) illustrates the density dependence of Q_imp for 0.1 MeV ≤ T ≤ 0.58 MeV. For the lowest temperature, Q_imp ≈ 0 only in the highly pure ¹⁴He layer. Over 10⁻⁹ fm⁻³ ≲ n_B ≲ 10⁻⁶ fm⁻³, Q_imp changes smoothly between ≈0.3 and ≈3; these values reflect the non-pure but rather stable composition, see Fig. 4. Over 10⁻⁴ fm⁻³ ≲ n_B ≲ 10⁻² fm⁻³, Q_imp fluctuates widely over fourteen orders of magnitude; this behavior reflects sudden and drastic changes in composition. With the increase of temperature, more species are present and the composition changes less drastically with density. Consequently, the Q_imp oscillations at high densities are dumped, while some oscillations appear at low densities; the lowest values correspond to the ¹⁴He layer. The upper panel addresses the effective interaction dependence of this quantity in the inner crust. It comes out that while the amplitude of the oscillations seems to be preserved, some offset (“phase shift”) exists.

Fig. 6. Impurity parameter Qimp in the NS crust. Bottom panel: Predictions of BBSk1 at various temperatures. Top panel: Effective interaction effects at T = 0.1 MeV.

While the oscillations of Q_imp are well known in the literature (see, e.g., Carreau et al. 2020b; Fantina et al. 2020; Potekhin & Chabrier 2021), our results are somewhat different. We have smoother behavior in the outer crust and strong oscillations in the outer layers of the inner crust. However, in the deeper layers of the inner crust, after the onset of ¹⁴He, Q_imp becomes almost density independent, which is not the case of other models (Carreau et al. 2020b; Potekhin & Chabrier 2021). Moreover, Q_imp becomes independent of the chosen effective interaction, which is in strong contrast with the results of Potekhin & Chabrier (2021). Thus, if the existence of the almost pure ¹⁴He or ⁷H layer is proven, it will be a step toward a universal (i.e., EOS independent) description of the transport properties of the deep layers of the inner crusts.

5.2. Finite temperature results

This subsection presents an overview of the composition and energetics of warm and dilute NM. A series of quantities corresponding to a fixed value of temperature (proton fraction) and different values of the proton fraction (temperature) are plotted as functions of density. The particular values of these quantities have been chosen to highlight the most interesting features. Also, they correspond to mesh points, which keeps the results free from inherent interpolation-related effects.

5.2.1. Composition

In NSE, the clusterized component of NM accommodates a wealth of nuclei, whose number virtually equals the number of nuclei included in the mass tables. Some of these nuclei have negligible yields, while others are abundant. The number of nuclei that appear in significant amounts depends on the thermodynamic conditions and intrinsic nuclear properties. Both aspects were extensively discussed by Hempel & Schaffner-Bielich (2010), Raduta & Gulminelli (2010, 2019), Gulminelli & Raduta (2015), Blinnikov et al. (2011), Furusawa et al. (2011, 2013, 2017a,b) and are not repeated here. The only aspect we wish to remind is that whenever the nuclear distributions are multi-modal, which is typically the case at relatively low temperatures and proton fractions exceeding a certain threshold, the average numbers of neutrons and protons of nuclei in a category might not be representative of any of the nuclei that appear in large amount.

With this caveat, the composition of NM is discussed here in terms of average mass and proton numbers,

$$\begin{array}{c}\hfill ⟨A_{\mathcal{C}}⟩=\frac{{\displaystyle {\sum }_{i\in \mathcal{C}}}{A}_i{Y}_i}{{\displaystyle {\sum }_{i\in \mathcal{C}}}{Y}_i},⟨Z_{\mathcal{C}}⟩=\frac{{\displaystyle {\sum }_{i\in \mathcal{C}}}{Z}_i{Y}_i}{{\displaystyle {\sum }_{i\in \mathcal{C}}}{Y}_i},\end{array}$$(54)

and mass fraction,

$$\begin{array}{c}\hfill X_{\mathcal{C}}=\frac{{\displaystyle {\sum }_{i\in \mathcal{C}}}{A}_i{Y}_i}{n_B},\end{array}$$(55)

of two generic species for which relative mass abundance is reported in the tables contributed to COMPOSE, i.e., “heavy” nuclei ($\mathcal{C}=\{A,Z\}$ with A ≥ 20) and “light” nuclei ($\mathcal{C}=\{A,Z\}$ with Z ≥ 3 and A < 20). As is Sect. 5.1, in the equations above, Y_i stands for the yield per unit volume of species i. The properties of the gas of unbound nucleons as well as the abundances of nuclides with Z = 1 and 2 and A ≥ 2, also reported in the tables contributed to COMPOSE, are addressed, too.

Figure 7 reports the average mass number and the ⟨Z⟩/⟨A⟩ value of “heavy” nuclei produced in neutron rich matter (Y_p = 0.1) with 0.1 ≤ T ≤ 5.75 MeV along with the mass fraction which this category binds. Shell effects are visible up to T = 0.58 MeV. As already discussed by Hempel & Schaffner-Bielich (2010), Gulminelli & Raduta (2015), Raduta & Gulminelli (2019) among others, “heavy” clusters’ size decreases with temperature and increases with density; deviations from the latter trend appear only at low temperatures and high densities because of structure effects. X_heavy features the same behavior as ⟨A_heavy⟩. At high densities and low (high) temperatures, its values are about 40% (10%). As the temperature increases, also ⟨Z_heavy⟩/⟨A_heavy⟩ increases; this is due to increasing abundance of lighter and more isospin symmetric species. Figure 8 illustrates the average mass and charge numbers of the nuclei with Z ≥ 3 and A < 20 as well as their mass fraction. It turns out that these species only appear for T ≳ 0.5 MeV. As the temperature increases, the density domain where they are present shrinks and moves toward higher densities, while X_light increases. Nonetheless, for all thermodynamic conditions considered here, X_light never exceeds a few percent. The suppression of light clusters at n_B = 10⁻⁶ fm⁻³ and T = 0.58 MeV is due to the onset of extra H and He isotopes. Abundant production at high densities of an increasing number of H and He isotopes, which compete with light clusters, also explains the double-peaked distribution of X_light for T = 3.63 and 5.75 MeV.

Fig. 7. Average mass numbers of nuclei with A ≥ 20, the corresponding ⟨Z⟩/⟨A⟩ ratio, and the mass fraction they bind (X) as functions of baryonic particle density for Yp = 0.1 and different temperatures. The considered effective interaction is BBSk1.

Fig. 8. Average mass and charge numbers of nuclei with Z ≥ 3 and A < 20 as well as the mass fraction they bind as functions of baryonic particle density for the cases considered in Fig. 7. We note that the X-range is not the same as in Fig. 7.

Figure 9 demonstrates the same quantities plotted in Fig. 7 but for T = 1 MeV and 0.01 ≤ Y_p ≤ 0.6. For all Y_p values, nuclear shell and odd-even effects make that ⟨A_heavy⟩ as a function of n_B presents a sequence of plateaus. The same is true for ⟨Z_heavy⟩ (not shown). This structure appears most clearly for Y_p = 0.5. For Y_p = 0.4, ⟨A_heavy⟩, ⟨Z_heavy⟩ and X_heavy reach the largest values; in particular, over 10⁻⁴ fm⁻³ ≲ n_B ≤ n_tr ≈ 6.7 ⋅ 10⁻² fm⁻³, X_heavy ≈ 1. For n_B ≲ 10⁻⁶ fm⁻³ and 0.1 ≤ Y_p ≤ 0.6, the clusters are almost symmetric; this is due to the fact that they are relatively light. For higher n_B, ⟨Z_heavy⟩/⟨A_heavy⟩ decreases with n_B and increases with Y_p. The former feature is due to the fact that ⟨A_heavy⟩ increases; the latter is the result of more protons being available. Remarkably, for high n_B and 0.4 ≤ Y_p ≤ 0.6, ⟨Z_heavy⟩/⟨A_heavy⟩≈Y_p. Figure 10 shows that the largest values of ⟨A_light⟩≈17 are obtained for 0.01 ≤ Y_p ≤ 0.25 and n_l ≲ n_B ≲ n_tr; for Y_p = 0.01 (0.25), n_l = 2.5 ⋅ 10⁻³ fm⁻³ (1.6 ⋅ 10⁻² fm⁻³). The largest values of ⟨Z_light⟩≈6 are obtained around n_B = 10⁻⁵ fm⁻³ for Y_p = 0.5,  0.6. Structure effects and/or competition among species lead to irregularities in the n_B-dependence of these quantities. Similarly to the situations considered in the left panel, X_light ≲ 2%; the highest values correspond to Y_p = 0.6.

Fig. 9. Same as in Fig. 7 but for T = 1 MeV and different values of Yp.

Fig. 10. Same as in Fig. 8 but for the cases considered in Fig. 9. We note that the X-range is not the same as in Fig. 9.

The relative abundance of unbound nucleons and their proton fraction are investigated in Figs. 11 and 12. Figure 11 shows that both n_gas/n_B and the density domain where the gas of unbound nucleons dominates increase with temperature. The onset of nuclei causes the gas to become depleted in protons. The transition to homogeneous matter manifests as a steep increase in both n_gas(n_B)/n_B and Y_p(n_B); for the considered value of Y_p, n_tr is independent of T. Figure 12 shows that for T = 1 MeV the onset density of the nuclear clusters does not depend on Y_p; the same holds for higher values of T. In contrast, the density for the transition to homogeneous matter depends on Y_p (even if this is barely visible on the plot). Almost pure neutron and proton gases exist over 10⁻⁵ fm⁻³ ≤ n_B ≤ n_tr for Y_p ≤ 0.4 and Y_p ≥ 0.5, respectively, but for 0.4 ≤ Y_p ≤ 0.5, the unbound nucleon gas represents only a tiny fraction of matter.

Fig. 11. Mass fraction of unbound nucleons (bottom panels) and their proton fraction (top panels) as functions of baryonic particle density for the cases considered in Fig. 7. We note that the X-range is not the same as in Fig. 7.

Fig. 12. Same as in Fig. 11 but for the cases considered in Fig. 9. We note that the X-range is not the same as in Fig. 9.

The case of hydrogen and helium isotopes, not included in the “light” nuclei, is analyzed in Fig. 13 for T = 1.44 MeV. The top panel illustrates the n_B − Y_p domains where the mass fraction of different hydrogen isotopes exceeds 0.2% as well as the values of these mass fractions. It comes out that only three species are populated abundantly enough: ²H, ³H and the very exotic ⁷H. Significant production of the isospin symmetric ²H occurs in matter with n_B ≈ 2 ⋅ 10⁻⁶ fm⁻³ and Y_p ≳ 0.25; the more neutron-rich ³H requires densities higher by a factor of ten and Y_p ≲ 0.35; neutron-rich dense matter favors the production of extremely neutron-rich ⁷H, which binds up to 20 times more mass than each of the lighter isotopes. The bottom panel illustrates the n_B − Y_p domains where the mass fractions of neutron-rich isotopes of helium exceed 1%. It turns out that dense neutron-rich matter also favors the production of ¹¹He, ¹²He, ¹³He, ¹⁴He. ¹⁴He (¹¹He) is present over the widest (narrowest) n_B − Y_p domain. The highest values of X_He obtained under these circumstances are the following: 1.5% (¹¹He), 3.3% (¹²He), 9.4% (¹³He) and 71.2% (¹⁴He).

Fig. 13. Domains of nB − Yp where at T = 1.44 MeV the mass fractions of H isotopes (top panel) and 11He (red), 12He (blue), 13He (green), 14He (cyan) (bottom panel) exceed 0.2% and 1%, respectively. The color scale in the top panel corresponds to the mass fraction bound in each of the considered H isotopes. The effective interaction is BBSk1.

Confirmation of abundant production of light exotic nuclei in dense matter requires a more elaborated treatment of in-medium effects (i.e., beyond excluded volume). So far only the case of ²H, ³H, ³He and ⁴He in symmetric matter has been microscopically addressed (Typel et al. 2010), which is not sufficient for a comprehensive understanding of stellar matter. However, Hempel et al. (2011) suggest that, whenever light clusters are abundant, accurate treatment of their energy shifts results in yields larger than those provided by the excluded volume approximation implemented in eNSE.

5.2.2. Energetics

The T-, Y_p-, and n_B-dependencies of the energetic state variables are addressed here as functions of density for constant values of T (Y_p) and different values of Y_p (T). Only the neutron and proton chemical potentials, (scaled) baryonic pressure, and baryonic energy per particle are considered. The latter two quantities are obtained upon subtraction of electron and photon gas contributions from the total energy density and total pressure and are preferred for better readability. We note that the Coulomb interaction among nuclei and electrons, as well as the one among electrons, is not affected by this subtraction.

Figures 14 and 15 deal with the T- and n_B-dependencies of matter with Y_p = 0.3. For T ≤ 0.58 MeV and n_B ≲ 10⁻² fm⁻³, both chemical potentials have low values and show little n_B-dependence. This behavior suggests that the unbound nucleon gas is dilute and that much of the density variation of μ_MF; n/p⁰ is annihilated by the cluster-induced offset $T{v}_0{\displaystyle \sum _{A,Z}}{N}_{A,Z}/V_{\mathrm{cl}}^{\mathrm{free}}$, see Eq. (18). For the two lowest values of T, the first conclusion is corroborated by the values of the baryonic energy per particle, e_B/n_B, which is consistent with a large fraction of matter bound in clusters, as well as its flatness as a function of n_B. In fact, the value of e_B/n_B − m_n ≈ −7 MeV is due to the binding energy of nuclei, the other terms in Eq. (37) being negligible. For T = 0.58 MeV, e_B(n_B)/n_B features a more complex structure that signals a transition from a gas-dominated composition at low densities, with e_B/n_B − m_n ≈ 0.5 MeV, to a cluster-dominated composition at n_B = 10⁻⁶ fm⁻³, where the values of e_B/n_B are very similar to those obtained at lower temperatures. For T ≥ 1.45 MeV and n_B ≲ 10⁻⁵ fm⁻³, μ_n and μ_p manifest strong T- and n_B-dependencies; this indicates that the gas component outsizes the clusterized one. Indeed, e_B/n_B − m_n is positive. For 10⁻³ fm⁻³ ≲ n_B ≲ 10⁻¹ fm⁻³, μ_p(n_B) decreases steeply at all temperatures. The same is the case of e_B(n_B)/n_B. This is the domain preceding the transition to homogeneous matter, where most of the matter is bound in clusters. Indeed, for T ≤ 3.63 MeV, e_B/n_B − m_n < 0. In the proximity of the transition to homogeneous matter, μ_n(n_B) presents a sudden drop of ≈5 MeV. Contrary to what happens in the case of μ_p(n_B), this is the artifact of the matching between the inhomogeneous and homogeneous phases, see Sect. 3. Higher densities correspond to the homogeneous matter regime; there, μ_n, p(n_B) increase steeply and show limited sensitivity to temperature.

Fig. 14. Neutron and proton chemical potentials as functions of baryonic particle density for Yp = 0.3 and 0.1 MeV ≤ T ≤ 5.75 MeV.

Fig. 15. PB/(nBT) and energy per particle (E/A = eB/nB) as functions of baryonic particle density for the same conditions as in Fig. 14. For better readability, the contributions of electron and photon gases have been removed.

With the exception of T = 0.58 MeV, e_B(n_B)/n_B and P_B(n_B)/(n_BT) feature wide plateaus over 10⁻¹² fm⁻³ ≲ n_B ≲ 10⁻⁶ − 10⁻³ fm⁻³. In the case of low-T, this result corresponds to clusterized matter, but the plot does not allow one to understand whether the ideal gas regime is reached or not. For T ≥ 1.45 MeV and low densities, P_B/(n_BT) = 1, which corresponds to the ideal gas of unbound nucleons. For higher densities, the behavior of P_B(n_B)/(n_BT) versus n_B is very similar to that of e_B(n_B)/n_B versus n_B. We note that, due to Coulomb, for low T and high n_B, P_B < 0, in agreement with Hempel & Schaffner-Bielich (2010), Raduta & Gulminelli (2019). Note that for both e_B(n_B)/n_B and P_B(n_B)/(n_BT) only the behavior of clusterized matter before the matching with homogeneous matter is depicted. The reason for dropping the higher density region resides in the difficulties of subtracting the electron contributions from the total quantities stored in the tables.

Complementary information is provided in Figs. 16 and 17, where we illustrate the Y_p-dependence of μ_n, p, e_B/n_B and P_B/(n_BT) as functions of n_B at T = 1 MeV. For Y_p = 0.01, where less than 0.5% of the matter is bound in clusters, see Fig. 9, the plateaus of e_B/n_B and P_B/(n_BT) as functions of n_B extend over the widest density domain, and both μ_n(n_B) and μ_p(n_B) are smooth. It turns out that the extremely neutron-rich NM behaves like an ideal gas up to n_B ≈ 10⁻³ fm⁻³. For 0.4 ≤ Y_p ≤ 0.6, where “heavy” clusters bind the largest fraction of matter, see Fig. 9, e_B/n_B and p_B/(n_BT) reach the lowest negative values. For Y_p ≥ 0.1, deviations from the ideal gas behavior start to manifest at n_B = 10⁻⁸ fm⁻³. As before, the sharp drops of μ_n/p around n_B ≈ 6 ⋅ 10⁻² fm⁻³ are artifacts of our matching procedure.

Fig. 16. Same as in Fig. 14 but for T = 1 MeV and 0.01 ≤ Yp ≤ 0.6.

Fig. 17. Same as in Fig. 15 but for the conditions considered in Fig. 16.

The EOS tables proposed in this paper use the same pool of nuclei and the same definition of clusters. As such, the model dependence between them is only due to different effective interactions. This dependence increases with the isospin asymmetry and barely manifests before the transition to homogeneous matter. It affects mostly the transition density and the composition of matter. In general, the main source of model dependence comes from the pool of nuclei and the definition of clusters in terms of binding energies, level density, maximum excitation energy, interactions with the medium, evolution of the shell effects with departure from the stability valley, and temperature. These aspects have been thoroughly addressed in the literature (Hempel et al. 2011; Buyukcizmeci 2013; Furusawa et al. 2013, 2017a,b; Raduta et al. 2021) and are beyond the scope of the present work. The behavior of suprasaturated matter, which is of utmost importance for the evolution of astrophysical phenomena, will be discussed elsewhere.

6. Outlook

In this paper, we have proposed a new set of general purpose EOS tables for simulations of CCSN and BNS mergers. The tables are available on COMPOSE, and they were derived within the non-relativistic mean-field model of NM and employ Brussels extended Skyrme interactions previously generated within a Bayesian inference of the EOS of dense matter (Beznogov & Raduta 2024). The underlying effective interactions were selected in such a way as to manifest the utmost variability in the behavior of Landau effective masses in the suprasaturation domain.

An eNSE model was developed and used to describe the composition and energetics of inhomogeneous matter present at subsaturation densities and for temperatures of up to a dozen mega electron volt. The pool of nuclei was represented by the nuclides present in the AME2020 (Wang et al. 2021) and DZ10 (Duflo & Zuker 1995) tables, from which the binding energies were taken. Internal partition functions were computed according to the back-shifted Fermi gas level density (Bucurescu & von Egidy 2005; von Egidy & Bucurescu 2005) while assuming that the upper limit of the excitation energy is equal to the binding energy. Interactions among nuclei and with the unbound nucleons were treated within the geometrical excluded-volume approximation. The isospin dependence of the intrinsic volumes of the clusters was implemented via the isospin dependence of the saturation density of the NM.

The highly efficient computer code developed to solve the eNSE equations allowed us to study how the transition to homogeneous matter occurs. As such, three mechanisms were identified. Two of them are possible due to the excluded volume approximation and correspond to the situations in which one of the two components, i.e., the ensemble of nuclei or unbound nucleons, fills the whole available volume. The third mechanism corresponds to the case where the homogeneous matter becomes energetically more favorable. The transition density values range from ≈n_sat/3 to almost n_sat, depending on temperature and Y_p. The highest values were obtained for temperatures close to the limiting temperature of the Coulomb instability and 0.4 ≲ Y_p ≲ 0.6, which roughly correspond to the situations where the transition occurs because the clusters fill the entire volume.

The ability of our eNSE code to approach n_sat made it possible to revisit the issue of thermodynamic instabilities related to the transition from inhomogeneous to homogeneous matter. Given that the domain of phase instability with respect to density fluctuations in NM shrinks with temperature, only the case of cold matter (T = 0.1 MeV) was considered. Our results indicate that no phase instability exists in clusterized matter with electrons.

The composition and energetics of inhomogeneous NM were systematically investigated over wide domains of temperature, density, and proton fraction. As expected, the same general features previously commented on by Hempel & Schaffner-Bielich (2010), Blinnikov et al. (2011), Furusawa et al. (2011, 2013, 2017a,b), Gulminelli & Raduta (2015), Raduta & Gulminelli (2019) were obtained. An original result here is the abundant production of neutron-rich isotopes of H and He in neutron-rich environments at finite temperatures, which is due to the inclusion of these species into the NSE pool of nuclei. Relying on the fact that the configurations, which minimize the densities of free energy and energy, coincide in the limit of vanishing temperatures, we also considered the composition of β-equilibrated matter at T = 0.1 MeV as a proxy for cold NS matter, even if this temperature is beyond the validity limit of NSE. We find that a relatively thick layer made of light neutron-rich nuclides, such as ¹⁴He or ⁷H, exists at the bottom of the inner crust. These layers are expected to challenge the occurrence of nuclear pasta phases and affect the transport properties of the inner crusts as well as the process of crystallization of NS crusts. However, it remains to be confirmed by dedicated microscopic studies.

The obvious drawback of our eNSE model is linked to the treatment of nuclear clusters at subsaturation densities. A more realistic implementation would require consideration of species with much larger mass numbers, typically of the order of a thousand, as well as temperature and in-medium modifications of the binding energies of all species. Temperature effects on the energy density functional of clusters would modify both their binding energy and volume, and they would also suppress shell effects. In-medium modifications of the binding energies are essential for the dissolution of nuclei in dense environments, and they guarantee a natural transition to homogeneous matter. A practical way to account for these effects is via microscopic-inspired parametrizations, as done by Furusawa et al. (2013, 2017a,b), Typel (2018).

Acknowledgments

We acknowledge support from the Ministry of Education and Research, CNCS/CCCDI-UEFISCDI, Project No. PN-IV-P1-PCE-2023-0324 and partial support from Project No. PN 23 21 01 02. The two authors contributed equally to this work.

References

All Tables

All Figures

Fig. 1. Selected properties of the mixture of unbound nucleons and nuclei while approaching the transition to homogeneous matter for T = 0.23 MeV (left column) and 5.75 MeV (right column) and various proton fractions, as indicated in the legend. Considered are the relative difference between the baryonic free energy densities of inhomogeneous and homogeneous matter (bottom row), intrinsic volume of nuclei with A ≥ 2, Vcl, relative to the full volume, V, (middle row) and relative density of unbound nucleons (top row) as functions of density. For better readability, the quantities in the bottom panel are multiplied with arbitrary factors shown in the legend. The considered interaction is BBSk1.

In the text

Fig. 2. Mechanisms (top panel) and the density (bottom panel) for the transition to homogeneous matter as a function of temperature and proton number. The color legend in the top panel is as follows. Mechanism I (blue): ntr is taken as nB where Vcl(n) reaches the maximum for fixed values of T and Yp; mechanism II (green): ntr corresponds to the lowest value of nB where Vcl/V ≥ 0.9; mechanism III (cyan): ntr corresponds to the lowest value of nB where finhomo/fhomo − 1 ≥ −10−6. Cases where the transition is considered to occur at the highest density where convergence of Eqs. (43) and (44) is reached are marked in red. In the bottom panel the color levels correspond to ntr in units of fm−3. The considered interaction is BBSk1.

In the text

Fig. 3. Thermodynamic (in)stability of cold clusterized matter. Green curves represent the properties of clusterized stellar matter with fixed lepton chemical potentials (expressed in MeV) in terms of the baryon chemical potential, μB(=μn), as a function of baryonic number density, nB, (left panel) and total pressure, P, (right panel). Magenta curves represent the properties of clusterized baryonic matter with fixed proton chemical potentials (expressed in MeV) in terms of neutron chemical potential, μn, as a function of neutron number density, nn, (left panel) and baryonic pressure, PB, (right panel). The considered interaction is BBSk1. Here, μn and μp do not contain rest-mass contributions.

In the text

Fig. 4. Composition of an NS crust as a function of baryon number density in terms of numbers of neutrons (bottom panel) and protons (top panel) of nuclei. The eNSE results obtained in this work for β-equilibrated matter with T = 0.1 MeV are compared with the results of Pearson et al. (2018). The considered effective interactions are BSk22 and BSk25 (Goriely et al. 2013). We note that BSk25 results of Pearson et al. (2018) are scaled by a factor of 1/5. The vertical dashed lines respectively mark the transition to homogeneous matter for BSk22, ntr = 5.70 ⋅ 10−2 fm−3 and 6.35 ⋅ 10−2 fm−3, for the eNSE results and the results of Pearson et al. (2018). Horizontal dotted lines correspond to magic numbers: Z = 8, 20, 28 and N = 8, 20, 28, 50, 82.

In the text

Fig. 5. Mass fractions of various species present in the inner crust and their dependence on the effective interaction as functions of baryon number density. The eNSE results are for T = 0.1 MeV. Top panel: eNSE predictions confronted with the zero-temperature results of Pearson et al. (2018). Bottom panel: Mass fractions of 14He, 7H and unbound nucleons, as predicted by eNSE. The considered effective interactions are BSk22 and BSk25 from (Goriely et al. 2013) and BBSk1 and BBSk2 proposed here. As in Fig. 4, vertical dashed lines mark the transition to homogeneous matter.

In the text

	Fig. 6. Impurity parameter Qimp in the NS crust. Bottom panel: Predictions of BBSk1 at various temperatures. Top panel: Effective interaction effects at T = 0.1 MeV.
In the text

	Fig. 7. Average mass numbers of nuclei with A ≥ 20, the corresponding ⟨Z⟩/⟨A⟩ ratio, and the mass fraction they bind (X) as functions of baryonic particle density for Yp = 0.1 and different temperatures. The considered effective interaction is BBSk1.
In the text

	Fig. 8. Average mass and charge numbers of nuclei with Z ≥ 3 and A < 20 as well as the mass fraction they bind as functions of baryonic particle density for the cases considered in Fig. 7. We note that the X-range is not the same as in Fig. 7.
In the text

	Fig. 9. Same as in Fig. 7 but for T = 1 MeV and different values of Yp.
In the text

	Fig. 10. Same as in Fig. 8 but for the cases considered in Fig. 9. We note that the X-range is not the same as in Fig. 9.
In the text

	Fig. 11. Mass fraction of unbound nucleons (bottom panels) and their proton fraction (top panels) as functions of baryonic particle density for the cases considered in Fig. 7. We note that the X-range is not the same as in Fig. 7.
In the text

	Fig. 12. Same as in Fig. 11 but for the cases considered in Fig. 9. We note that the X-range is not the same as in Fig. 9.
In the text

	Fig. 13. Domains of nB − Yp where at T = 1.44 MeV the mass fractions of H isotopes (top panel) and 11He (red), 12He (blue), 13He (green), 14He (cyan) (bottom panel) exceed 0.2% and 1%, respectively. The color scale in the top panel corresponds to the mass fraction bound in each of the considered H isotopes. The effective interaction is BBSk1.
In the text

	Fig. 14. Neutron and proton chemical potentials as functions of baryonic particle density for Yp = 0.3 and 0.1 MeV ≤ T ≤ 5.75 MeV.
In the text

	Fig. 15. PB/(nBT) and energy per particle (E/A = eB/nB) as functions of baryonic particle density for the same conditions as in Fig. 14. For better readability, the contributions of electron and photon gases have been removed.
In the text

	Fig. 16. Same as in Fig. 14 but for T = 1 MeV and 0.01 ≤ Yp ≤ 0.6.
In the text

	Fig. 17. Same as in Fig. 15 but for the conditions considered in Fig. 16.
In the text