© The Authors 2024

1. Introduction

Extensive efforts have been directed toward the dynamical simulations of core-collapse supernovae explosions (Pons et al. 1999; Liebendörfer et al. 2005; Oertel et al. 2017; Burgio & Fantina 2018) and the subsequent formation and evolution of a proto-neutron star (PNS). A realistic equation of state (EOS) of asymmetric nuclear matter over a wide range of densities and temperatures is one of the most vital inputs to these calculations (Bethe 1990). Current insights suggest three distinct phases in this process (Pons et al. 1999; Barrère et al. 2022): (i) Roughly one second after the supernova core bounce, a comparatively cool central region, bordered by a hotter mantle, rapidly collapses, emitting neutrinos while accreting material. (ii) Over the next 20 s or so, a slowly developing state of the PNS can be identified. The system first deleptonizes and heats up the interior parts of the star, and then begins to cool down further through neutrino diffusion. (iii) After several minutes, the final state of the neutron star (NS) emerges, initially cooling by neutrino emission and later by photon emission from its surface (Pons & Viganò 2019).

A cold NS becomes stratified in three primary regions: the outer crust, the inner crust, and a uniform core. Each region is distinguished by its unique set of physical properties. In the NS core, matter forms a homogeneous liquid comprising primarily neutrons, together with some specific fractions of protons, electrons, and muons, ensuring the system remains in charge and β equilibrium. At even greater densities deeper inside the core, hyperons, other strange particles, and deconfined quarks can emerge (Shapiro & Teukolsky 1983; Haensel et al. 2007). Transitioning from the core toward the outer layers, the density decreases, causing positive charges to cluster individually, identified by charge Z, and arrange into a solid lattice. This organization minimizes Coulomb repulsion between them. The lattice is embedded in a gas of neutrons and a background of electrons, ensuring overall charge neutrality for the system. This region is known as the inner crust. In the deepest layers of the inner crust, nuclear structures may take on nonspherical forms, commonly referred to as “nuclear pasta,” driven by energy minimization (Baym et al. 1971a; Lorenz et al. 1993). Above the inner crust, in regions of lower densities, all neutrons become confined within nuclear clusters, giving rise to a lattice structure enriched with nuclei interspersed with a degenerate electron gas. This layer is defined as the outer crust (Baym et al. 1971b) and it extends from the interior, characterized by the neutron drip density, to the exterior, encompassing the envelope layer, which plays an important role in NS cooling (Potekhin et al. 2015; Dehman et al. 2023). The crust of the NS, although it comprises a small fraction of its total mass and radius, is decisive in various observable signals from NSs, such as, glitches in pulsar NSs (Piekarewicz et al. 2014), bursts and outbursts in intensely magnetized NSs (Beloborodov & Li 2016; Coti Zelati et al. 2018; Dehman et al. 2020), and NS asteroseismology phenomena (Steiner & Watts 2009; Sotani et al. 2012; Neill et al. 2023).

The first observation of a NS merger event, GW170817, by Abbott et al. (2017) has opened new avenues for studying the properties of matter under extreme conditions of densities and temperatures (Baiotti 2019). This event, resulting from the collision of two NSs, creates conditions that can lead to either the formation of a black hole or a notably massive NS. However, the final outcome of the GW170817 merger is still a subject of ongoing debate (Pooley et al. 2018). The pre-merger phase, known as the inspiral phase, leaves a unique imprint on the gravitational wave (GW) signal due to the tidal deformation of the stars involved, which is influenced by the EOS at zero or near zero temperature. This effect becomes noticeable early in the binary dynamics (Flanagan & Hinderer 2008; Hinderer et al. 2010; Shibata 2015; Krastev & Li 2019). The post-merger phase is expected to significantly heat the remnant, raising its temperature to tens of mega-electronvolt. This heating affects post-merger characteristics, such as the remnant’s lifespan, the GW spectrum, and the ejected mass, all of which are crucially dependent on the finite temperature EOS (Oechslin et al. 2007; Sekiguchi et al. 2011; Bauswein et al. 2013; Soma & Bandyopadhyay 2020).

A unified EOS for the crust and the core of NSs at zero temperature was introduced in Sharma et al. (2015), based on the BCPM nuclear energy density functional (Baldo et al. 2008, 2010, 2013; Baldo et al. 2017). In the BCPM functional, the bulk part of the energy density is entirely given by ab initio Brueckner–Hartree–Fock (BHF) calculations of nuclear matter. The bulk term is supplemented with surface and Coulomb terms, allowing one to describe the inhomogeneous nuclear structures of the NS crust and the homogeneous matter of the NS core simultaneously in a microscopic-based and consistent approach. Despite significant progress in developing the EOS for cold dense matter, the EOS at finite temperature has not been as thoroughly investigated. Our study aims to introduce a suitable method for examining thermal effects on the BCPM EOS. We particularly focus on the critical role of the hot inner crust.

Together with the growing interest in the high-energy astrophysical phenomena associated with compact stars, in recent years there has been much progress in the formulation of hot microscopic EOSs for the core of these systems (see, e.g., Lu et al. 2019; Logoteta et al. 2021; Figura et al. 2021, and references quoted therein). Few calculations of hot EOSs are available, however, that include the crust at finite temperature computed with the same model as the core, and these models generally are phenomenological ones (see Oertel et al. 2017; Raithel et al. 2021 for recent reviews of hot EOSs). Hence, in the present work we address a first exploration of a thermal EOS including the hot inner crust calculated on a microscopic basis by using the BCPM functional. We develop the study within a neutrino-free β-stable framework, assuming an isothermal profile of hot nuclear systems consisting of neutrons (n), protons (p), electrons (e), and muons (μ). In this context, the hot inner crust is present at relatively low densities and moderate temperatures. The mean free path of neutrinos in this regime is expected to be large enough for neutrinos not to be trapped (Reddy et al. 1999). For a broader scope, though, calculations with trapped neutrinos as well as configurations at constant entropy per baryon should also be addressed. The present thermal EOS is particularly relevant for compact objects in warm environments where neutrinos are not trapped and matter is β equilibrated. This includes scenarios such as the post-merger phase of binary NS events or the last stages of the evolution of PNSs (Oertel et al. 2017; Kumar & Bošnjak 2020; Bethe 1990), assuming that the temporal evolution is slow enough to justify that matter approaches the indicated conditions (Lu et al. 2019).

