Open-source nuclear equation of state framework based on the liquid-drop model with Skyrme interaction

The equation of state (EOS) of dense matter is an essential ingredient for numerical simulations of core-collapse supernovae and neutron star mergers. The properties of matter near and above nuclear saturation density are uncertain, which translates into uncertainties in astrophysical simulations and their multimessenger signatures. Therefore, a wide range of EOSs spanning the allowed range of nuclear interactions are necessary for determining the sensitivity of these astrophysical phenomena and their signatures to variations in input microphysics. We present a new set of finite temperature EOSs based on experimentally allowed Skyrme forces. We employ a liquid-drop model of nuclei to capture the nonuniform phase of nuclear matter at subsaturation density, which is blended into a nuclear statistical equilibrium EOS at lower densities. We also provide a new, open-source code for calculating EOSs for arbitrary Skyrme parametrizations. We then study the effects of different Skyrme parametrizations on thermodynamical properties of dense astrophysical matter, the neutron star mass-radius relationship, and the core collapse of 15 and 40 solar mass stars.

Figure 1

Density dependence of the symmetry energy $\mathcal{S}\left(n\right)$ for all considered Skyrme parametrizations. The thick curves show $\mathcal{S}\left(n\right)$ for uniform nuclear matter (neutrons and protons only) obtained from Eq. (49) with ${\epsilon }_B(n,y)=E_B(n,y)/n$ given by Eq. (7). The thin curves correspond to $\mathcal{S}\left(n\right)$ for the full high-density EOS at zero temperature, allowing for nonuniform and uniform nuclear matter. At densities below the transition to uniform matter, $\mathcal{S}\left(n\right)$ is obtained from Eq. (49) with ${\epsilon }_B(n,y)$ replaced by $F/n$ with $F$ from Eq. (2). Note that the high-density ($n\gg n_0$) behavior of $\mathcal{S}\left(n\right)$ is highly uncertain. The high-density shape of $\mathcal{S}\left(n\right)$ for the Skyrme parametrizations shown here is a mere artifact of the expansion about $n_0$ and is not necessarily physical. At low density, the binding energy of nearly symmetric nuclei increases the value of the symmetry energy for nonuniform matter.

Figure 2

Comparison with L&S. We show the ratio of nuclear pressure for a range of proton fractions obtained with our $\mathrm{LS}{220}^{†}$ implementation and with the original L&S implementation $\mathrm{LS}{220}^{*}$. The solid black curve delineates where the heavy-nuclei number fraction $X_i=u{n}_i/n$ changes from zero to a nonzero value. Below and to the left of the curve, matter is nonuniform, while above and to the right of the line it is uniform. The dashed line shows where the nuclear pressure is zero. Differences between $\mathrm{LS}{220}^{†}$ and $\mathrm{LS}{220}^{*}$ are largest near this line. The wide horizontal band at the bottom of the panels marks the region where the original L&S implementation does not converge for nonuniform nuclear matter and assumes that the system is uniform.

Figure 3

Total pressure per nucleon $P/n$ in $\mathrm{MeV}{\mathrm{baryon}}^{-1}$ for the NRAPR Skyrme parametrization [59] and proton fractions $y=0.01$, 0.10, 0.30, and 0.50. The solid black curves denote, from top to bottom, the adiabats at entropies $s={10}^n{k}_B{\mathrm{baryon}}^{-1}$ for $n=2$ to $n=-3$ in $-0.5$ increments. The dashed black curves correspond, from right to left, to the isoergs for specific energies $\epsilon =100$, 30, 10, 3, 0, and $-3\text{MeV}{\mathrm{baryon}}^{-1}$. Note that only the $y=0.30$ and $y=0.50$ panels contain the $\epsilon =-3\text{MeV}{\mathrm{baryon}}^{-1}$ isoerg. The pressure per nucleon is dominated by electrons, positrons, and photons in large portions of the density-temperature space. Only at the highest densities, at and above saturation density, is the pressure dominated by the nucleon contributions and the impact of strong interactions.

Figure 4

Temperature along adiabats with specific entropy $s=10$, 1.0, 0.1, and $0.01k_B{\mathrm{baryon}}^{-1}$ in the single-nucleus approximation and for all considered Skyrme parametrizations. Note that electrons, positrons, and photons are included. The adiabats differ mostly in regions dominated by nucleonic pressure around and above saturation density $n_0\sim 0.16{\text{fm}}^{-3}$ and for $y\simeq 0.30$ and $T\lesssim 1$ MeV, where the number of free neutrons varies significantly between parametrizations.

Figure 5

