A new equation of state for core-collapse supernovae based on realistic nuclear forces and including a full nuclear ensemble

The nuclear equation of state (EOS) is an important ingredient to determine the dynamics of core-collapse supernovae [1–4]. It provides information on thermodynamical quantities such as pressure and composition of high-density matter realized in the supernova core, which is either a mixture of nuclei or uniformly distributed nucleons, for given baryon number density n_B, temperature T and proton fraction Y_p. The pressure is obviously a quantity that has a direct influence on the dynamics of supernova, determining, e.g., the strength of the shock wave at core bounce and the evolution of proto-neutron stars (PNS) in the post-bounce phase [5–9]. The abundance of various nuclei is important, affecting the neutronization of matter through electron captures on nuclei and the opacity for neutrino via the coherent scattering on heavy nuclei [10–12]. Note that the population of light nuclei such as deuteron, triton, helion and alpha particles may have some implications for shock revival [13–16].

One of the critical elements in the construction of the nuclear EOS is the theoretical approach to nuclear interactions, on which the bulk properties of individual nuclei as well as the free energy of homogeneous nuclear matter depend strongly. In the previous EOSs used in supernova simulations, Skyrme type interactions [17] or meson-exchange models were adopted [18–26]. Both of them are based not on realistic interactions in vacuum that are known experimentally but on phenomenological medium-dependent interactions, whose model parameters are determined so that some nuclear properties such as the saturation density are reproduced. Recently, Togashi et al [27] have constructed an EOS based on a realistic two-body nuclear potential, Argonne v18 [28], employing the variational method [29]. This is the first EOS of such kind, meant for supernova simulations.

The calculation of the EOS for homogeneous nuclear matter is not sufficient if the EOS were to be used in supernova simulations. We need in addition to treat inhomogeneous nuclear matter composed of various nuclei and dripped nucleons, which exists generically at sub-nuclear densities. The common choice is the so-called single nucleus approximation (SNA), in which the ensemble of heavy nuclei is represented just by a single, supposedly the most abundant nucleus. It was adopted in the two standard EOSs, widely employed in the simulations of core-collapse supernovae [17–20] and also in the latest EOS [27]. Either the compressible liquid drop model (LDM) or the Thomas–Fermi approximation was employed in those EOSs. In such approximations, in-medium effects on the representative nucleus, such as compression, can be taken into account easily. Note that the abundance of various nuclei is important for the estimation of neutralization rates [30] and that it should be determined consistently with the modifications of nuclear properties in-medium, since nuclear masses depends on the population of other nuclei whereas the nuclear abundance is dictated by the nuclear masses. It was demonstrated that the average mass number and the total mass fraction of heavy nuclei do not coincide with those of the representative nucleus in general [31, 32].

In the last few years, some research groups have constructed the multi-nuclear species EOSs [21, 33–35], in which the full ensemble of nuclei are taken into account with various finite-density and temperature effects, such as the formation of nuclear pastas [36–40], smearing of shell effects in heavy nuclei [41–44] and self- and Pauli-energy shifts in light nuclei [45, 46], being only partially considered. We have also constructed such an EOS (FYSS EOS) [24–26], in which the mass formula was extended for heavy nuclei and took a quantum approach for light clusters to accommodate all the in-medium effects mentioned above. Although the FYSS EOS may be better in the treatment of in-medium effects, it is sometimes pointed out that the relativistic mean field (RMF) with the TM1 parameter set [47], on which FYSS EOS is based, gives too high a symmetry energy  36.9 MeV compared with the conventional values of  29.0–33.7 MeV [48–50].

The purpose of the present study is to construct a new, multi-nuclear species EOS in the framework of the FYSS EOS but with a different EOS for uniform nuclear matter as a base. We use the EOS by Togashi et al [29], who applied the variational method, one of the most popular techniques for the construction of the EOS of cold neutron star matter, using realistic nuclear forces expressed with the Argonne v18 potential [28]. Note that in their EOS the symmetry energy is  30.0 MeV, and with the inclusion of three-body forces, the saturation properties of uniform nuclear matter are satisfactorily reproduced. This is the first multi-nucleus EOS based on realistic nuclear forces. This article is organized as follows. In section 2, we describe the formulation adopted in this paper to construct the EOS. The results are shown in section 3 in comparison with the FYSS EOS derived from the RMF with the TM1 parameter set. We also compare the applications of these two EOSs to the one-dimensional supernova simulation in section 4. The paper is wrapped up with a summary and some discussion in section 5.

The calculation below is based on the formulation of the FYSS EOS to describe the nuclear ensemble [24–26] and the variational principle method for uniform nuclear matter [27, 29]. For details, we refer the reader to Furusawa et al [26] and Togashi et al [29]. The free energy of the FYSS EOS is represented as

$$\begin{eqnarray}&&f={f}_{p,n}+\displaystyle \sum _{j}{n}_{j}({E}_{j}^{t}+{M}_{j})+\displaystyle \sum _{i}{n}_{i}({E}_{i}^{t}+{M}_{i}),\end{eqnarray} \tag{ 1 }$$

where ${f}_{p,n}$ is the free energy density of the free nucleons outside nuclei, ${n}_{j/i}$ are the number densities of individual light nuclei j with the proton number ${Z}_{j}\leqslant 5$ and heavy nuclei i with the proton number $6\leqslant {Z}_{i}\leqslant 1000$, ${E}_{i/j}^{t}$ are the translational energies of heavy and light nuclei and ${M}_{i/j}$ are masses of heavy and light nuclei. In this work, the variational calculation takes the place of the RMF calculation in the previous model of the FYSS EOS, which provides the free energy density of nucleons ${f}_{p,n}$, the bulk energies of heavy nuclei in M_i and the self-energy shifts of light nuclei in M_j.

2.1. Free energy of free nucleons

The variational calculations for the free energy of free nucleons are based on [29, 51, 52]. We start from the nuclear Hamiltonian composed of a two-body potential V_ij and three-body potentials V_ijk, as in the Fermi hypernetted chain variational calculations [53]:

$$\begin{eqnarray}&&H=-\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}+\displaystyle \sum _{i\lt j}^{N}{V}_{{ij}}+\displaystyle \sum _{i\lt j\lt k}^{N}{V}_{{ijk}},\end{eqnarray} \tag{ 2 }$$

where m is set to be the mass of the neutron.

The free energy derived from two-body interactions is obtained by using an extension of the variational method by Schmidt and Pandharipande [54, 55] with AV18 two-body potential [28] and the healing distance condition, which reproduces the internal energy per baryon of symmetric nuclear matter and neutron matter at zero temperature obtained by Akmal et al [53]. The internal energy for three-body interactions is based on the UIX three-body potential [56, 57] with a phenomenological correction that reproduces empirical saturation properties. The entropy is expressed with average occupation probabilities of single quasi-nucleon states as in the case of a non-interacting Fermi gas [54]. The total free energy per baryon is minimized with respect to effective masses for nucleons. The optimized free energies agree reasonably with those of Akmal et al [53] at zero temperature and those of Mukherjee [58] at finite temperature.

We account for the excluded-volume effect. The local number densities of free protons and neutrons are defined as ${{n}}_{{p}/{n}}^{\prime }=({{N}}_{{p}/{n}})/(V-{V}_{{N}})$ where V is the total volume, V_N is the volume of all nuclei, and ${N}_{{p}/{n}}$ are the numbers of free protons and neutrons. The free energy density of free nucleons is given by

$$\begin{eqnarray}&&{f}_{{\rm{p}},{n}}=\displaystyle \frac{(V-{V}_{{N}})}{V}({{n}}_{{p}}^{\prime }+{{n}}_{{n}}^{\prime })\omega (T,{{n}}_{{p}}^{\prime }+{{n}}_{{n}}^{\prime },{{n}}_{{p}}^{\prime }/({{n}}_{{p}}^{\prime }+{{n}}_{{n}}^{\prime })),\end{eqnarray} \tag{ 3 }$$

where $\omega (T,{n}_{{\rm{B}}},x)$ is the free energy density per baryon in the unoccupied volume, $V-{V}_{{N}}$, for nucleons obtained from the variational calculation at temperature T, the local number density of baryon ${{n}}_{{\rm{B}}}={{n}}_{{p}}^{\prime }+{{n}}_{{n}}^{\prime }$ and the charge fraction $x={{n}}_{{p}}^{\prime }/({{n}}_{{p}}^{\prime }+{{n}}_{{n}}^{\prime })$.

2.2. Masses of heavy nuclei

The masses of heavy nuclei are assumed to be the sum of the bulk, Coulomb, surface and shell energies: ${M}_{i}={E}_{i}^{{\rm{B}}}+{E}_{i}^{{\rm{C}}}+{E}_{i}^{{\rm{Sf}}}+{E}_{i}^{{\rm{Sh}}}$. We define the saturation density of nuclei n_si(T) as the baryon number density, at which the free energy per baryon $\omega (T,{n}_{{\rm{B}}},{Z}_{i}/{A}_{i})$ given by the variational method becomes a minimum, which is the same as that for the free energy density of free nucleons. Thus, ${n}_{{si}}(T)$ depends on the temperature T and the charge fraction in each nucleus ${Z}_{i}/{A}_{i}$. At high temperatures $T\geqslant {T}_{{ci}}$, the free energy, $\omega (T,{n}_{{\rm{B}}},{Z}_{i}/{A}_{i})$, has no minimum because of finite entropy. Then the saturation density ${n}_{{si}}(T)$ above T_ci is assumed to be equal to ${n}_{{si}}({T}_{{ci}})$. When the saturation density n_si so obtained is lower than the baryon number density of the whole system n_B, we reset the saturation density as the baryon number density ${n}_{{si}}={n}_{{\rm{B}}}$. This treatment of the saturation density is needed to obtain reasonable bulk energies at high temperatures and densities [25].

The bulk, surface, Coulomb and shell terms in M_i are expressed as