This paper is organized as follows. Section 2 describes the BCPM EOS. Sections 2.1 and 2.2 detail the EOS at zero and finite temperatures, respectively. The results obtained in this study are illustrated in Sect. 3, where we highlight the significant influence of the hot inner crust. The conclusions are drawn in Sect. 4.

2. BCPM and equation of state of β-stable matter

2.1. Zero-temperature equation of state

The microscopic BHF calculations can be directly employed to obtain the EOS of the liquid core of NSs, where the nuclei have dissolved into their constituent protons and neutrons. However, a BHF calculation of finite nuclei and nuclear structures in a NS crust is not yet feasible. To describe finite nuclei, retaining as much information from the ab initio BHF calculations as possible, the BCPM nuclear energy density functional was developed (Baldo et al. 2008, 2010, 2013, 2017). Initially designed for characterizing the ground state of finite nuclei, the BCPM functional comprises a bulk component given by the BHF results in symmetric nuclear matter and neutron matter through the local density approximation, supplemented with a finite-range term to account for the surface properties. The Coulomb, spin-orbit, and pairing contributions are also included (Baldo et al. 2008, 2013). Having in total four adjustable parameters (for the surface and spin-orbit terms only), BCPM describes the properties of finite nuclei with the same success as usual nuclear functionals that contain many more free parameters. The bulk component of the model was obtained using the BHF approach, where the calculations employed the Argonne v₁₈ nucleon-nucleon (NN) interaction, with the inclusion of three-body forces reduced to a two-body density-dependent term (Wiringa et al. 1995, chapter 1 of Baldo 1999, and references therein, Taranto et al. 2013). The resulting EOS for both symmetric and asymmetric nuclear matter satisfies several criteria set by heavy ion collisions and recent astrophysical observations.

In the BCPM functional, the bulk contribution to the energy per particle is the sum of the kinetic energy per particle (e_kin) and the potential energy per particle (v_int):

$$\begin{array}{c}\hfill e(n,\beta )=e_{\mathrm{kin}}(n,\beta )+v_{\mathrm{int}}(n,\beta ).\end{array}$$(1)

The kinetic energy per particle at zero temperature is that of a noninteracting cold Fermi gas with degeneracy factor of 2 and isospin asymmetry β:

$$\begin{array}{c}\hfill e_{\mathrm{kin}}(n,\beta )=\frac{1}{2}\frac{3}{5}\frac{{ħ}^2}{2m}(\frac{3{\pi }^{2}n}{2}{)}^{2/3}[(1+\beta {)}^{5/3}+(1-\beta {)}^{5/3}],\end{array}$$(2)

where n = n_n + n_p (n_n and n_p are, respectively, the neutron and proton number densities) and β = (n_n − n_p)/n is the isospin asymmetry parameter. The potential energy per particle in the bulk, v_int(n, β), is represented as a quadratic interpolation between the interaction energy per particle in symmetric nuclear matter (SNM), v_int(n, 0), and pure neutron matter (PNM), v_int(n, 1):

$$\begin{array}{c}\hfill v_{\mathrm{int}}(n,\beta )=v_{\mathrm{int}}(n,0)+(v_{\mathrm{int}}(n,1)-v_{\mathrm{int}}(n,0)){\beta }^2.\end{array}$$(3)

Here, v_int(n, 0) and v_int(n, 1) are given by the microscopic BHF calculations at various densities. For computational efficiency, in BCPM an accurate polynomial fit of the discrete BHF points is performed as a function of the density n (Baldo et al. 2010; Sharma et al. 2015). This fit remains valid up to a density n ≈ 0.4 fm⁻³. For bulk matter at higher density values, as the ones that can be found in NS cores, we utilize functional forms that provide excellent parametrizations of the BHF results for SNM and PNM at T = 0 (Burgio & Schulze 2010; Sharma et al. 2015).

The BCPM functional was employed in Sharma et al. (2015) to establish a unified EOS for the outer crust, inner crust and core of NSs at zero temperature. Many-body calculations of the inhomogeneous structures in the NS crust currently fall beyond the scope of the BHF approach, which is used for modeling the homogeneous core. The calculation of the NS crust with the BCPM functional maintains a consistent microscopic approach in describing the entire stellar structure. The NS crust is modeled within the Wigner–Seitz (WS) cell approximation, dividing space into noninteracting cells, each containing a single nuclear cluster in charge and β equilibrium. In the outer crust, matter is composed of fully ionized atomic nuclei, forming a solid lattice to minimize Coulomb repulsion, permeated by a degenerate electron gas. The critical data for constructing the outer crust EOS are the nuclear masses, sourced from the AME2012 evaluation (Audi et al. 2012), or, if they are unknown, calculated using the Hartree–Fock–Bogoliubov (HFB) method with the BCPM functional (Sharma et al. 2015). As the star’s average density increases, nuclei become increasingly neutron-rich, eventually leading to the inner crust, where neutrons begin to drip. The inner crust thus consists of nuclear clusters immersed in a gas of dripped neutrons and a background of electrons. In the bottom of the inner crust, surrounding the uniform core, a thin layer of pasta phases (nuclear clusters that adopt nonspherical shapes) is predicted by the BCPM calculations (Sharma et al. 2015).

A fully quantum calculation of the inner crust is a challenging task due to the presence of the neutron gas. To describe the cold inner crust in Sharma et al. (2015), self-consistent Thomas-Fermi (TF) calculations with BCPM were performed. This approach offers a significant advantage, because the EOS in the inner crust is primarily influenced by the neutron gas, which means that shell and pairing effects have a marginal impact on the EOS and the proton fraction (Pearson & Chamel 2022). Using this formalism, the EOS of the crust and the core are both determined through the BCPM energy density functional. The generalization of this formalism to finite temperature for the core and the inner crust is presented in Sect. 2.2.

2.2. Equation of state at finite temperature

In this work we extend the BCPM EOS to finite temperature. We introduce the thermal effects through the Fermi occupation numbers in the energy density functional, while keeping the interactions the same as at T = 0, as usually done in nuclear functionals. This approach is often known as the frozen-correlations approximation (Baldo & Ferreira 1999; Burgio et al. 2007), as it assumes the nucleon single-particle potential at finite T to be the same as at T = 0. It has been checked in hot BHF calculations that this assumption is verified with good accuracy, at least for temperatures up to about 30 MeV (Baldo & Ferreira 1999; Burgio et al. 2007; Burgio & Schulze 2010; Lu et al. 2019).