Adiabatic index $\mathrm{\Gamma }$ along the $s=1k_B{\mathrm{baryon}}^{-1}$ adiabat for full EOSs in the single-nucleus approximations. At low densities and proton fractions $y\gtrsim 0.1$, electrons dominate and $\mathrm{\Gamma }\sim 4/3$. At high densities and $y=0.50$, $\mathrm{\Gamma }$ is roughly the same for all EOSs reflecting the well constrained properties of symmetric nuclear matter. The sharp discontinuity is due to the transition between nonuniform and uniform nuclear matter at $n\simeq 0.1{\text{fm}}^{-3}$. It becomes smoother at lower proton fraction due to large free neutron contributions. Also, as the proton fraction decreases, differences between parametrizations increase due largely to variations in the density-dependent symmetry energy (cf. Fig. 1).

Figure 6

Average heavy nucleus mass number $\overline{A}$ along the $s=1k_B{\mathrm{baryon}}^{-1}$ adiabat as a function of density for the considered Skyrme parametrizations in the single-nucleus approximation (SNA) and, at low densities, for nuclear-statistical equilibrium (NSE) with 23, 807, and 3 335 nuclides (we discuss our NSE treatment and matching to SNA in Sec. 7a). For the SNA curves, differences in nuclear sizes result from differences in symmetry energy and the surface properties obtained for each parametrization. The 3335 nuclide NSE network exhibits large oscillations in $\overline{A}$. These are due to nuclear shell effects included implicitly in the nuclear masses.

Figure 7

Neutron-proton chemical potential difference $\widehat{\mu }={\mu }_n-{\mu }_p$ along the $s=1.0k_B{\mathrm{baryon}}^{-1}$ adiabat as a function of density for the considered Skyrme parametrizations at proton fractions $y=0.01$, 10, 0.30, and 0.50. Note that we multiply $\widehat{\mu }$ by a factor of $-25$ in the bottom right panel showing the $y=0.5$ case. $\widehat{\mu }$ is sensitive to the density dependence of the symmetry energy and the differences between parametrizations seen here correlate with those in Fig. 1.

Figure 8

Pressure $P$ (top panel), specific internal energy $\epsilon$ (plus an additive constant ${\epsilon }_0=2\times {10}^{19}\text{erg}{\text{g}}^{-1}$; center panel), and proton fraction $y$ for low temperature neutrinoless beta equilibrated matter (bottom panel). We show results for five select Skyrme parametrizations that span the range of maximum neutron star masses shown in Fig. 9. The $\mathrm{LS}{220}^{*}$ (LS220) curve uses L&S's (our) implementation of the L&S $K_0=220$ MeV parametrization. Differences between LS220 and $\mathrm{LS}{220}^{*}$ are due to the L&S implementation limit of proton fractions $y\ge 0.035$ and a small difference in the proton masses used (cf. Sec. 4a). The proton fraction $y$ at high densities, $\rho \gtrsim {10}^{14.5}\text{g}{\text{cm}}^{-3}$, mirrors the high-density behavior of the symmetry energy $\mathcal{S}\left(n\right)$ (cf. Fig. 1).

Figure 9

Mass-radius curves for cold, beta-equilibrated neutron stars (NSs) obtained by solving the Tolman-Oppenheimer-Volkoff equations for the considered Skyrme parametrizations. We summarize NS properties in Table 4. The gray strip represents the mass of the NS $\mathrm{PSR}\mathrm{J}0348+0432$, $M_{\mathrm{J}0348+0432}=2.01\pm 0.04M_{\odot }$ [79]. The yellow region indicates the NS mass-radius constraints from model A of Nättilä et al. [12]. Besides the LS220 parametrization, only SLy4 and (barely) KDE0v1 and NRAPR satisfy the $M_{\max }\gtrsim 2M_{\odot }$ constraint. Differences between our implementation of LS220 and the original L&S implementation ($\mathrm{LS}{220}^{*}$) are due to the lower limit of $y\ge 0.035$ in the latter (cf. Fig. 8).

Figure 10

Radial rest-mass density (top panel) and proton fraction profiles (bottom panel) of cold, beta-equilibrated $1.4M_{\odot }$ neutron stars (NSs) obtained with the considered Skyrme parametrizations. Note that LS220 is an outlier, yielding the lowest central density, the largest radius, and the highest central proton fraction. This is due primarily to its large $L$ parameter and the linear behavior of its density-dependent symmetry energy, which results in the smallest symmetry energies below saturation density and the largest symmetry energies above saturation density of all the EOS considered here (cf. Fig. 1). We summarize key NS quantities in Table 4.

Figure 11