$$\begin{eqnarray}&&{E}_{i}^{{\rm{B}}}={A}_{i}\{\omega (T,{n}_{{si}},{Z}_{i}/{A}_{i})\},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}{E}_{i}^{{\rm{C}}}=\left\{\begin{array}{ll}\displaystyle \frac{3}{5}{\left(\displaystyle \frac{3}{4\pi }\right)}^{-1/3}{e}^{2}{n}_{{si}}^{2}{\left(\displaystyle \frac{{Z}_{i}-{n}_{{\rm{p}}}^{\prime }{V}_{i}^{{N}}}{{A}_{i}}\right)}^{2}{{V}_{i}^{{N}}}^{5/3}D({u}_{i}) & ({u}_{i}\leqslant 0.3),\\ \displaystyle \frac{3}{5}{\left(\displaystyle \frac{3}{4\pi }\right)}^{-1/3}{e}^{2}{n}_{{si}}^{2}{\left(\displaystyle \frac{{Z}_{i}-{n}_{{\rm{p}}}^{\prime }{V}_{i}^{{N}}}{{A}_{i}}\right)}^{2}{{V}_{i}^{{\rm{B}}}}^{5/3}D(1-{u}_{i}) & ({u}_{i}\geqslant 0.7),\end{array}\right.\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}{E}_{i}^{{\rm{Sf}}}=\left\{\begin{array}{ll}4\pi {r}_{{N}{i}}^{2}\,{\sigma }_{i}{\left(1-\displaystyle \frac{{{n}}_{{p}}^{\prime }+{{n}}_{{n}}^{\prime }}{{{n}}_{{si}}}\right)}^{2}{\left(\displaystyle \frac{{T}_{{cs}}^{2}-{T}^{2}}{{T}_{{cs}}^{2}+{T}^{2}}\right)}^{5/4} & ({u}_{i}\leqslant 0.3),\\ 4\pi {r}_{{\rm{B}}{i}}^{2}\,{\sigma }_{i}{\left(1-\displaystyle \frac{{{n}}_{{p}}^{\prime }+{{n}}_{{n}}^{\prime }}{{{n}}_{{si}}}\right)}^{2}{\left(\displaystyle \frac{{T}_{{cs}}^{2}-{T}^{2}}{{T}_{{cs}}^{2}+{T}^{2}}\right)}^{5/4} & ({u}_{i}\geqslant 0.7),\end{array}\right.\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}{E}_{i}^{{\rm{Sh}}}=\left\{\begin{array}{ll}{E}_{i0}^{{\rm{Sh}}}\displaystyle \frac{{\tau }_{i}}{{\rm{\sinh }}{\tau }_{i}} & (\rho \leqslant {10}^{12}\,{\rm{g}}\,{{\rm{cm}}}^{-3}),\\ {E}_{i0}^{{\rm{Sh}}}\left(\displaystyle \frac{{\tau }_{i}}{{\rm{\sinh }}{\tau }_{i}}\right)\left(\displaystyle \frac{{\rho }_{0}-\rho }{{\rho }_{0}-{10}^{12}\,{\rm{g}}\,{{\rm{cm}}}^{-3}}\right) & (\rho \gt {10}^{12}\,{\rm{g}}\,{{\rm{cm}}}^{-3}).\end{array}\right.\end{eqnarray} \tag{ 7 }$$

The bulk energy E_i^B is evaluated from the free energy per baryon of uniform nuclear matter $\omega (T,{n}_{{\rm{B}}},x)$ given by the variational method at ${n}_{{\rm{B}}}={n}_{{si}}$ and $x={Z}_{i}/{A}_{i}$ for a given temperature T.

The Coulomb energy E_i^C is obtained by integrating the Coulomb potential over the Wigner–Seitz cell containing nucleus i, dripped free protons and uniformly distributed electrons [24]. The cell is set to satisfy charge neutrality within the volume as ${V}_{i}=({Z}_{i}-{n}_{{p}}^{\prime }{V}_{{N}i})/({n}_{{e}}-{n}_{{p}}^{\prime })$, where V_Ni is the volume of the nucleus in the cell and can be calculated as ${V}_{{N}{i}}={A}_{{i}}/{n}_{{si}}$ and ${n}_{{e}}={Y}_{{p}}{n}_{{\rm{B}}}$ is the number density of electrons. The vapor volume and nucleus volume fraction in the cell are given by ${V}_{i}^{{\rm{B}}}={V}_{i}-{V}_{i}^{{N}}$ and ${u}_{i}={V}_{i}^{{N}}/{V}_{i}$, respectively. Then E_i^C is expressed as equation (5) where $D({u}_{i})=1-\tfrac{3}{2}{u}_{i}^{1/3}+\tfrac{1}{2}{u}_{i}$ and e is the elementary charge. Each nucleus is assumed to enter the nuclear pasta phase individually when the volume fraction reaches u_i = 0.3 and that the bubble shape is realized when it exceeds 0.7.

The surface energy E_i^Sf is obtained from surface area times surface tension ${\sigma }_{i}={\sigma }_{0}-{A}_{i}^{2/3}/{(4\pi {r}_{{N}{i}}^{2})[{S}_{s}(1-2({Z}_{i}/{A}_{i})}^{2})]$, where the values of the constants, ${\sigma }_{0}=1.15\,\mathrm{MeV}$ fm⁻³ and ${S}_{s}=45.8\,\mathrm{MeV},$ are adopted from Lattimer and Swesty [17]. We may have to choose the values that are consistent with the EOSs for uniform nuclear matter. The simpler estimate, however, is adopted in this work, since there is no concrete way to derive the surface coefficients only from the EOS for uniform nuclear matter. The last two factors in equation (6) describe the reduction of the surface energy when the density contrast between the nucleus and the nucleon vapor decreases and when the temperature grows. The critical temperature is set to be ${T}_{{cs}}=18\,\mathrm{MeV}$ [59] and the radii of nuclei and bubbles are expressed as ${r}_{{N}{i}}={(3/4\pi {V}_{i}^{{N}})}^{1/3}$ and ${r}_{{\rm{B}}{i}}={(3/4\pi {V}_{i}^{{\rm{B}}})}^{1/3}$. We use cubic polynomials of u_i for the interpolation of Coulomb and surface energies between the droplet and bubble phases to ensure continuous and smooth connections at u_i = 0.3 and u_i = 0.7. We referred Watanabe et al [36] for the criterion of intermediate states, although the thermodynamic quantities are hardly affected by its choice. The interpolation is mainly required to ensure the smooth change in mass fractions of nuclei [25].

The shell energy at zero density and temperature is assumed to be the difference between the experimental or theoretical mass data [60, 61] and our liquid drop mass formula in the vacuum limit without the shell term: ${E}_{i0}^{{\rm{Sh}}}={M}_{i}^{{\rm{data}}}-{[{E}_{i}^{{\rm{B}}}+{E}_{i}^{{\rm{C}}}+{E}_{i}^{{\rm{Sf}}}]}_{T,{{n}}_{{p}/{n}}^{\prime },{{n}}_{{e}}=0}$. The factor ${\tau }_{i}/\sinh {\tau }_{i}$ in equation (7) expresses the washout of shell effects approximately with ${\tau }_{i}=2{\pi }^{2}T/(41{A}_{i}^{-1/3})$, which reproduces qualitatively the feature that shell effects disappear around $T\,\sim \,\mathrm{2.03.0}$ MeV [26]. The linear interpolation, $({\rho }_{0}-\rho )/({\rho }_{0}-{10}^{12}\,{\rm{g}}\,{\mathrm{cm}}^{-3})$, accounts for the disappearance of shell effects at high densities where ρ and ${\rho }_{0}$ are the density and the saturation density, respectively [24].

2.3. Masses of light nuclei

We describe light nuclei as quasi-particles outside heavy nuclei, whose masses are assumed to be given by ${M}_{j}={M}_{j}^{{\rm{data}}}+{\rm{\Delta }}{E}_{j}^{{\rm{Pa}}}+{\rm{\Delta }}{E}_{j}^{{\rm{SE}}}+{\rm{\Delta }}{E}_{j}^{{\rm{C}}}$. The Pauli-energy shift ${\rm{\Delta }}{E}_{j}^{{\rm{Pa}}}$ for deuteron (d), triton (t), helion (h) and alpha particle (α), which is fitted to the result of quantum statistical calculations [46], is given by [45]

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{j}^{{\rm{Pa}}}({{n}}_{{p}{l}},{{n}}_{{n}{l}},T)=-{\tilde{{n}}}_{j}\left[1+\displaystyle \frac{{\tilde{{n}}}_{j}}{2{\tilde{{n}}}_{j}^{0}(T)}\right]\delta {B}_{j}(T),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\delta {B}_{j}(T)=\left\{\begin{array}{ll}{a}_{j,1}/{T}^{3/2}[1/\sqrt{{y}_{j}}-\sqrt{\pi }{a}_{j,3}\exp ({a}_{j,3}^{2}{y}_{j})\mathrm{erfc}({a}_{j,3}\sqrt{{y}_{j}})] & {\rm{for}}\ j=d,\\ {a}_{j,1}/{({{Ty}}_{j})}^{3/2}\ \ & {\rm{for}}\ j=t,h,\alpha ,\end{array}\right.\end{eqnarray} \tag{ 9 }$$

where ${\tilde{n}}_{j}=2({Z}_{j}\ {n}_{{pl}}+{{N}}_{j}\ {n}_{{nl}})/{A}_{j}$, ${y}_{j}=1+{a}_{j,2}/T$, ${n}_{{pl}/{nl}}$ are the local proton and neutron number densities that include light nuclei as well as free nucleons and ${\tilde{n}}_{j}^{0}(T)={B}_{j}^{0}/\delta {B}_{j}(T)$ with the binding energy in vacuum, B_j⁰. The Pauli-energy shifts for the light nuclei other than $d,t,h$ and α are calculated in the same way as that for α.

The self-energy shift ${\rm{\Delta }}{E}_{j}^{{\rm{SE}}}$ is the sum of the self-energy shifts of individual nucleons in light nuclei ${\rm{\Delta }}{E}_{{n}/{p}}^{{\rm{SE}}}={{\rm{\Sigma }}}_{{n}/{p}}^{0}(T,{{n}}_{{p}}^{\prime },{{n}}_{{n}}^{\prime })-{{\rm{\Sigma }}}_{{n}/{p}}(T,{{n}}_{{p}}^{\prime },{{n}}_{{n}}^{\prime })$ with ${{\rm{\Sigma }}}^{0}$ and Σ being the vector and scalar potentials of nucleon and the contribution from their effective masses ${\rm{\Delta }}{E}_{j}^{\mathrm{eff}.\mathrm{mass}}=(1-{m}^{* }/m){s}_{j}$ with ${m}^{* }={m}_{{\rm{B}}}-{{\rm{\Sigma }}}_{{n}/{p}}(T,{{n}}_{{p}}^{\prime },{{n}}_{{n}}^{\prime })$:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{j}^{{\rm{SE}}}({n}_{{p}}^{\prime },{n}_{{n}}^{\prime },T)=({A}_{j}-{Z}_{j}){\rm{\Delta }}{E}_{{n}}^{{\rm{SE}}}+{Z}_{j}{\rm{\Delta }}{E}_{{p}}^{{\rm{SE}}}+{\rm{\Delta }}{E}_{j}^{\mathrm{eff}.\mathrm{mass}}.\end{eqnarray} \tag{ 10 }$$

For the potentials ${{\rm{\Sigma }}}^{0}$ and Σ, we employ the parametric formula for RMF with the DD2 parameter set, equations (A1) and (A2) in [45]. In the previous EOS [26], the potentials are evaluated consistently with the RMF with the TM1 parameter set employed for free nucleons. The self-energy shifts for the light nuclei other than $d,t,h$ and α are set to be zero. The parameters ${a}_{j/1}$, ${a}_{j/2}$, ${a}_{j/3}$, and s_j are given in table 1 in [45]. The Coulomb energy shift is obtained from the same calculation for Coulomb energy of heavy nuclei, equation (5), by subtracting its vacuum limits although the nuclear pasta phases are not considered here. This mass evaluation for light clusters is stitched together by using a variety of different model assumptions, whereas mass fractions of light clusters in the FYSS EOS agree well with some experimental predictions [50, 62].

2.4. Translational energies of nuclei and thermodynamical quantities

The translational energy of nuclei is based on that for an ideal Boltzmann gas and is given by

$$\begin{eqnarray}&&{F}_{i/j}^{t}=T\left\{\mathrm{log}\left(\displaystyle \frac{{n}_{i/j}}{{g}_{i/j}(T){({M}_{i/j}T/2\pi {{\hslash }}^{2})}^{3/2}}\right)-1\right\}\left(1-\displaystyle \frac{{n}_{B}}{{n}_{s}}\right).\end{eqnarray} \tag{ 11 }$$

The last factor takes account of the excluded-volume effect in the same way as in Lattimer and Swesty [17]; each nucleus can move in the free space that is not occupied by baryons and its translational motion is suppressed at high densities. The internal degree of freedom ${g}_{i/j}(T)$ is defined as ${g}_{i}(T)=({g}_{i}^{0}-1){\tau }_{i}/\sinh {\tau }_{i}+1$ for heavy nuclei and ${g}_{j}(T)={g}_{j}^{0}$ for light nuclei where the spin degree of freedom of the ground state is defined by ${g}_{i/j}^{0}$. The factor ${\tau }_{i}/\sinh {\tau }_{i}$ is incorporated washout of shell effects as in equation (7).

The abundances of nuclei as a function of n_B, T and Y_p are obtained by minimizing the model free energy in equation (1) with respect to the number densities of nuclei and nucleons under baryon and charge conservations. The nuclear abundance of each nucleus is determined by the nucleon chemical potentials ${\mu }_{p/n}$ and the nuclear mass ${M}_{i/j}$ under chemical equilibrium. Since both ${\mu }_{p/n}$ and ${M}_{i/j}$ depend on the local proton and neutron densities ${n}_{p/n}^{\prime }$ in our model, we solve the equations relating ${\mu }_{p/n}$ and ${n}_{p/n}^{\prime }$ as well as the two conservation equations to determine the four variables: ${\mu }_{p}$, ${\mu }_{n}$, ${n}_{p}^{\prime }$ and ${n}_{n}^{\prime }$.

Thermodynamical quantities other than the free energy are derived by partial differentiations of the optimized free energy density, in which all terms for in-medium effects are properly taken into account to ensure thermodynamical consistency. For example, the baryonic pressure and the entropy are obtained by the first derivative with respect to the baryonic density and the temperature, respectively.

For comparison we present the results of the new EOS based on the variational method together with those of the previous FYSS EOS based on the RMF with the TM1 parameter set [26]. The free energy of free nucleons, bulk energy of heavy nuclei and self-energy shifts of light nuclei are different between the two models. The shell, surface and Coulomb energies of heavy nuclei are also affected by the replacement of bulk properties such as the saturation densities of heavy nuclei. Table 1 summarizes the following bulk properties of nuclear matter at zero temperature: the nuclear saturation density of symmetric nuclear matter n_s0, internal energy per baryon of symmetric matter at the saturation point E₀, compressibility K, symmetry energy J, and slope parameter in the density dependence of symmetry energy L. The variational EOS is softer with smaller values of K, J and L [27] although both EOSs can support the cold neutron stars of masses ∼2 M_⊙ limit set by recent observations [63, 64]. On the other hand, radii of 1.4 M_⊙ neutron stars are quite different and ∼14.2 km and 11.5 km, for the RMF and variational EOSs, respectively [27].

Figure 1 compares the saturation densities n_si of heavy nuclei as a function of the charge fraction ${Z}_{i}/{A}_{i}$ at zero and critical temperatures, the latter of which correspond to the temperatures where the saturation disappears. The larger values of J and L in the FYSS EOS lead to much lower saturation densities for neutron-rich nuclei compared with the new EOS. Not only the bulk energy but also the Coulomb and surface energies depend on n_si.

Zoom In Zoom Out Reset image size

Figure 2 displays the individual values of the bulk, shell, Coulomb, surface and binding energies per baryon of nuclei with A_i = 100 at two temperatures. Note that the sign is changed in the former two and the baryon rest mass is subtracted from the bulk energy, which is defined as ${E}_{i}^{^{\prime} {\rm{B}}}$, and the binding energy per baryon is expressed as ${B}_{i}=-({E}_{i}^{^{\prime} {\rm{B}}}+{E}_{i}^{{\rm{C}}}+{E}_{i}^{{\rm{Sf}}}+{E}_{i}^{{\rm{Sh}}})/{A}_{i}$. The larger n_si in the new EOS results in smaller surface energies and larger Coulomb energies of heavy nuclei. Neutron-rich nuclei in the FYSS EOS have smaller absolute values of the bulk energy due to the larger values of J and L. At the higher temperature ($T\gtrsim 3$ MeV), absolute values of the bulk energy become larger due to the entropy contribution from nucleons in the nucleus, which increases as the saturation density is reduced, especially for the FYSS EOS.

Zoom In Zoom Out Reset image size

Below we present the results for T = 1, 3 and 5 MeV and ${Y}_{p}=0.2$ and 0.4, typical values in the supernova core. Figure 3 shows the mass fractions of free protons and neutrons, and of light and heavy nuclei. For the neutron-rich conditions of Y_p = 0.2, the mass fraction of heavy nuclei in the new EOS is larger than that in the FYSS EOS. This is because neutron-rich nuclei have larger binding energies or smaller masses. As for Y_p = 0.4, the difference between the two EOSs is small compared with the neuron-rich conditions, since the bulk properties in the asymmetric nuclear matter are not much different between them. We find that the mass fraction of light nuclei decreases at some density and the heavy nuclei take their place at higher densities in both EOSs. The mass fraction of free neutrons in the new EOS tends to be smaller than in the FYSS EOS due to the smaller value of J in the former, which leads to large populations of neutron-rich nuclei. The mass fractions of free protons, on the other hand, are almost the same between the two EOSs.

Zoom In Zoom Out Reset image size

Figure 4 displays the critical lines, where the mass fractions of heavy nuclei $({Z}_{i}\geqslant 6)$, ${X}_{A}={\sum }_{i}{A}_{i}{n}_{i}/{n}_{{\rm{B}}}$ and of light nuclei $({Z}_{j}\lt 6)$, ${X}_{a}={\sum }_{j}{A}_{j}{n}_{j}$, become 10⁻⁴. At low temperatures ($T\lesssim 3\,\mathrm{MeV}$ ), heavy nuclei appear at lower densities in the new EOS than in the FYSS EOS whereas, at high temperatures, the order of the critical densities becomes opposite, since the RMF we use yields larger absolute values of the bulk energy, which originate from the entropy contributions of nucleons in nuclei. On the other hand, both EOSs give similar critical densities in figure 4, at which light nuclei show up. Note also the densities, at which nuclear pastas melt into uniform nuclear matter, depend on the saturation densities and those in the FYSS EOS are smaller.

Zoom In Zoom Out Reset image size

Figure 5 shows the mass fractions of the isotopes of Ni (Z_i = 28) and Sn (Z_i = 50). In all cases considered there, larger amounts of neutron-rich isotopes appear in the new EOS than in the FYSS EOS due to the smaller symmetry energy. Figure 6, on the other hand, displays the mass number distributions of heavy nuclei for four different conditions. At ${n}_{{\rm{B}}}={10}^{-3}$ fm⁻³, T = 1 MeV and Y_p = 0.2, nuclei with large mass numbers are more abundant in the new EOS, since the smaller symmetry energy allows the nuclei in the vicinity of the neutron magic number N = 126 to be more populated, which can also be confirmed in figure 5. For the other three conditions, the mass distribution in the new EOS shifts slightly leftward compared with the FYSS EOS owing to the smaller surface and larger Coulomb energies that in turn result from the larger saturation densities in the former EOS. The difference between the two EOSs is particularly small at T = 1 MeV and Y_p = 0.4, which can be confirmed both in figures 5 and 6. Note that zig-zag features remarkable at T = 1 MeV, for which the shell effects are responsible, are almost washed out at T = 3 MeV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 exhibits the average mass numbers of heavy nuclei as a function of density and indicates that they are smaller in the new EOS than in the FYSS EOS in most cases due to the larger saturation densities in the former. At ${n}_{{\rm{B}}}\lesssim {10}^{-2}$ fm⁻³, T = 1 MeV and Y_p = 0.2, the result is the other way round because of the larger abundance of neutron-rich magic nuclei for the new EOS.

Zoom In Zoom Out Reset image size

Finally we compare thermodynamical quantities of relevance. Figures 8–10 present the free energy per baryon, baryonic pressure and entropy per baryon, respectively. Around the saturation density, the FYSS EOS based on the RMF gives much larger free energies and baryon pressures for all combinations of density, temperature and proton fractions considered there due to its stiffer bulk property. On the other hand, the difference of the free energies is not significant at low densities. At low temperatures and high densities, the pressures become negative, since the Coulomb energies of heavy nuclei is decreased by the rise of electron density [24]. The densities, at which the pressure becomes negative, are influenced by the Coulomb energy and the nuclear composition and are hence a bit different between the two EOSs. Note that the total pressure of supernova matter is always positive owing to electron's contribution [21]. The entropy per baryon is slightly affected by the difference in the mass fractions of nucleons and nuclei as displayed in figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Implementing the new EOS in our neutrino radiation-hydrodynamics code, we have conducted spherically symmetric simulations of core-collapse for a $11.2\ {M}_{\odot }$ supernova progenitor [65]. We have paid particular attention to the difference from the FYSS EOS [26], with which we perform a comparable simulation with the same setting. We designate the simulation with the new EOS as the 'VM' model and the one with the FYSS EOS as the 'RMF' model.

Our code solves simultaneously the Boltzmann equations for neutrino transport, hydrodynamic equations for matter motions and Poisson equations for Newtonian gravity, since they are all coupled with one another. We follow the dynamics up to $150\ \mathrm{ms}$ after core bounce. Neutrino-matter interactions are treated in the same way as in Nagakura et al [66, 67], which include updated sets of weak interactions and of electron capture (EC) rates for heavy nuclei among other things [68–70]. They are calculated based on the nuclear populations obtained from a nuclear statistical equilibrium EOS and with the previous version of the FYSS EOS [25] and tabulated as functions of ρ, T and Y_p. Note that the same table is used for the simulations with the new and FYSS EOSs. It is true that the use of the EC rates based on the nuclear abundances derived with each EOSs is more consistent in these simulations, but we are concerned with the difference not in the EC rates but in the EOSs themselves in this paper. We therefore decided to use the same EC table in this study. The results of fully consistent core-collapse simulations will be reported elsewhere soon. For more technical details of our code, we refer readers to Nagakura et al [66].

It turns out that the overall dynamics are not much different between the two simulations. In particular, the shock radius is surprisingly similar to each other as a function of time even though the distributions of ρ, T and Y_p are not so identical as described in detail below. This seems to be an accident, with different effects canceling each other. The use of the same EC rates may also be contributing to this similarity. In neither case, we obtain an explosion as expected for spherical models. If one looks into details, however, some differences become apparent, to which we shift our attention hereafter.

In figure 11, we display radial profiles of several important quantities as a function of mass coordinate at three different times in the collapse phase, i.e., when the central density (${\rho }_{c}$) reaches ${10}^{11},{10}^{12}$ and ${10}^{13}\ {\rm{g}}\,{\mathrm{cm}}^{-3}$. In these figures, we can immediately recognize the differences in the temperature and entropy profiles in the central region between the two models, which are particularly significant after the central density reaches ${10}^{12}\ {\rm{g}}\,{\mathrm{cm}}^{-3}$. This higher entropy for the new EOS originates from smaller mass fractions of heavy nuclei in the core matter with ${Y}_{p}\sim 0.46$ at the beginning of the collapse. The reason is as follows: in the matter with heavy nuclei and α particles dominating over free nucleons, the mass fraction of heavy nuclei is approximately given as ${X}_{A}\sim (0.5-{Y}_{p})/(0.5-\bar{Z}/\bar{A})$, which is in turn obtained from the baryon number conservation ${X}_{A}+{X}_{\alpha }\sim 1$ and the charge neutrality $(\bar{Z}/\bar{A}){X}_{A}+0.5{X}_{\alpha }\sim {Y}_{e}={Y}_{p}$, where $\bar{Z}$ and $\bar{A}$ are the average atomic and mass numbers of heavy nuclei. For a given Y_p, neutron-rich nuclei are more abundant, implying a smaller $\bar{Z}/\bar{A}$, in the new EOS, because of the smaller symmetry energy as discussed in the previous section. This is borne out in the bottom left panel of figure 11. The matter around the mass coordinate of $1.0\ {M}_{\odot }$ is similar to the matter at the center at the initial time. It is evident that the VM model gives a smaller X_A and a larger ${X}_{a}\sim {X}_{\alpha }$, and hence a higher entropy per baryon although ρ, T and Y_p are almost identical between the two models. In the early phase of core-collapse (${\rho }_{c}\sim {10}^{10}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$), the temperature of the central region increases more rapidly by adiabatic compression in the VM model with the higher entropy than in the RMF model. The higher temperature accelerates electron captures and deleptonizes the core more quickly (see the radial profile of Y_e in figure 11). The smaller electron fraction decreases the degeneracy pressure of electrons, accelerating the collapse, raising the temperature further, and resulting in rather large deviations in temperature and entropy (see the corresponding panels). Note that the entropy is changed by electron captures and accompanied emissions of neutrinos.

Zoom In Zoom Out Reset image size

It is also important to note that such high temperatures and small electron fractions reduce the mass fraction of heavy nuclei further, enhancing the mass fractions of nucleons and light nuclei instead. The matter in the central region ($\lesssim 0.6\ {M}_{\odot }$) consists of heavy nuclei with $\bar{A}\gtrsim 60$, light nuclei and dripped neutrons when ${\rho }_{c}$ exceeds ∼10¹² g cm⁻³. For such matter with $T\lesssim 3\,\mathrm{MeV}$ and ${Y}_{p}(={Y}_{e})\lesssim 0.4$, the new EOS, in contrast to the previous case, tends to give greater values to X_A than the FYSS EOS, as shown in the previous section. Owing to the higher T and lower Y_p ($={Y}_{e}$) in the VM model than in the RMF model, however, the order is reserved, keeping X_A lower in the dynamical simulations of the VM model. The smaller mass fractions of heavy nuclei imply lower opacities from neutrino–nucleus coherent scatterings, the dominant source of opacity during collapse, meaning that the deleptonization is further accelerated in the VM model. This in turn makes the mass of the inner core smaller. As is well known, the small mass of the inner core is disadvantageous for the outward propagation of the shock wave generated by core bounce, since the gravitational energy liberated and hence the initial shock energy tend to be smaller, and the outer core, the obstacle for shock propagation, becomes larger [71, 72].

Figure 12 shows the radial profiles of the same quantities as in figure 11 but at the time of core bounce, which we define in this study as the time when the shock wave reaches the mass coordinate of $0.6\ {M}_{\odot }$. As discussed in the previous paragraph, the stronger deleptonization for the VM model gives the central electron fraction ∼0.25 at bounce, which is $\sim 10 \%$ smaller than that for the RMF model. We find, however, that the shock wave is stronger in the VM model than in the RMF, which can be understood from the greater jump of entropy at the shock wave for the VM model. This result may seem at odds with our previous statement that the shock energy tends to be lower for the smaller inner core. This is true only when the EOS is identical, though. In fact, the new EOS is softer than the FYSS EOS as mentioned earlier. Then the inner core becomes not only less massive but also more compact, which is corroborated in the top left panel of the figure. The central density at core bounce is indeed higher in the VM model. The smaller radius of the inner core means a larger liberation of gravitational energy. Indeed, this is sufficient to compensate for the smaller mass of inner core and leads to the similar trajectory of shock wave at later times in the two models (see the right panel in figure 13).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The left panel in figure 13 compares the time evolutions of the central density in the post-bounce phase between the two models. The new EOS has consistently higher values than the RMF model, an indication that the former is softer and the inner core is more compact. In spite of the rather large difference in the central density between the two models, the time evolutions of the shock radii are quite similar to each other (see the right panel in the same figure). It should be noted, however, that this coincidence might be an artifact of our use of the same EC rates for heavy nuclei, which is certainly inconsistent in the VM and RMF models, not correctly reflecting the differences in the nuclear composition. In fact, the outer cores have very similar deleptonization histories during the collapsing phase between the two models (see the panel for Y_e in figure 11), which means that the shock waves propagate through very similar outer cores. It is also important to note that the luminosity and mean energy of neutrinos are also almost identical between the two models (see figure 14), meaning that neutrino coolings are almost the same, which also contributes to the similar evolutions of shock radius.

Zoom In Zoom Out Reset image size

It is interesting to note that substantial amounts of heavy nuclei exist around the PNS (see the left panel in figure 15) and they survive for more than 100 ms after bounce (see the upper left panel in figure 16). We expect that they might affect the diffusion of neutrinos, since the coherent scattering by heavy nuclei gives far greater opacities than scatterings on free nucleons or light nuclei. Incidentally, there is a non-negligible mass fraction of light nuclei just behind the shock wave (see the right panel of figure 15). This may indicate that weak interactions between neutrinos and these light nuclei should not be ignored in the post-shock region and may play an important role in the timing of the neutrino burst. We will address this issue by incorporating the relevant weak interactions with light nuclei in simulations [14–16]. It should be noted that non-negligible mass fractions of light nuclei appear just outside the PNS. In fact, their mass fraction sometimes overwhelms that of nucleons. These light nuclei, both behind the shock wave and near the PNS, may hence have the potential to change the neutrino cooling/heating (see also [14]) and affect shock revival eventually. The detailed study on these effects will be presented in our forthcoming paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have constructed a new EOS for use in supernova simulations. It includes the full ensemble of nuclei and is derived from the EOS for uniform nuclear matter that employed realistic nuclear forces in the variational method. This should be contrasted with standard EOSs, which include a single representative nucleus and treat the nuclear interactions phenomenologically with either Skyrme type interactions or meson-exchange models. In the calculation of nuclear masses, we have taken into consideration various in-medium effects such as the smearing of shell effects in heavy nuclei and the self- and Pauli-energy shifts in light nuclei at high densities and/or temperatures.

For some representative conditions in supernovae, we have compared the new EOS and our previous one: the FYSS EOS, in which the RMF EOS with the TM1 parameter set is employed. At T = 1 MeV and Y_p = 0.4, typical values at the center in the initial phase of collapse, the difference between the nuclear compositions in the two EOSs is not so large, since normal nuclear statistical equilibrium is expected to prevail with negligible in-medium effects. For more neutron-rich conditions (${Y}_{p}\sim 0.2$), the smaller symmetry energy in the new EOS provides larger mass fractions of heavy nuclei than the FYSS EOS. At higher temperatures ($T\gtrsim 3\,\mathrm{MeV}$), the larger entropy contribution leads to larger total mass fraction of heavy nuclei in the FYSS EOS than in the new EOS. It is also found that the new EOS tends to have larger mass fractions of nuclei with small mass numbers and charge fractions owing to the larger saturation densities and the smaller symmetry energy. The abundance of light nuclei, on the other hand, is not much changed by the replacement of the previous EOS with the new one although it is somewhat affected by the difference in the mass fractions of heavy nuclei. This is because in the new EOS the self- and Pauli-energy shifts are not evaluated consistently with the EOS for uniform nuclear matter.

Employing the two EOSs, we have performed dynamical simulations of core-collapse under the assumption of spherical symmetry. We have found that the new EOS gives smaller mass fractions of heavy nuclei and higher entropies for conditions typical in the early phase of core-collapse. This is due to the lower symmetry energy in the new EOS, which leads to slightly lower values of the charge fraction. Owing to higher entropies, the temperature rises more quickly in the simulation with the new EOS. The higher temperatures then enhance electron captures, resulting in lower electron fractions and even higher temperatures as the collapse proceeds. As a result, the mass fraction of heavy nuclei continues to be lower at later times in the collapsing phase, where Y_e becomes much lower with the new EOS than with the FYSS EOS in spite of the fact that the new EOS tends to give larger mass fractions of heavy nuclei in neutron-rich matter. The reduction of electron fraction implies a smaller mass of the inner core.

The lower electron fractions obtained in the simulation with the new EOS imply that the mass of the inner core is smaller. On the other hand, the inner core is more compact for the new EOS, with the central density being higher. This is due to the fact that the new EOS is softer than the FYSS EOS around nuclear density. These two effects, i.e., the smaller mass and the higher density of the inner core tend to work oppositely for the initial strength of the shock wave generated by core bounce, with the former (latter) reducing (enhancing) the shock strength. As a consequence, the shock is a bit stronger initially for the new EOS than for the FYSS EOS. It turns out, however, that the ensuing outward propagations of the shock wave in the outer core are not much different between the two EOSs. This may be an artifact, though, as the same EC rates for heavy nuclei were used in spite of different compositions in the two simulations. Indeed the deleptonization may be more suppressed for the new EOS, since neutron-rich nuclei, which have smaller EC rates in general, are more abundant in the new EOS. The results of more consistent simulations will be reported in our forthcoming paper.

Very recently, three of the authors of this paper published a paper on another EOS for supernova simulations, which we refer to as the Togashi EOS [27]. It is derived from the same EOS obtained with the variational method for uniform nuclear matter by Togashi et al [29] and employs the SNA based on the Thomas–Fermi approximation. As expected, the differences between our new EOS and the Togashi EOS are very similar to those observed between the FYSS and STOS EOSs [18–20]. Details are discussed in Furusawa et al [24]. As a matter of fact, the latter two EOSs employ the same RMF EOS for uniform nuclear matter but at sub-nuclear densities the former handles the ensemble of nuclei just as the new EOS does whereas the latter resorts to the SNA that the Togashi EOS adopts. The multi-nuclear species EOSs, i.e., the new EOS and the FYSS EOS, commonly predict smaller atomic and mass numbers than SNA EOSs, or the Togashi EOS and the STOS EOS [24, 31, 32]. In the new EOS and the FYSS EOS, shell effects are taken into account and, as a result, the atomic and mass numbers grow with density in a step-wise fashion while they increase smoothly in the STOS and Togashi EOSs, in which shell effects are neglected. Furthermore, the entropy per baryon is lower in the latter two EOSs because they ignore the translational motion of the representative nucleus [21, 24].

Although the major differences between the new EOS and the Togashi EOS are mostly ascribed to the SNA adopted in the latter, there are certainly other differences. In fact, the differences between the STOS and Togashi's EOSs are not always consistent with those between the FYSS EOS and the new EOSs. For instance, the Togashi EOS predicts larger mass numbers than the STOS EOS for neutron-rich conditions (${Y}_{p}\sim 0.1$) because of the larger saturation density in the former, which leads to greater surface energy described by gradient terms in the energy functional employed in the Thomas–Fermi calculation. On the other hand, the new EOS gives smaller mass numbers than the FYSS EOS more often than not although there is no general trend in part due to the shell effects. This happens because we use the same coefficient for surface tensions, ${\sigma }_{0}$, in the new EOS and the FYSS EOS, which was originally adopted in the LDM [17]. Then, the surface energy becomes smaller in the new EOS with the larger saturation density. It is true that the surface tension should be treated in a consistent manner with the EOS for uniform nuclear matter but the difference that the surface energy would make is far smaller than that stemming from the SNA approximation discussed in the previous paragraph.

Note that the mass models of heavy and light nuclei employed in the new EOS still have some room for improvement, since there are uncertainties in our phenomenological treatment of in-medium effects. In addition, in the new EOS we have adopted the same parameters in the formulae for the shell and surface energies of heavy nuclei as well as for the Pauli energies of light nuclei as those employed in the FYSS EOS. They should be changed according to the theory used to describe uniform nuclear matter. The light clusters, especially deuterons, in the new EOS may have to be described in the framework of the variational method by using realistic nuclear forces. Nevertheless, we believe that the new EOS is one of the most elaborate EOSs constructed so far for supernova simulations. The tabulated new EOS is available in the public domain¹⁰ .

SF and HN were supported by Japan Society for the Promotion of Science Postdoctoral Fellowships for Research Abroad. Some numerical calculations were carried out on the PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work was supported by the RIKEN iTHES Project and in part by the usage of supercomputer systems through the Large Scale Simulation Program (Nos. 15/16/-08,16/17-11) of High Energy Accelerator Research Organization (KEK) and Post-K Projects (hp 150225, hp160071, hp160211,hp170031,hp170230) at K-computer, RIKEN AICS as well as the computational resources provided by RCNP at Osaka University, YITP at Kyoto University, University of Tokyo and JLDG. This work was supported by a Grant-in-Aid for the Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan (24103006, 24244036, 16H03986, 15K05093, 24105008, 25400275, 26104006).