2.2.1. The liquid core

To determine the self-consistent occupation numbers in asymmetric nuclear matter at a given temperature, we minimize the thermodynamic potential, given by:

$$\begin{array}{cc}\hfill \mathrm{\Omega }& ={\displaystyle \sum _{q=n,p}\sum _k}{E}_q(n_n(k),n_p(k))-TS-{\displaystyle \sum _{q=n,p}\sum _k}{\mu }_q{n}_q(k),\hfill \\ \hfill & ={\displaystyle \sum _{q=n,p}\sum _k}{F}_q(n_n(k),n_p(k))-{\displaystyle \sum _{q=n,p}\sum _k}{\mu }_q{n}_q(k).\hfill \end{array}$$(4)

Here, E represents the internal energy, T is the temperature, F is the free energy, S is the entropy, and μ_q and n_q(k) are, respectively, the chemical potential and the occupation number of each type of nucleon, with q = n, p. The entropy per particle for asymmetric nuclear matter with proton and neutron fraction y_q = n_q/n at finite temperature is expressed as:

$$\begin{array}{c}\hfill s(n,\beta ,T)={\displaystyle \sum _{q=n,p}}{y}_q{s}_q(n,\beta ,T),\end{array}$$(5)

with s_q being the entropy per particle of each component:

$$\begin{array}{c}\hfill s_q(n,\beta ,T)=-{\displaystyle \sum _k}(n_q(k)ln{n}_q(k)+[1-n_q(k)]ln[1-n_q(k)]).\end{array}$$(6)

In Eqs. (4) and (6), n_q(k) represent the Fermi occupation numbers, which are the solutions of the variational equations obtained by applying the variational principle to the grand potential (Eq. (4)). These occupation numbers read as:

$$\begin{array}{c}\hfill n_q(k)=\frac{1}{1+e^{({\epsilon }_q(k)-{\mu }_q)/T}},\end{array}$$(7)

where ε_q(k) denotes the single-particle spectrum given by:

$$\begin{array}{c}\hfill {\epsilon }_q(k)=\frac{{ħ}^2{k}^2}{2m}+V_q(n,\beta ).\end{array}$$(8)

In the frozen-correlations approximation, the single-particle potential V_q depends on the baryon number density n and the isospin asymmetry β and is independent of temperature. It is given by:

$$\begin{array}{c}\hfill V_q(n,\beta )=\frac{\partial (n{v}_{\mathrm{int}}(n,\beta ))}{\partial n_q}=v_{\mathrm{int}}(n,\beta )+n\frac{\partial v_{\mathrm{int}}(n,\beta )}{\partial n_q},\end{array}$$(9)

with v_int(n, β) of the cold calculation, Eq. (3). The number density n_q of neutrons and protons reads as:

$$\begin{array}{c}\hfill n_q=\frac{2}{{(2\pi )}^3}{\displaystyle {\int }_0^{\infty }}{n}_q(k)d^{3}k=\frac{1}{2{\pi }^2}(\frac{2mT}{{ħ}^2}{)}^{3/2}{J}_{1/2}({\eta }_q),\end{array}$$(10)

where J_1/2(η_q) is the Fermi integral of index ν = 1/2:

$$\begin{array}{c}\hfill J_{1/2}({\eta }_q)={\displaystyle {\int }_0^{\infty }}\frac{z^{1/2}\mathrm{d}z}{1+e^{(z-{\eta }_q)}}·\end{array}$$(11)

Here, we have z = [ε_q(k)−V_q]/T and the fugacity reads η_q = [μ_q − V_q]/T.

At finite temperature, the total energy per particle of asymmetric nuclear matter becomes

$$\begin{array}{c}\hfill e(n,\beta ,T)=e_{\mathrm{kin}}(n,\beta ,T)+v_{\mathrm{int}}(n,\beta ),\end{array}$$(12)

where the interaction term, v_int(n, β), is defined in Eq. (3), and the nucleonic kinetic energy per particle reads as:

$$\begin{array}{c}\hfill \frac{1}{n}{\displaystyle \sum _{q=n,p}}\frac{{ħ}^2{\tau }_q}{2m},\end{array}$$(13)

where

$$\begin{array}{c}\hfill {\tau }_q=\frac{2}{{(2\pi )}^3}{\displaystyle {\int }_0^{\infty }}{k}^2{n}_q(k)d^{3}k=\frac{1}{2{\pi }^2}(\frac{2mT}{{ħ}^2}{)}^{5/2}{J}_{3/2}({\eta }_q),\end{array}$$(14)

and J_3/2(η_q) is the Fermi integral of index ν = 3/2:

$$\begin{array}{c}\hfill J_{3/2}({\eta }_q)={\displaystyle {\int }_0^{\infty }}\frac{z^{3/2}\mathrm{d}z}{1+e^{(z-{\eta }_q)}}·\end{array}$$(15)

The proton and neutron chemical potentials are expressed in terms of the temperature T, fugacity η_q (obtained by inverting Eq. (10) for a given nucleon density n_q), and the single-particle potential V_q:

$$\begin{array}{c}\hfill {\mu }_q(n,\beta ,T)={\eta }_{q}T+V_q(n,\beta ).\end{array}$$(16)

The pressure is given by:

$$\begin{array}{c}\hfill P(n,\beta ,T)=n[{\displaystyle \sum _{q=n,p}}{y}_q{\mu }_q(n,\beta ,T)-f(n,\beta ,T)],\end{array}$$(17)

where f(n, β, T) is the nucleonic free energy per particle.

In homogeneous neutrino-free stellar matter containing nucleons and leptons (electrons and muons), under the constraint of charge neutrality (n_p = n_e + n_μ), the conditions for β equilibrium are:

$$\begin{array}{c}\hfill {\mu }_n-{\mu }_p={\mu }_e,{\mu }_{\mu }={\mu }_e.\end{array}$$(18)

In our approach, nucleons are treated as nonrelativistic particles, while the leptons are considered as relativistic free particles. With the inclusion of leptons, the energy in Eq. (12) also contains the leptonic contribution:

$$\begin{array}{c}\hfill \frac{1}{n}\frac{2}{{(2\pi )}^3}{\displaystyle {\int }_0^{\infty }}{\epsilon }_l(k)n_l(k)d^{3}k,l=e,\mu ,\end{array}$$(19)