Radial rest-mass density (top panel) and proton fraction (bottom panel) profiles of cold, beta-equilibrated maximum-mass neutron stars (NSs) as predicted by the considered Skyrme parametrizations. Note that the maximum mass varies between parametrizations. Table 4 summarizes key NS quantities for all parametrizations. As in the $1.4M_{\odot }$ NS case shown in Fig. 10, the LS220 parametrization is an outlier and yields the lowest central density, the highest central proton fraction, and the largest radius for its maximum-mass NS configuration.

Figure 12

Mass-radius curves for cold, beta-equilibrated neutron stars (NSs) obtained with the SkT1 parametrizations and its various high-density modifications. The gray strip represents the mass $\mathrm{PSR}\mathrm{J}0348+0432M_{\mathrm{J}0348+0432}=2.01\pm 0.04M_{\odot }$ [79] and the yellow region indicates the NS mass-radius constraints from model A of Nättilä et al. [12].

Figure 13

Relative error in the specific entropy as a function of density during adiabatic compression tests with the full EOS based on the SLy4 parametrization and the single-nucleus approximation at all densities. The different panels show results for proton fractions $y=0.01,0.10,0.3$, and 0.5. All curves start at $T=0.01$ MeV and we choose the initial densities to obtain starting values of the specific entropy of $s_0=0.1$, 0.2, 0.5, and $1.0k_B{\mathrm{baryon}}^{-1}$. Thick (thin) curves correspond to tests with the high (standard) resolution tables (cf. Table 6). In most cases the thick lines (obtained from high resolution tables) are very close to 1 at all densities and lie on top of each other, which may make their visualization difficult. For $s_0\gtrsim 0.5{k_B\text{baryon}}^{-1}$ and proton fractions $y\gtrsim 0.3$, the standard-resolution tables perform very well. In stellar collapse, the specific entropy always stays higher than $0.5{k_B\text{baryon}}^{-1}$ and proton fractions below $\sim 0.3$ are not reached until the final phase of collapse. Errors at lower $s_0$ and $y$ are largely numerical and are reduced by employing the high resolution table. However, large changes in entropy can still occur near the first-order phase transition between nonuniform and uniform matter at $n\simeq 0.1{\mathrm{fm}}^{-3}$.

Figure 14

Results from core collapse simulations with the $15M_{\odot }$ progenitor and EOSs generated with various Skyrme parametrizations and a 3335-nuclide NSE EOS at low densities. From top to bottom, we show $1.2\mathrm{s}$ of postbounce evolution of the central density ${\rho }_c$, central temperature $T_c$, and central specific entrop $s_c$. Note that the entropy stays roughly constant (modulo mild numerical oscillations) throughout the postbounce evolution, as it should for thermodynamically consistent EOSs. As postbounce accretion adds mass to the proto-NS, it contracts, which is marked by an increase in ${\rho }_c$ and softer EOSs result in a steeper increase. The splitting of the $T_c$ evolutions into two groups of parametrizations can be understood by considering that those resulting in lower temperatures have a larger effective nucleon mass (see Sec. 2a and Table 2).

Figure 15

Results of core collapse simulations with the $15M_{\odot }$ progenitor and the LS220 Skyrme parametrization. We compare results obtained for pure SNA (at all densities) with results from simulations that use SNA at high densities, smoothly matched to an NSE EOS with varying number of nuclides at low densities (see Sec. 7a for details). From top to bottom, we plot the postbounce evolution of shock radius $R_s$, accretion rate ${\dot{M}}_{300}$ at a radius of 300 km, average nuclear mass number 5 km above the shock, and difference in specific internal energy, $\mathrm{\Delta }\epsilon$ 5 km above and 5 km below the shock. Note that the pure-SNA simulation predicts a slightly larger shock radius than the SNA+NSE simulations between $\sim 50$ and $\sim 80$ ms after bounce. This is a consequence of the SNA simulation having a slightly lower accretion rate and slightly less bound nuclei crossing the shock in that interval.

Figure 16

From top to bottom we plot the postbounce time evolution of the central density ${\rho }_c$, central temperature $T_c$, and central specific entropy $s_c$ for the black hole (BH) formation simulations with the $40M_{\odot }$ progenitor. Proto-NS collapse and BH formation are marked by a sudden extreme steepening of the ${\rho }_c$ slope. The different graphs correspond to simulations with different Skyrme parametrizations. We employ a 3335-nuclide NSE EOS at low densities (cf. Sec. 7a). Note that the specific entropy stays, as it should, roughly constant (modulo numerical noise that can be reduced with higher-resolution EOS tables) throughout the postbounce evolution and up to BH formation. Thermal pressure support in the proto-NS mantle plays an important role in supported proto-NS masses that are $0.2\text{–}0.6M_{\odot }$ higher than the maximum cold NS mass. Thermal contributions are largest for those parametrizations that result in low effective nucleon masses and higher proto-NS temperatures.