where ${\epsilon }_l(k)=\sqrt{{ħ}^2{k}^2{c}^2+m_l^2{c}^4}$ and n_l(k) is the occupation number of leptons:

$$\begin{array}{c}\hfill n_l(k)=\frac{1}{1+e^{({\epsilon }_l(k)-{\mu }_l)/T}}·\end{array}$$(20)

Likewise, the total entropy and the total pressure of the system are obtained by adding the contribution of free leptons to that of nucleons.

2.2.2. The hot inner crust

We have employed the TF method at finite temperature to calculate the EOS of the inhomogeneous matter of the crust using the BCPM functional. The nuclear clusters in the hot crust are computed inside spherical WS cells of radius R_c. Each cell is electrically neutral, and interactions between cells are neglected. At nuclear densities, the electrons exhibit highly relativistic behavior, as their Fermi momenta far exceed their rest mass. Therefore, they can be assumed to be uniformly distributed within the WS cell. Inside the cell, we impose β equilibrium, leading to the condition μ_n = μ_p + μ_e when neutrinos have left the star. Similarly to the core, the properties of the hot system within the WS cell are determined by minimizing the thermodynamic potential (Eq. (4)), which, in the crust, also includes nonuniform contributions (see below). It is important to note that at finite T, nuclei become unstable against nucleon evaporation, giving rise to a surrounding gas of evaporated nucleons. Thus, within the WS cell, a coexistence between the nuclear cluster plus gas and the gas alone occurs. Consequently, in the hot TF calculations in a WS cell there exist two solutions of the TF equations in equilibrium: one corresponding to the liquid-plus-gas (LG) phase and the other one corresponding to the gas (G) phase alone (Suraud 1987; Sil et al. 2002, and references therein).

In the context of the crust, the internal energy contributing to the thermodynamic potential is expressed as follows within the WS cell of volume V_c:

$$\begin{array}{c}\hfill E={\displaystyle {\int }_{V_c}}[\mathcal{H}(n_n,n_p)+{\mathcal{E}}_{\mathrm{el}}+{\mathcal{E}}_{\mathrm{coul}}+{\mathcal{E}}_{\mathrm{ex}}+m_p{n}_p+m_n{n}_n]\mathrm{d}\mathit{r}.\end{array}$$(21)

Here, ℋ(n_n, n_p) represents the nuclear energy density, where n_n = n_nr) and n_p = n_p(r) are the neutron and proton number densities, respectively, which are position dependent within the cell. In the hot TF approach, ℋ incorporates the kinetic energy density of protons and neutrons at finite temperature (see Eq. (13)). Additionally, it incorporates the cold interacting part 𝒱(n_n, n_p), as determined by the BCPM functional, covering both bulk (Eq. (3)) and surface contributions (Sharma et al. 2015). The surface term is taken as at zero temperature:

$$\begin{array}{cc}\hfill {\mathcal{E}}_{\mathtt{surf}}(n_n,n_p)& =\frac{1}{2}{\displaystyle \sum _{q,q^{\prime }}}{n}_q(\mathit{r}){\displaystyle \int v_{q{q}^{\prime }}}(\mathit{r}-{\mathit{r}}^{\mathbf{\prime }})n_{q^{\prime }}({\mathit{r}}^{\mathbf{\prime }})\mathrm{d}{\mathit{r}}^{\mathbf{\prime }}\hfill \\ \hfill & -\frac{1}{2}{\displaystyle \sum _{q,q^{\prime }}}{n}_q(\mathit{r})n_{q^{\prime }}(\mathit{r}){\displaystyle \int v_{q{q}^{\prime }}}({\mathit{r}}^{\mathbf{\prime }})\mathrm{d}{\mathit{r}}^{\mathbf{\prime }}.\hfill \end{array}$$(22)

The second term in Eq. (22) is subtracted to avoid contamination of the bulk part derived from the microscopic nuclear matter calculations. For the finite-range form factors, we use the same Gaussian shape v_qq′(r) as in Baldo et al. (2010), Sharma et al. (2015). The term ℰ_el in Eq. (21) represents the energy density arising from the motion of electrons at finite temperature, with their energy per particle given in Eq. (19). For the densities and temperatures of our calculations of the crust, muons did not occur in this region of the star.

At finite temperature, the direct Coulomb contribution in a WS cell corresponding to the LG and G phases, as discussed in Sil et al. (2002), is expressed as:

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{LG}}^{\mathrm{coul}}(n_{\mathrm{LG}}^p,n_e)=\frac{1}{2}(n_{\mathrm{LG}}^p(\mathit{r})-n_e)(V_{\mathrm{LG}}^{\mathrm{coul},p}(\mathit{r})-V_e(\mathit{r})),\end{array}$$(23)

and

$$\begin{array}{cc}\hfill {\mathcal{E}}_{\mathrm{G}}^{\mathrm{coul}}(n_{\mathrm{G}}^p,n_e)& =\frac{1}{2}(n_{\mathrm{G}}^p(\mathit{r})-n_e)(V_{\mathrm{G}}^{\mathrm{coul},p}(\mathit{r})-V_e(\mathit{r}))\hfill \\ \hfill & +n_L^p(\mathit{r})(V_{\mathrm{G}}^{\mathrm{coul},p}(\mathit{r})-V_e(\mathit{r})),\hfill \end{array}$$(24)

where $n_L^p=n_{\text{LG}}^p-n_{\text{G}}^p$ is the proton density of the nuclear cluster. The potentials $V_{\mathrm{LG}(\mathrm{G})}^{\mathrm{coul},p}(\mathit{r})$ and V_e(r) are given by:

$$\begin{array}{c}\hfill V_{\mathrm{LG}(\mathrm{G})}^{\mathrm{coul},p}(\mathit{r})={\displaystyle \int \frac{e^2{n}_{\mathrm{LG}(\mathrm{G})}^p({\mathit{r}}^{\mathbf{\prime }})}{|\mathit{r}-{\mathit{r}}^{\mathbf{\prime }}|}}\mathrm{d}{\mathit{r}}^{\mathbf{\prime }},V_e(\mathit{r})={\displaystyle \int \frac{e^2{n}_e}{|\mathit{r}-{\mathit{r}}^{\mathbf{\prime }}|}}\mathrm{d}{\mathit{r}}^{\mathbf{\prime }}.\end{array}$$(25)

The direct part of the single-particle Coulomb potential, obtained by performing the functional derivatives of the LG and G direct Coulomb energies with respect to n_LG and n_G, respectively, is the same for both phases and reads:

$$\begin{array}{c}\hfill V_{\mathrm{coul}}^d(\mathit{r})={\displaystyle \int \frac{e^2(n_{\mathrm{LG}}^p({\mathit{r}}^{\mathbf{\prime }})-n_e)}{|\mathit{r}-{\mathit{r}}^{\mathbf{\prime }}|}}\mathrm{d}{\mathit{r}}^{\mathbf{\prime }}.\end{array}$$(26)

The total Coulomb energy in the LG and G phases is comprised of the direct contributions, as described in Eqs. (23) and (24), respectively, along with the exchange contribution from protons and electrons, which is calculated at the Slater level:

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{LG}(\mathrm{G})}^{\mathrm{ex}}(n_{\mathrm{LG}(\mathrm{G})}^p,n_e)=-\frac{3}{4}(\frac{3}{\pi }{)}^{1/3}{e}^2(n_{\mathrm{LG}(\mathrm{G})}^p{(\mathit{r})}^{4/3}+n_e^{4/3}).\end{array}$$(27)

Taking functional derivatives of the thermodynamical potential (Eq. (4)) with the energy given by Eq. (21), which includes Coulomb effects in both LG and G phases, with respect to $n_{\text{LG}}^q$ and $n_{\text{G}}^q$ for neutrons and protons, we obtain the following set of coupled equations, as described in Sil et al. (2002):

$$\begin{array}{cc}& T{\eta }_{\mathrm{LG}}^q(\mathit{r})+V_{\mathrm{LG}}^q(\mathit{r})+V_{\mathrm{LG}}^{\mathrm{coul}}(\mathit{r})={\mu }_q,\hfill \end{array}$$(28)

$$\begin{array}{cc}& T{\eta }_{\mathrm{G}}^q(\mathit{r})+V_{\mathrm{G}}^q(\mathit{r})+V_{\mathrm{G}}^{\mathrm{coul}}(\mathit{r})={\mu }_q.\hfill \end{array}$$(29)

Here, $V_{\text{LG}(\text{G})}^q$ represents the nuclear part of the single-particle potential in the LG(G) phases, and $V_{\mathrm{LG}(\mathrm{G})}^{\mathrm{coul}}$ is the total Coulomb potential, which is the sum of the direct term from Eq. (26) and the exchange part obtained from derivatives with respect to $n_{\text{LG}}^p$ ($n_{\text{G}}^p$) of the exchange energy as given in Eq. (27).

For a certain temperature T and average baryon density n in the WS cell, assuming a radius R_c (equivalently, a baryon number A) for the cell, the set of variational equations is solved self-consistently using the method described in Sil et al. (2002). This method allows for the determination of the composition (A, Z) of minimal free energy per baryon in β equilibrium. Next, we search for the optimal cell size (i.e., optimal baryon number) for the given baryon density n by repeating the calculation for various values of R_c. In the present study at finite T, we restrict ourselves to spherical WS cells. This choice is motivated by the fact that in the BCPM calculations at T = 0, nuclear pasta shapes appear just in a small density region near the transition to the core, cf. Fig. 4 of Sharma et al. (2015). Furthermore, as the relevant energy surfaces are very flat near the optimal configurations and they present tiny energy differences between spherical and nonspherical cells, the different geometrical shapes have a negligible impact on the EOS (pressure vs density) (Sharma et al. 2015; Pearson et al. 2020). Therefore, for simplicity, here we consider spherical cells only. However, the consideration of pasta phases would be important in applications where the detailed structure of the crust is required, as pasta may strongly influence elasticity, transport and other properties of the NS crust (Haensel et al. 2007).

We note that at finite T, the distinctive nuclear shell effect is eroded, as thermal excitation disrupts the well-defined energy levels within the nucleus present at T = 0. It is a well-known consequence of the smearing of the Fermi surface caused by the Fermi occupation numbers. In hot nuclear systems with temperatures higher than T ≈ 2 − 3 MeV the shell effects have vanished (Barranco & Buchler 1981; Brack et al. 1985; Pi et al. 1986), meaning that these systems are optimal objects for TF calculations like ours.

We are focused on calculating the EOS in various regions of the NS, with a specific emphasis on determining the structure and EOS in the crust using the BCPM model at a finite T. To obtain the EOS in the inner crust, it is necessary to compute the pressure, which is determined by taking appropriate derivatives of the energy with respect to the size of the WS cell. As shown in Appendix A of Sharma et al. (2015) (also see Pearson et al. 2018), at T = 0 the pressure in the inner crust is a result of the neutron and electron gases within which the nuclear structures are embedded. At finite T, it consists of the contributions of the nucleons (neutrons and protons) and the free electrons in the gas phase, plus the Coulomb exchange pressure of the charged particles in the gas.

The exploration of the presence of hyperons and other exotic degrees of freedom in hot supranuclear matter is extremely interesting and a major topic of ongoing research (Oertel et al. 2017; Marques et al. 2017). We recall, however, that since we are largely focused on the regime of densities of the crust, where hyperons are not expected to occur in any significant amount (Menezes & Providência 2017), we consider the picture where the baryonic composition of the star corresponds to purely nucleonic matter. Recent Bayesian analyses of the possible behaviors of the EOS of dense matter indicate that the existing astrophysical measurements of NSs are compatible with the fully nucleonic hypothesis for the composition of dense matter (Thi et al. 2021).

3. Impact of hot crust

To study the EOS in the NS crust at finite temperature, we start by examining the transition between crust matter and uniform matter across different temperature ranges. Although it is possible to estimate the crust-core transition density from the core side by looking at the threshold for the instability of uniform matter against clustering (Kubis 2004; Xu et al. 2009; Moustakidis 2012; Gonzalez-Boquera et al. 2019), in this study, we seek the transition density from the perspective of the crust to emphasize the behavior of the crust at finite temperature. In the inner crust, we search for the nuclear composition that provides the minimal free energy per baryon in β equilibrium using the hot TF method discussed in the previous section. This calculation of the crust includes self-consistently the Coulomb and surface effects, which are absent in uniform matter. As in the T = 0 case (Sharma et al. 2015), the transition from the crust to the uniform matter is determined by an energy criterion. This transition occurs when the free energy per baryon of the homogeneous phase is lower than that obtained for the clustered phase within the WS cell. At T = 0 the BCPM spherical-cell calculations predict that the crust-core interface is located at baryon density n = 0.08 fm⁻³.

In Fig. 1, we show the regions where the clustered and the uniform matter are the most stable phases as a function of the temperature and the average baryon density, derived from hot TF calculations with BCPM. This figure features a thick blue line in the density-temperature plane that separates these two regions. The figure also illustrates that the crust-uniform matter transition boundary, as a function of temperature, is a bivalued function that decreases with increasing temperature until it converges to a single endpoint. This endpoint occurs at a specific limiting temperature (T_lim), predicted to be of 7.21 MeV, beyond which crust-uniform matter transitions are not feasible. The upper and lower transition densities illustrated in Fig. 1 are collected in Table 1 for specific temperature values. The concept of the upper transition density is analogous to the transition density between the inner crust and the core of NSs at zero temperature, whereas the lower transition density does not have a T = 0 counterpart.

Fig. 1. Upper and lower branches of the transition density between the uniform matter and the clustered matter computed at different temperature values. The blue shaded area corresponds to the presence of nuclear clusters, determined using TF calculations within a WS cell.

Figure 2 displays the difference in free energy per baryon between crust matter and uniform matter for T = 0, 2, 5, and 7.5 MeV. At each temperature, the crust phase is favored while the result is negative. Starting from T = 0, the gap between the free energy of the crust phase and the uniform phase is seen to shrink rapidly as T increases. We note that in this figure we have divided the T = 0 result by 10 and the T = 2 result by 4 to display the four temperatures with the same vertical scale. It can also be observed that the T = 0 curve crosses the zero axis of Fig. 2 only once, at n = 0.08 fm⁻³, which is the crust-core transition density of the cold case. However, with raising T, the free energy gap between crust and uniform matter presents an inverted-bell shape, such that for T < T_lim, as illustrated by the T = 2 and 5 MeV cases, the gap vanishes at two different density points, corresponding to the lower ($n_t^l$) and the upper ($n_t^u$) transition density for the given temperature. The separation between $n_t^l$ and $n_t^u$ closes with higher T, until $n_t^l=n_t^u$ for the limiting temperature T_lim = 7.21 MeV. Above T_lim, the free energy of the clustered matter remains higher than for the uniform matter, as shown in the plot by the T = 7.5 MeV case.

Fig. 2. Difference between the calculated free energy per particle of the clustered matter and of the uniform matter for T = 0, 2, 5, and 7.5 MeV. Notice that the T = 0 result is divided by a factor 10 and the T = 2 MeV result by a factor 4.

Some differences are noticed between the zero and the finite temperature scenarios. At zero T, the inner crust of a NS consists of nuclear clusters that encapsulate all the protons present in the WS cell, amid a gas of dripped neutrons. As the average baryon density rises, the contribution from the neutron gas also increases. A threshold density exists, above which uniform matter becomes the most stable phase. With thermal excitation, more neutrons evaporate from the clusters into the surrounding gas, resulting in an increased neutron number in the gas phase. Protons also evaporate at finite T, contributing to the formation of a charged gas. In this scenario, a lower transition density appears, below which uniform matter becomes again the most stable phase. Our calculations have shown that the reemergence of the uniform matter in the lower transition density branch is due to both the proton Coulomb contribution and the interaction between protons and neutrons in the gas phase. As depicted in Fig. 1 and seen also in Fig. 2, the density range at a given T in which clustered matter is the most stable phase diminishes with increasing T. The upper and lower boundaries defining this region converge, ultimately merging into a single point at the limiting temperature. This scenario is akin to the phenomena observed in nuclear matter when analyzing instabilities within a homogeneous medium (Barranco & Buchler 1981; Lattimer & Swesty 1991; Hempel & Schaffner-Bielich 2010).

It is noteworthy that at the transition density between the crust and uniform matter for both the lower and upper branches, our findings indicate that the proton fraction at these points is approximately 0.03. Additionally, we have computed the plasma parameter, defined as Γ = (Z_cle)²/(R_cT), at the transition density points illustrated in Table 1 and Fig. 1. We find that Γ is consistently less than 175 at all transition densities, indicating that the clusters are in the liquid phase (Pi et al. 1986). This result is in agreement with previous findings reported in the literature (Aguilera et al. 2008).

Certain studies of the NS EOS at finite temperature assume a cold crust because the EOS of the hot crust for the considered nuclear interaction is not available. Consequently, these EOSs are composed of a hot core starting from a specified value of the crust-core transition density computed at zero T, supplemented by a cold crust contribution obtained from existing literature (see, e.g., Bombaci 1996; Lu et al. 2019, and references therein). In the left panel of Fig. 3, we illustrate this type of NS EOSs using the BCPM density functional at various temperatures while considering a cold crust. In this case, we have maintained the zero-temperature density transition point between the inner crust and the core at n (T = 0) = 0.08 fm⁻³, as reported in Sharma et al. (2015). For temperatures T ≲ 2 MeV, the discontinuity between the cold crust and the hot core EOS is acceptable. However, as the temperature increases, for example T = 5, 10, and 20 MeV, the discontinuity becomes more and more pronounced.

Fig. 3. EOSs of β-stable neutrino-free matter at different temperature values of T = 0, 2, 5, 10, and 20 MeV. In the left panel, we use a cold crust assuming the transition between the cold crust and the uniform core fixed at n = 0.08 fm−3. In the right panel, we use a hot inner crust and the crust–uniform matter transition is determined as described in Sect. 3, consistently with Fig. 1.

When a constant T is considered (isothermal profile), the radius of the star cannot be defined by the usual condition of vanishing pressure at the surface of the star. This is because at finite T, the pressure does not vanish even at very low densities (Buchler & Coon 1977), and, consequently, the thermal effect causes isolated NSs to expand. To address this, a temperature drop from the hot interior to the surface of the star is necessary. A proper determination of the temperature profile of the star would require dynamical simulations of thermal transport coupled with the T-dependent EOS. Given the difficulty and existing uncertainties, other alternatives have been explored in the literature, such as using a neutrino sphere to ensure the temperature drops to zero at low density (Gondek et al. 1997). For this reason, in our subsequent calculations for modeling NSs, we coupled the finite-temperature EOS of the core and the crust at n = 10⁻⁴ fm⁻³ with the cold outer crust to obtain a cold surface. The rationale for selecting the value of 10⁻⁴ fm⁻³ aligns with the outer-to-inner crust transition density of the BCPM functional at T = 0 (Sharma et al. 2015). Moreover, this value is compatible with the range of densities predicted in the literature where the neutrino sphere is employed (Gondek et al. 1997; Strobel et al. 1999; Fischer et al. 2009). Therefore, our EOS at finite T for NS calculations is built up as follows. For temperatures below the limiting temperature of 7.21 MeV, we compute the EOS taking into account the hot inner crust and the uniform matter according to the different density regions shown in Figs. 1, 2 and Table 1. When the temperature exceeds 7.21 MeV, the existence of the hot inner crust ceases; consequently, we employ calculations for hot uniform matter throughout. For all T values, we rely on the cold outer crust for densities below 10⁻⁴ fm⁻³. It is worth mentioning that we have compared our NS mass-radius (M–R) relation results for T = 15 MeV (see below) with those from Gondek et al. (1997), which were obtained by considering an isothermal temperature profile at T = 15 MeV and a neutrino sphere, instead of a cold outer crust. We found that our results are consistent with theirs.

Plots of our EOSs at different temperatures with the hot inner crust prescription are depicted in the right panel of Fig. 3. The impact of the hot temperature EOS in the crustal region can be appreciated at baryon densities below 0.08 fm⁻³ by comparing the left and right panels of this figure. In Fig. 4, we plot a wider range of the EOSs at finite temperatures, extending up to T = 50 MeV. A consistent trend of growth is observed for all temperatures.

Fig. 4. EOSs of β-stable neutrino-free matter at different temperature values ranging from T = 0 to 50 MeV. For T lower than Tlim = 7.21 MeV, we display the NS EOS containing the hot inner crust using the BCPM functional. At T > 7.21 MeV, the NS is composed of homogeneous dense matter, thus we illustrate the results of the bulk part of the hot BCPM functional.

The M–R relations of NSs composed of neutrino-free β-stable matter are illustrated in Fig. 5 at different temperature values. The colorbar reflects the central density of the star. The M–R relations have been calculated by solving the Tolman–Oppenheimer–Volkoff equations (Shapiro & Teukolsky 1983). The observational masses used to depict the constraint of M ∼ 2 M_sun in the M–R graph are 1.97 ± 0.04 M_sun from the PSR J1614-2230 pulsar (Demorest et al. 2010) and 2.01 ± 0.04 M_sun from the PSR J0348+0432 pulsar (Antoniadis et al. 2013). The T = 0 EOS satisfies this limit for the maximum mass, as well as the 1.4 M_sun radius constraint of R_1.4 = 11.9 ± 1.4 km extracted by the LIGO-Virgo collaboration from the analysis of GW170817 (Abbott et al. 2018), and is also compatible with the mass and radius constraints of 2.08 ± 0.07 M_sun measured in the PSR J0740+6620 pulsar (Fonseca et al. 2021) and R_1.4 = 12.45 ± 0.65 km deduced from NICER data on PSR J0740+6620 (Miller et al. 2021). We see in Fig. 5 that the maximum mass of a NS slightly increases as T grows, compared to the cold T = 0 case, but remains almost constant with temperature, in consonance with previous findings (Bombaci 1996; Burgio & Schulze 2010; Lu et al. 2019). As T increases, there is, however, a strong effect on the stellar radius, which for the same mass shifts toward a larger value when the star is hot. As a consequence the compactness of the NS, 𝒞 = GM/(Rc²), decreases with higher temperatures. For instance, if we consider the maximum mass configurations, we have 𝒞 = 0.30 for T = 0, whereas 𝒞 = 0.23 for T = 20 MeV. This temperature dependence of the compactness may have observable consequences in the simulations of core-collapse supernovae, PNSs, and hot remnants of binary NS mergers (Oechslin et al. 2007). On the other hand, as evident from the colorbar, for a given T the central density takes the highest value for the maximum mass and decreases as the NS mass decreases and the NS radius increases. Furthermore, for a given value of the mass, specially for not too small masses, the central density of the star is rather independent of temperature, showing a mild reduction with higher T. For example, for M = 1.5 M_sun (2 M_sun) the central density experiences a slight decrease of 5% (2.5%) from T = 0 to 20 MeV.

Fig. 5. Computed mass–radius relations of NSs with an isothermal profile at temperatures of T = 0, 2, 5, 15, and 20 MeV. The observational mass constraints are 1.97 ± 0.04 Msun (gray band) from the PSR J1614-2230 pulsar (Demorest et al. 2010) and 2.01 ± 0.04 Msun (cyan band) from the PSR J0348+0432 pulsar (Antoniadis et al. 2013).

To emphasize the impact of the hot inner crust, we present a comparison of the M–R relations in Fig. 6. On the one hand, we consider a full cold crust below a density of 0.08 fm⁻³ (dotted lines), and on the other hand, a hot inner crust (solid lines, the same of Fig. 5) at T = 5 and 15 MeV. We include the M–R relation at T = 0 as a reference. The two arrows indicate the shift in the M–R relation for T = 5 and T = 15 MeV when considering the hot inner crust. The EOSs with a cold crust predict M–R relations that are similar to those at T = 0, particularly for low temperature values, for instance T = 5 MeV, and around the maximum mass, that is, M ∼ 2 M_sun. However, with the hot inner crust, the influence of temperature on the M–R relation becomes much more pronounced. Clearly, when one attaches a cold crust below the cold crust-core transition density (0.08 fm⁻³) to the hot core, the effect of temperature on the radius of the star is largely underestimated. That is, even if a hot core EOS is implemented, with a cold crust the NS radius increases very moderately compared with the T = 0 EOS. When comparing the central densities predicted by considering a full cold crust and a hot inner crust (refer to Fig. 6), we notice that, for a given NS radius, the central density shifts to larger values when the hot crust is considered.

Fig. 6. Mass–radius relations of NSs with an isothermal profile at temperatures of T = 0, 5, and 15 MeV. The dotted lines correspond to EOSs with cold crust and the solid lines correspond to EOSs with hot crust (the same as in Fig. 5). Note that the horizontal scale goes up to 20 km, whereas in Fig. 5 it goes up to 30 km. The observational constrains are the same as in Fig. 5. The two arrows indicate the shift in the M–R relation for T = 5 and T = 15 MeV.

4. Conclusion

In Sharma et al. (2015), the BCPM energy density functional was used to derive a unified EOS for cold NSs based on ab initio BHF calculations. Astrophysical studies of core-collapse supernovae, formation of PNSs, or the remnants of binary NS mergers, furthermore demand modeling hot compact objects that reach high temperatures. Here, we predicted the thermal properties of stellar matter by generalizing at finite temperature the study with the BCPM functional. To ascertain the EOS of the hot crust across various temperature values, we conducted hot self-consistent TF calculations using the BCPM functional involving spherical droplet configurations within WS cells. In this first exploration of a unified hot EOS with the BCPM functional, we were largely driven by investigating the temperature effects on the baryonic matter and developed our study for isothermal β equilibrated, neutrino-transparent configurations.

A prominent feature of our findings is that in hot stars there is a strong reduction of the value of the crust-core transition density. The microscopic calculations predict that the inner crustal region completely dissolves at temperatures exceeding 7.21 MeV, a point we referred to as the limiting temperature. Moreover, our analysis of the hot inner crust has unveiled two distinct branches of the interface between uniform and clustered matter below the limiting temperature. The branch with higher density values represents the upper transition density, marking the transition between the inner crust and the core; similar to the crust-core transition at zero temperature. Conversely, the branch with lower density values indicates a lower transition density, which does not occur at zero temperature, suggesting the reappearance of uniform matter in the outer layers of the hot NS.

Furthermore, we conducted a study on the M–R relation at finite temperature. Here, the hot EOS obtained in our calculations was complemented by the cold BCPM EOS for the outer crust (Sharma et al. 2015), which extends up to a density of 10⁻⁴ fm⁻³. This step was taken to ensure vanishing pressure at the surface of the star. A more accurate, and desirable, exploration needs a temperature profile of the star, but the latter is still uncertain; hence, we left the formulation of the calculations under a NS temperature profile for a future study. Due to the existence of the limiting temperature, above T_lim = 7.21 MeV the star consists entirely of uniform matter at densities greater than those corresponding to the border of the cold outer crust (10⁻⁴ fm⁻³). Our findings highlight the growing relevance of thermal effects on the M–R relation. On the one hand, for a given NS mass, the radius increases with temperature, and this effect is more pronounced at higher temperatures and at low stellar masses. On the other hand, for a given NS radius, the mass also increases with temperature but saturates to a value close to the maximum mass of a cold NS for sufficiently high temperatures. When comparing our results with those obtained with a cold crust attached to the hot core (see Fig. 6), we observed that the thermal effects from the hot core on the M–R relation are moderate, while the influence of the hot crust on the M–R relation is much more substantial.

To summarize, in this study, we successfully extended the BCPM EOS, for neutrino-free β-stable matter, to finite temperatures, thereby ensuring a meticulous treatment of the inner crust. The BCPM energy density functional offers a unified microscopic EOS for NSs under both cold and hot conditions. Furthermore, we determined M–R relations and central densities for hot compact objects at various temperature levels in an isothermal picture. Our results suggest that for NSs in hot environments, such as the end products of NSs merger events, it is essential to incorporate the hot inner crust.

Acknowledgments

This investigation was initiated in collaboration with Prof. Artur Polls (†12 August 2020), who made significant contributions to this study. We would like to dedicate this work in honor of his memory. We acknowledge useful discussions with J. N. De and J. A. Pons. C.D. is supported by the Prometeo Excellence Programme grant CIPROM/2022/13 from the Generalitat Valenciana, ERC Consolidator Grant “MAGNESIA” No. 817661, and has partial support from NORDITA and grant SGR2021-01269 while working on this article. This work is also partially supported by the María de Maeztu CEX2020-001058-M Excellence Unit Programme. M.C. and X.V. acknowledge partial support from Grants No. PID2020-118758GB-I00 and No. CEX2019–000918-M (through the “Unit of Excellence María de Maeztu 2020-2023” award to ICCUB) from the Spanish MCIN/AEI/10.13039/501100011033.

References

All Tables

All Figures

	Fig. 1. Upper and lower branches of the transition density between the uniform matter and the clustered matter computed at different temperature values. The blue shaded area corresponds to the presence of nuclear clusters, determined using TF calculations within a WS cell.
In the text

	Fig. 2. Difference between the calculated free energy per particle of the clustered matter and of the uniform matter for T = 0, 2, 5, and 7.5 MeV. Notice that the T = 0 result is divided by a factor 10 and the T = 2 MeV result by a factor 4.
In the text

	Fig. 3. EOSs of β-stable neutrino-free matter at different temperature values of T = 0, 2, 5, 10, and 20 MeV. In the left panel, we use a cold crust assuming the transition between the cold crust and the uniform core fixed at n = 0.08 fm−3. In the right panel, we use a hot inner crust and the crust–uniform matter transition is determined as described in Sect. 3, consistently with Fig. 1.
In the text

	Fig. 4. EOSs of β-stable neutrino-free matter at different temperature values ranging from T = 0 to 50 MeV. For T lower than Tlim = 7.21 MeV, we display the NS EOS containing the hot inner crust using the BCPM functional. At T > 7.21 MeV, the NS is composed of homogeneous dense matter, thus we illustrate the results of the bulk part of the hot BCPM functional.
In the text

	Fig. 5. Computed mass–radius relations of NSs with an isothermal profile at temperatures of T = 0, 2, 5, 15, and 20 MeV. The observational mass constraints are 1.97 ± 0.04 Msun (gray band) from the PSR J1614-2230 pulsar (Demorest et al. 2010) and 2.01 ± 0.04 Msun (cyan band) from the PSR J0348+0432 pulsar (Antoniadis et al. 2013).
In the text

Fig. 6. Mass–radius relations of NSs with an isothermal profile at temperatures of T = 0, 5, and 15 MeV. The dotted lines correspond to EOSs with cold crust and the solid lines correspond to EOSs with hot crust (the same as in Fig. 5). Note that the horizontal scale goes up to 20 km, whereas in Fig. 5 it goes up to 30 km. The observational constrains are the same as in Fig. 5. The two arrows indicate the shift in the M–R relation for T = 5 and T = 15 MeV.

In the text