Finite-temperature Extension for Cold Neutron Star Equations of State

Many dynamical phenomena—including core-collapse supernovae, the formation and subsequent cooling of proto-neutron stars, and both the electromagnetic and gravitational signals from neutron star mergers—depend sensitively on the neutron star equation of state (EOS) at densities where the EOS is not well understood. In addition, for these dynamical phenomena, there are two further complications. First, temperatures may range from below the Fermi temperature, for which "cold" EOS would suffice, to temperatures of up to 10–100 MeV in neutron star mergers (e.g., Oechslin et al. 2007). Second, the composition may range from nearly pure neutron matter to symmetric matter, with some dynamical timescales that are shorter than the timescale required to establish β-equilibrium. While astrophysical observations of stationary neutron stars probe the cold EOS in β-equilibrium and laboratory experiments constrain the hot EOS of symmetric matter, extrapolations between the two regimes remain difficult. (For a schematic representation of these various regimes, see Figure 1. For recent reviews, see e.g., Lattimer & Prakash 2016; Özel & Freire 2016.) Such extrapolations to arbitrary proton fraction and temperature add further uncertainty to the EOS and complicate numerical simulations of these phenomena.

Zoom In Zoom Out Reset image size

A large number of EOS have been calculated in the zero-temperature limit, ranging from purely nucleonic models (e.g., Baym et al. 1971; Friedman & Pandharipande 1981; Akmal et al. 1998; Douchin & Haensel 2001) to models incorporating quark degrees of freedom using state-of-the-art results from perturbative QCD (e.g., Fraga et al. 2014). Laboratory experiments and neutron-star observations do not yet have sufficient power to distinguish between these models. Furthermore, it is likely that these EOS do not span the full range of possible physics. This possibility has motivated the creation of a large number of parametric EOS, which were first introduced by Read et al. (2009) and Özel & Psaltis (2009). These parametric models do not require a priori knowledge of the high-density nuclear physics governing the EOS and, hence, can be used to probe unknown physics from neutron star observations.

A much smaller number of EOS that self-consistently incorporate finite-temperature effects have been calculated to date. Among the most well-known of these are the LS model, which is based on finite-temperature compressible liquid drop theory with a Skyrme nuclear force (Lattimer & Swesty 1991) and Shen et al.'s (1998) EOS, which was calculated using relativistic mean field (RMF) theory with a Thomas–Fermi approximation. More recently, the statistical model developed by Hempel & Schaffner-Bielich (2010) has been applied to an additional ∼10 combinations of RMF models and nuclear mass tables.

Just as parameterizations of the cold EOS have proven to be useful in representing a broader range of physics, so too would a parametric finite-temperature EOS be useful for incorporating EOS effects into supernova and merger calculations. To this end, many authors have employed so-called "hybrid EOS," in which a thermal component for an ideal fluid is added to an arbitrary cold EOS to account for heating (Janka et al. 1993). The ideal-fluid thermal component is parameterized in terms of a simple adiabatic index as ${P}_{\mathrm{th}}={\epsilon }_{\mathrm{th}}({{\rm{\Gamma }}}_{\mathrm{th}}-1)$, where ${P}_{\mathrm{th}}$ and ${\epsilon }_{\mathrm{th}}$ are the thermal pressure and energy density and Γ_th is the adiabatic index, which is assumed to be constant. Although this approach is computationally simple, it neglects the effect of degeneracy on the thermal pressure. At high densities and finite temperatures, part of the available energy acts to lift degeneracy rather than contributing additional thermal support. This causes a net reduction in the thermal pressure at high densities, compared to the prediction for an ideal fluid.

The density-dependence of these thermal effects depends directly on the density-dependence of the nucleon effective mass, as has been shown for many EOS (Constantinou et al. 2014, 2015). Constantinou et al. (2015) performed a Sommerfeld expansion to approximate the thermal properties at next-to-leading order and showed that the expansion terms require both the effective mass and its derivatives. Given a complete expression for the density-dependence of the effective mass, they showed that this formalism can be used to accurately approximate the thermal properties of a wide variety of EOS. Constantinou et al. (2017) later expanded this work and showed that the formalism can even be used to recreate models beyond mean field theory, such as Zhang & Prakash's (2016) two-loop exchange model.

The strong dependence of thermal properties on the effective mass can also be seen in the behavior Γ_th. For example, Constantinou et al. (2015) compared two EOS with similar zero-temperature properties but with different single-particle potentials, and hence different density-dependencies in their nucleon effective masses. They found substantially different thermal properties for the two EOS and that a constant Γ_th model failed to describe either EOS. Zhang & Prakash (2016) also found a strong density-dependence in Γ_th for their two-loop exchange model. These results indicate that Γ_th has a significant density-dependence for a diverse range of analytic models, which is not captured in the constant Γ_th approximation of the hybrid EOS.

Neglecting the effect of degeneracy on the thermal pressure also has important consequences for dynamical simulations. For example, Bauswein et al. (2010) compared the properties of a neutron star–neutron star merger that would be predicted by a hybrid EOS and by more realistic EOS. Specifically, they compared Shen et al.'s (1998) and Lattimer & Swesty's (1991) EOS with hybrid EOS that were constructed from the zero-temperature versions of these same EOS, with either ${{\rm{\Gamma }}}_{\mathrm{th}}=1.5$ or 2. They found that using the hybrid EOS predicts post-merger frequencies from a hypermassive neutron star that are 50–250 Hz smaller than those found with a realistic finite-temperature EOS. Moreover, the lifetime of the hypermassive remnant can deviate by a factor of two from the more realistic value. In addition, the post-collapse accretion disk mass around the resulting black hole can differ by up to 30% when the simplified thermal effects are used (Bauswein et al. 2010). These results suggest that it is indeed important to account for the effect of degeneracy on the thermal pressure when simulating neutron star mergers.

Constantinou et al.'s (2015) Sommerfeld expansion results can be used to explicitly correct a hybrid EOS to include degenerate effects, as long as the particle interactions and potentials of the cold EOS are known. However, requiring knowledge of the potentials of the cold EOS renders these corrections inapplicable to piecewise-polytropic EOS or other parametric forms of the EOS that are agnostic in their descriptions of the microphysics.

The goal of this paper is to develop a physically motivated framework to incorporate the thermal pressure that maintains the wide applicability of the hybrid EOS approach. With such a model, it will be possible to robustly add thermal effects to any cold EOS in β-equilibrium, without having to make the simplifying assumptions of an ideal fluid at all densities. Although the framework that we present in this paper is specific to neutron–proton–electron (n–p–e) matter, it could be generalized to include more exotic particles. We also include a symmetry-energy dependent correction that extrapolates the proton fraction away from β-equilibrium. The complete model thus allows us to build an EOS at finite-temperature and arbitrary proton fraction from any cold n–p–e EOS in neutrinoless β-equilibrium, including piecewise-polytropic EOS. Moreover, the model is analytic and in closed-form, and thus can be calculated efficiently in dynamical simulations.

We start in Section 2 with a brief review of existing finite-temperature EOS and a discussion of the regimes in which thermal effects become important. In Section 3, we outline our model. We provide the symmetry-energy dependent extrapolation to arbitrary proton fraction in Section 4. In Section 5, we introduce our M*-approximation of the thermal effects. We summarize the model in Section 6, in which all of the relevant equations can be found in Boxes I and II. Finally, we quantify the performance of our model in Section 7. We find that with a relatively small set of parameters, our complete model is able to recreate existing finite-temperature EOS with introduced errors of ≲20% for densities above the nuclear saturation density.

Before introducing our new approximation for the pressure at arbitrary proton fraction and temperature, we will briefly review the finite-temperature EOS that have been previously developed.

Two of the most widely used finite-temperature EOS are the models of Lattimer & Swesty (1991, hereafter LS), which is based on a finite-temperature liquid drop model with a Skyrme nuclear force, and Shen et al. (1998, hereafter STOS), which is a RMF model that is extended with the Thomas–Fermi approximation. An additional eight EOS have been calculated with Hempel & Schaffner-Bielich's (2010, hereafter, HS) framework, which is a statistical model that consists of an ensemble of nuclei and interacting nucleons in nuclear statistical equilibrium and, hence, goes beyond the single nucleus approximation that both LS and STOS assume. Each HS EOS represents the nucleons with a RMF model and additionally includes excluded volume effects. Of the RMF models that have been used with the HS method, six are nucleonic: TMA (Toki et al. 1995), TM1 (Sugahara & Toki 1994), NL3 (Lalazissis et al. 1997), FSUGold (Todd-Rutel & Piekarewicz 2005), IUFSU (Fattoyev et al. 2010), DD2 (Typel et al. 2010); while the models BHBΛϕ and BHBΛ include hyperons with and without the repulsive hyperon–hyperon interaction mediated by the ϕ meson, respectively (Banik et al. 2014). Additionally, Steiner et al. (2013) created a set of two finite-temperature EOS, SFHo/x, that also used the statistical method of HS but which included new RMF parameterizations and constraints from neutron star observations. There are also the EOS of G. Shen, which are based on a virial expansion and nuclear statistical equilibrium calculations at low densities and RMF calculations at high densities, using the FSUGold (Shen et al. 2011a) and NL3 (Shen et al. 2011b) models. Tables of these various EOS can be found on M. Hempel's website ¹ , at stellarcollapse.org, and/or the CompOSE database. ²

More recently, several new finite-temperature EOS have been added to the CompOSE database. These include the SLY4-RG model, which is calculated in nuclear statistical equilibrium using a Skyrme energy functional (Gulminelli & Raduta 2015; Raduta & Gulminelli 2019); chiral mean-field theory models, which include hyperons as additional degrees of freedom (e.g., Dexheimer 2017); generalized relativistic density functional models (e.g., Typel 2018); and models calculated using a variational method applied to two- and three-body nuclear potentials (e.g., Togashi et al. 2017).

For the sake of simplicity, in the following analysis we will focus on a subset of these EOS and will include only models that are nucleonic. In particular, our sample will include STOS and the eight nucleonic EOS that are calculated with the HS method, to represent the models based on RMF theory. We will also include LS (with a compression modulus K = 220 MeV) and SLY4-RG to represent non-relativistic models with Skyrme nuclear forces.

In spite of the increasing number of finite-temperature EOS that have been calculated, they nevertheless span a relatively limited range of physics, especially when compared to the diversity of cold EOS models. In order to span a broader range of possible physics, many authors have used the so-called "hybrid EOS," which assume that the thermal pressure is given simply by an ideal-fluid term that can be added to any cold EOS. The hybrid EOS were first introduced by Janka et al. (1993) and have been used in many subsequent works (for recent reviews, see Shibata & Taniguchi 2011; Faber & Rasio 2012; Baiotti & Rezzolla 2017; Paschalidis & Stergioulas 2017). In these hybrid EOS, the thermal pressure is written as

$$\begin{eqnarray}&&{P}_{\mathrm{th},\mathrm{hybrid}}(n,T)={{nE}}_{\mathrm{th},\mathrm{hybrid}}(n,T)({{\rm{\Gamma }}}_{\mathrm{th}}-1),\end{eqnarray} \tag{ 1 }$$

where ${E}_{\mathrm{th},\mathrm{hybrid}}(n,T)$ is the thermal contribution to the energy per baryon, n is the baryon number density, and Γ_th is the thermal adiabatic index and is constrained to be $1\leqslant {{\rm{\Gamma }}}_{\mathrm{th}}\leqslant 2$. In the hybrid approximation, Γ_th is assumed to be constant.

Following Etienne et al. (2008), the hybrid temperature-dependence of ${E}_{\mathrm{th},\mathrm{hybrid}}$ is included as an ideal fluid plus a contribution from relativistic particles; i.e.,

$$\begin{eqnarray}&&{E}_{\mathrm{th},\mathrm{hybrid}}(n,T)=\displaystyle \frac{3}{2}{k}_{{\rm{B}}}T+\displaystyle \frac{4\sigma }{c}\displaystyle \frac{{f}_{s}}{n}{T}^{4},\end{eqnarray} \tag{ 2 }$$

where k_B is the Boltzmann constant, T is the temperature, and $\sigma \equiv {\pi }^{2}{k}_{{\rm{B}}}^{4}/[60{{\hslash }}^{3}{c}^{2}]$ is the Stefan–Boltzmann constant, with ℏ the Planck constant and c the speed of light. The parameter f_S represents the number of ultra-relativistic species that contribute to the thermal pressure. For ${k}_{{\rm{B}}}T\ll 2{m}_{e}{c}^{2}$, where m_e is the mass of an electron, photons will dominate and f_S = 1. For ${k}_{{\rm{B}}}T\gg 2{m}_{e}{c}^{2}$, electrons and positrons become relativistic as well and yield ${f}_{S}=1+2\times (7/8)=11/4$. Finally, for ${k}_{{\rm{B}}}T\,\gtrsim 10\,\mathrm{MeV}$, thermal neutrinos and anti-neutrinos appear, rendering ${f}_{s}=11/4+3\times (7/8)=43/8$. If right-handed neutrinos were to exist, then this would become ${f}_{s}\,=11/4+3\times 2\times (7/8)=8$.

We note that all 12 EOS discussed above neglect neutrinos in their calculations. The STOS EOS additionally neglects leptons and photons, which we add wherever we use STOS in this paper. For the STOS thermal lepton and photon contribution, we use Equation (2) with the appropriate lepton density. For the cold lepton energy, we add the contribution for a degenerate gas of relativistic electrons. Because all of the EOS neglect neutrinos, we will also neglect neutrinos in our comparisons and thus we will calculate f_S only as

$$\begin{eqnarray}{f}_{S}=\left\{\begin{array}{ll}1, & {k}_{{\rm{B}}}T\lt 1\,\mathrm{MeV},\\ 11/4, & {k}_{{\rm{B}}}T\geqslant \,1\,\mathrm{MeV}.\end{array}\right.\end{eqnarray} \tag{ 3 }$$

We, therefore, account for the degrees of freedom introduced by the possible presence of ultra-relativistic positrons. However, throughout this paper, we will assume that the population of positrons is small and that their contribution to the pressure or energy at higher densities is negligible. If there were a scenario in which the population of positrons were significant compared to the electrons, then one would have to explicitly account for the positrons in particle-counting as well as in imposing charge neutrality.

In order to highlight the regimes where a realistic finite-temperature EOS and the hybrid approximation differ, we show a phase diagram in Figure 2. In this plot, we show various regions calculated for the EOS STOS at a fixed proton fraction of Y_p  = 0.1. The total pressure, P_total, is thus calculated at Y_p  = 0.1 and a given temperature. The cold contribution, P_cold, is calculated at the same Y_p and at zero-temperature. ³ Finally, the thermal contribution, P_th, is defined as ${P}_{\mathrm{total}}-{P}_{\mathrm{cold}}$ for the same proton fraction.

Zoom In Zoom Out Reset image size

In this figure, the blue-shaded region shows the regime where the total pressure is dominated by the cold pressure; there, the thermal pressure of STOS contributes <1% of the total pressure. The red-shaded region represents the regime where the thermal pressure can be approximated by the ideal-fluid pressure (${P}_{\mathrm{th},\mathrm{ideal}}={{nK}}_{{\rm{B}}}T$), to within 1%. The white region between these two extremes represents the range of parameter space in which the thermal pressure is important but the ideal-fluid approximation does not yet apply. In this white region, the effects of degeneracy on the thermal pressure cannot be neglected.

For comparison, Figure 2 also shows the projected temperature–density profiles from three different simulations of relevant astrophysical phenomena. The green line shows the profile of a hypermassive neutron star remnant 12.1 ms after the merger of two 1.35 M_⊙ neutron stars, as simulated using the EOS STOS (Sekiguchi et al. 2011). The orange and purple lines both come from numerical simulations of the evolution of a proto-neutron star using a bulk version of the LS EOS. The orange line gives the profile of the proto-neutron star at 200 ms after the core bounce, while the purple line shows the profile of the proto-neutron star at the end of the de-leptonization phase (Camelio et al. 2017). We note that these profiles are not necessarily calculated at Y_p  = 0.1, but we include them nevertheless to show the approximate relevant temperatures and densities for these phenomena.

In order to further explore the dependence on the proton fraction, we also calculated the regime where degeneracy dominates for increasing values of Y_p . We find that as the proton fraction increases toward Y_p  = 0.5, the white degeneracy region in Figure 2 shrinks but still largely encompasses the shown profiles. We thus find that all of these simulations primarily probe the phase space where degenerate thermal effects are important. This suggests that using the hybrid approximation instead of the full thermal pressure may bias the outcomes from such simulations.

In order to construct a finite-temperature EOS at arbitrary proton fraction, our model must be able to extrapolate from β-equilibrium to an arbitrary Y_p , as well as from cold matter to an arbitrary temperature. This will naturally introduce dials into our model that can be adjusted to represent a wide range of physics, based on the symmetry energy, its slope, and the strength of particle interactions that we wish to include. Moreover, we will show that with a small set of parameters, the EOS that are currently in use in the literature can be replicated to high accuracy.

We start with our model in general terms, for which we will derive analytic expressions in the following sections. Our final model will be for the complete energy per baryon, $E(n,{Y}_{p},T)$, separated into analytic, physically motivated terms. A summary of the final equations can be found in Boxes I and II in Section 6.

We can expand the energy per particle of nuclear matter, E_nucl, about the neutron excess parameter, $(1-2{Y}_{p})$, to second order as

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{nucl}}(n,{Y}_{p},T) & = & {E}_{\mathrm{nucl}}(n,{Y}_{p}=1/2,T)\\ & & +\,{E}_{\mathrm{sym}}(n,T){\left(1-2{Y}_{p}\right)}^{2},\end{array}\end{eqnarray} \tag{ 4 }$$

where ${E}_{\mathrm{nucl}}(n,{Y}_{p}=1/2,T)$ represents the energy of symmetric nuclear matter and

$$\begin{eqnarray}&&{E}_{\mathrm{sym}}(n,T)\equiv \displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}{E}_{\mathrm{nucl}}(n,{Y}_{p},T)}{\partial {\left(1-2{Y}_{p}\right)}^{2}}{| }_{{Y}_{p}=1/2}\end{eqnarray} \tag{ 5 }$$

is the symmetry energy. The proton fraction is related to the overall baryon density, n, according to

$$\begin{eqnarray}&&{Y}_{p}=\displaystyle \frac{{n}_{p}}{n}=\displaystyle \frac{{N}_{p}}{{N}_{n}+{N}_{p}},\end{eqnarray} \tag{ 6 }$$

where n_p is the proton density, N_p is the total number of protons, and N_n is the total number of neutrons. Throughout this paper, we enforce charge neutrality, which requires that the proton and electron densities balance. Thus, the electron density, n_e , can be written as

$$\begin{eqnarray}&&{n}_{e}={Y}_{p}n.\end{eqnarray} \tag{ 7 }$$

Finally, by requiring that the baryonic components combine to give the total density n, we can write the neutron density as

$$\begin{eqnarray}&&{n}_{n}=(1-{Y}_{p})n.\end{eqnarray} \tag{ 8 }$$

We can further expand Equation (4) by separating the energy of cold, symmetric matter from its thermal contribution; i.e.,

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{nucl}}(n,{Y}_{p},T) & = & {E}_{\mathrm{nucl}}(n,{Y}_{p}=1/2,T=0)\\ & & +\,{E}_{\mathrm{nucl},\mathrm{th}}(n,{Y}_{p}=1/2,T)\\ & & +\,{E}_{\mathrm{sym}}(n,T){\left(1-2{Y}_{p}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 9 }$$

Here and throughout this paper, we use the subscript "th" to indicate the thermal contribution to a variable, after the cold component has been subtracted.

In order to write the energy with respect to a cold EOS in β-equilibrium, as is often most relevant to start from in the study of neutron stars, we eliminate the cold, symmetric term in Equation (9) to yield

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{nucl}}(n,{Y}_{p},T) & = & {E}_{\mathrm{nucl}}(n,{Y}_{p,\beta },T=0)\\ & & +\,{E}_{\mathrm{nucl},\mathrm{th}}(n,{Y}_{p}=1/2,T)\\ & & +\,{E}_{\mathrm{sym}}(n,T){\left(1-2{Y}_{p}\right)}^{2}\\ & & -\,{E}_{\mathrm{sym}}(n,T=0){\left(1-2{Y}_{p,\beta }\right)}^{2},\end{array}\end{eqnarray} \tag{ 10 }$$

where ${Y}_{p,\beta }$ represents the proton fraction of a zero-temperature system in β-equilibrium. We note that the proton fraction depends on the density, i.e., ${Y}_{p,\beta }={Y}_{p,\beta }(n)$, but for simplicity we suppress this in our notation.

Finally, we must add the contribution of leptons and photons to this expression. The zero-temperature energy from relativistic degenerate electrons is given by

$$\begin{eqnarray}&&{E}_{\mathrm{lepton}}(n,{Y}_{p},T=0)=3{{KY}}_{p}{\left({Y}_{p}n\right)}^{1/3},\end{eqnarray} \tag{ 11 }$$

where the extra factor of Y_p comes from our definition of E as the energy per baryon, combined with Equations (6) and (7). Here, $K\equiv {(3{\pi }^{2})}^{1/3}({\hslash }c/4)$. Additionally, there will also be a thermal contribution, ${E}_{\mathrm{lepton},\mathrm{th}}(n,{Y}_{p},T)$, which we derive in Section 5.

Thus, our skeletal model for the total energy is given by the following set of equations:

$$\begin{eqnarray}&&E(n,{Y}_{p},T)=E(n,{Y}_{p},T=0)+{E}_{\mathrm{th}}(n,{Y}_{p},T)\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&\begin{array}{l}E(n,{Y}_{p},T=0)=E(n,{Y}_{p,\beta },T=0)\\ \,+\,{E}_{\mathrm{sym}}(n,T=0)\left[{\left(1-2{Y}_{p}\right)}^{2}-{\left(1-2{Y}_{p,\beta }\right)}^{2}\right]\\ \,+\,3K\left({Y}_{p}^{4/3}-{Y}_{p,\beta }^{4/3}\right){n}^{1/3}\end{array}\end{eqnarray} \tag{ 12b }$$

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{th}}(n,{Y}_{p},T) & = & {E}_{\mathrm{nucl},\mathrm{th}}(n,{Y}_{p}=1/2,T)\\ & & +\,{E}_{\mathrm{lepton},\mathrm{th}}(n,{Y}_{p},T)\\ & & +\,{E}_{\mathrm{sym},\mathrm{th}}(n,T){\left(1-2{Y}_{p}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 12c }$$

From these relations, we can derive the pressure via the standard thermodynamic relation,

$$\begin{eqnarray}&&P\equiv -{\left.\displaystyle \frac{\partial U}{\partial V}\right|}_{{N}_{q},S}{\left.={n}^{2}\left[\displaystyle \frac{\partial E(n,T=0)}{\partial n}\right]\right|}_{{Y}_{p},S}\end{eqnarray} \tag{ 13 }$$

where U is the total energy, V is the volume, N_q is the number of each species q, and S is the total entropy. From Equation (6), it is clear that evaluating these derivatives at constant N_q is equivalent to evaluating them at constant Y_p . In this paper, we will mainly plot results in terms of pressure. We summarize the complete expressions for pressure in Box II of Section 6.

While this set of expressions may seem to have a large number of terms, this separation allows these terms to be represented analytically. Moreover, as we will show later, the parameters of each term are linked directly to physics on which there are experimental constraints and of which further constraints are the motivation of many observations of astrophysical neutron stars: namely, the value of the symmetry energy at the saturation density, the slope of the symmetry energy, and the strength of interactions between particles.

We turn first to the symmetry energy correction term, ${E}_{\mathrm{sym}}(n,T)$ of Equation (4). The symmetry energy is defined as the per-nucleon difference in energy between symmetric matter and pure neutron matter. In other words, the symmetry energy represents the excess energy of matter with unequal numbers of protons and neutrons. In nuclear models, the symmetry energy is typically calculated as an expansion around the nuclear saturation density, for matter with Y_p  = 1/2. In Equation (4), we perform the expansion with respect to the proton fraction, and in the following section, we will introduce a density-dependence to extrapolate beyond the saturation density, where the coefficients of our approximation are experimentally constrained. In this section, we will provide the approximation for ${E}_{\mathrm{sym}}(n,T)$ at zero-temperature. For the thermal contribution to the symmetry energy, which turns out to be negligible, see Section 5.

It is particularly useful to parameterize the symmetry energy in terms of its separate kinetic and potential components at zero-temperature (e.g., Tsang et al. 2009; Steiner et al. 2010), modified by a parameter η to account for short-range correlations due to the tensor force acting between a spin-triplet or isospin-singlet proton–neutron pair. These correlations can significantly reduce the kinetic symmetry energy to even a negative value at the saturation density, compared to the kinetic energy of an uncorrelated Fermi gas model (Lovato et al. 2011; Vidaña et al. 2011; Xu & Li 2011; Carbone et al. 2012; Rios et al. 2014; Hen et al. 2015). In this framework, we parameterize the symmetry energy of Equation 12(b) as

$$\begin{eqnarray}&&{E}_{\mathrm{sym}}(n,T=0)=\eta {E}_{\mathrm{sym}}^{\mathrm{kin}}(n)+\left[{S}_{0}-\eta {E}_{\mathrm{sym}}^{\mathrm{kin}}({n}_{\mathrm{sat}})\right]{\left(\displaystyle \frac{n}{{n}_{\mathrm{sat}}}\right)}^{\gamma },\end{eqnarray} \tag{ 14 }$$

as in Li et al. (2015). Here, ${E}_{\mathrm{sym}}^{\mathrm{kin}}(n)$ is the "kinetic" symmetry energy, arising from the change in the Fermi energy of a gas at density n as the relative proton/neutron fraction changes, n_sat = 0.16 fm⁻³ is the nuclear saturation density, ⁴ and the second term represents the "potential" symmetry energy, which accounts for the interactions between particles. Because the exact form of the potential symmetry energy is not well known, it is anchored at the saturation density by the magnitude of the overall symmetry energy, ${S}_{0}\equiv {E}_{\mathrm{sym}}({n}_{\mathrm{sat}})$, and is given an arbitrary density-dependence through the constant γ.

In contrast, the kinetic energy term can be calculated directly from the nuclear momentum distribution. The kinetic energy of a free Fermi gas is given simply by

$$\begin{eqnarray}&&\displaystyle \frac{{\varepsilon }_{k,q}}{n}=\displaystyle \frac{3}{5}{E}_{{\rm{f}}}({n}_{q})\end{eqnarray} \tag{ 15 }$$

where ${\varepsilon }_{k,q}$ is the kinetic energy per particle, q represents the particle (either a neutron or proton), and E_f(n) is the Fermi energy,

$$\begin{eqnarray}&&{E}_{{\rm{f}}}({n}_{q})=\displaystyle \frac{{{\hslash }}^{2}}{2m}{\left(3{\pi }^{2}{n}_{q}\right)}^{2/3},\end{eqnarray} \tag{ 16 }$$

in which m is the mass of the relevant particle. For our approximation, we will neglect the small difference between the proton and neutron mass, and simply take m ≈ m_n , where m_n is the neutron mass.

By taking the difference between symmetric matter and pure neutron matter, the kinetic symmetry energy as a function of the total density is then

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{sym}}^{\mathrm{kin}}(n) & = & \displaystyle \frac{3}{5}\left[2{E}_{{\rm{f}}}\left({n}_{p}={n}_{n}=\displaystyle \frac{1}{2}n\right)-{E}_{{\rm{f}}}({n}_{n}=n)\right]\\ & = & \displaystyle \frac{3}{5}\left({2}^{1/3}-1\right){E}_{{\rm{f}}}(n).\end{array}\end{eqnarray} \tag{ 17 }$$

We can also eliminate the parameter η in Equation (14) by introducing the constant L, which is related to the overall slope at the saturation density via,

$$\begin{eqnarray}&&L\equiv 3{n}_{\mathrm{sat}}\left[\displaystyle \frac{\partial {E}_{\mathrm{sym}}(n,T=0)}{\partial n}\right]{| }_{{n}_{\mathrm{sat}}}.\end{eqnarray} \tag{ 18 }$$

Combining Equations (14) and (18), we can solve for η in terms of the quantities S₀ and L, which are constrained by nuclear physics experiments for matter near Y_p  = 1/2 (Lattimer & Lim 2013). We find

$$\begin{eqnarray}&&\eta =\displaystyle \frac{5}{9}\left[\displaystyle \frac{L-3{S}_{0}\gamma }{\left({2}^{1/3}-1\right)\left(2/3-\gamma \right){E}_{{\rm{f}}}({n}_{\mathrm{sat}})}\right],\end{eqnarray} \tag{ 19 }$$

thereby leaving one free parameter, γ, which is constrained by nuclear experiments to lie in the range ∼0.2 to 1.2 (see, e.g., Figure 2 of Li et al. 2015; Tsang et al. 2009).

We thus have a complete expression for the symmetry energy that depends only on the three parameters γ, S₀, and L, which, in principle, can be constrained by nuclear experiments. We can now use this functional form to fit for γ, by combining it with the following relationship between the symmetry energy and Y_p,β for charge-neutral n–p–e matter in neutrinoless β-equilibrium,

$$\begin{eqnarray}&&\displaystyle \frac{{Y}_{p,\beta }}{{\left(1-2{Y}_{p,\beta }\right)}^{3}}=\displaystyle \frac{64}{3{\pi }^{2}n}{\left[\displaystyle \frac{{E}_{\mathrm{sym}}(n,T=0)}{{\hslash }c}\right]}^{3}\end{eqnarray} \tag{ 20 }$$

(for a derivation of this relation, see, e.g., Blaschke et al. 2016, or Appendix A). When solved for ${Y}_{p,\beta }$, this becomes

$$\begin{eqnarray}\begin{array}{rcl}{Y}_{p,\beta } & = & \displaystyle \frac{1}{2}+\displaystyle \frac{{\left(2{\pi }^{2}\right)}^{1/3}}{32}\displaystyle \frac{n}{\xi }\\ & & \times \,\left\{\Space{0ex}{1.70em}{0ex}{\left(2{\pi }^{2}\right)}^{1/3}\right.-\left.\displaystyle \frac{{\xi }^{2}}{n}{\left[\displaystyle \frac{{\hslash }c}{{E}_{\mathrm{sym}}(n,T=0)}\right]}^{3}\right\},\end{array}\end{eqnarray} \tag{ 21 }$$

where, for simplicity, we have introduced the auxiliary quantity ξ, defined as

$$\begin{eqnarray}\begin{array}{rcl}\xi & \equiv & {\left[\displaystyle \frac{{E}_{\mathrm{sym}}(n,T=0)}{{\hslash }c}\right]}^{2}\\ & & \times \,{\left\{24n\left[1+\sqrt{1+\displaystyle \frac{{\pi }^{2}n}{288}{\left(\displaystyle \frac{{\hslash }c}{{E}_{\mathrm{sym}}(n,T=0)}\right)}^{3}}\right]\right\}}^{1/3}.\end{array}\end{eqnarray} \tag{ 22 }$$

For each of the EOS in our sample, we stitch together a complete cold EOS at β-equilibrium from the publicly available tables at fixed Y_p , by requiring that ${\mu }_{e}+{\mu }_{p}\,-{\mu }_{n}=0$, where μ_i is the chemical potential of each species. We then use the corresponding density-dependent proton fraction, Y_p,β , to fit for γ using Equations (14)–(20) and keeping S₀ and L fixed for each EOS. We perform the fits using a standard least-squares method and limit the density range to n ≥ 10⁻² fm⁻³. In principle, Equations (14)–(20) apply only to n–p–e matter, which will be uniform only above 0.5n_sat. However, in practice, we find a very small difference in the fits for γ whether we include densities above $0.5{n}_{\mathrm{sat}}=0.08$ fm⁻³ or whether we start the fits at a slightly lower but still astrophysically relevant cutoff of n = 10⁻² fm⁻³. We show the resulting fit values in Table 1.

We note that the range of EOS provided in Table 1 is intentionally broad. While the symmetry energy parameters of some of these EOS disagree with the combined set of experimental constraints (see Lattimer & Lim 2013 for a recent review) or are in disagreement with certain theoretical considerations, such as chiral effective field theory results for pure neutron matter (see, e.g., Krüger et al. 2013), they are all consistent with at least some experimental constraints on S₀ and L.

We find that γ spans roughly the range of experimentally allowed values, between 0.15 and 1.0, as expected. The only exception is SFHx, which has an extremely low value of L. This makes the result of the fit highly sensitive to the density range that is included. For consistency, we still constrain the densities to n ≥ 10⁻² fm⁻³ for the fit to this EOS; however, the inferred value for γ ranges from the reported value of −0.04 up to 0.18, depending on where the density cutoff is placed. Thus, the particular value for γ for SFHx should be taken with some caution.

We have here used Equation (20) to fit for γ from the β-equilibrium proton fractions of realistic EOS. We also wish to emphasize that Equation (20) can, of course, be used to calculate Y_p,β , given a choice of S₀, L, and γ. Once these three parameters are specified, Equations (21)–(22) can be used to calculate Y_p,β for any EOS. Consequently, all that is required of the cold EOS is knowledge of the run of pressure with density. This feature makes it possible to apply our model to piecewise polytropes or other families of parametric EOS that may not directly calculate Y_p,β .

We show an example of the performance of this model for E_sym(n, T = 0) in Figure 3 for the EOS NL3 (Lalazissis et al. 1997, 1999) and DD2 (Typel et al. 2010). We show these two EOS as representative samples, with NL3 representing the family of EOS with larger L values and DD2 representing the EOS with smaller symmetry energy slopes (see Table 1). The top panel of Figure 3 shows the zero-temperature pressure predicted by NL3 and DD2 at Y_p  = 0.1 as blue and orange diamonds, respectively. The colored lines show our model: starting with the corresponding EOS in β-equilibrium, adding the symmetry energy correction of Equations (14)–(19), and correcting for the leptons, all according to Equation 12(b). For these models, we take the values of S₀, L, and γ for each EOS from Table 1. We note that we are plotting pressures but we could have similarly shown the energy. We use Equation (13) to convert the equations of this section to pressures. For the complete set of pressure expressions, see Section 6 and Box II.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 3 shows the residuals between our model and the pressure predicted by each EOS at Y_p  = 0.1. We find that our model performs very well at densities above 0.5 n_sat, with errors ≲10%. At the highest densities, using our model compared to the full EOS introduces errors of only ∼1%. The residuals for the other EOS in our sample are comparably small.

For Y_p  = 0.3, we find that the residuals between our model and NL3 and DD2 are comparable to those shown in Figure 3. Therefore, we conclude that this model reasonably captures the Y_p -dependence of the cold EOS for a large range of L values.

We thus have an expression for the symmetry energy at zero-temperature that depends only on n, Y_p , S₀, L, and the narrowly constrained parameter γ. There are two possible routes to create a finite-temperature EOS with this framework. One possibility is to start from a cold, physically motivated EOS, which will provide predicted values for S₀, L, and Y_p,β . In this case, Equation (20) can be used to fit for γ. We have provided such fits for the EOS in our sample in Table 1. Alternatively, a cold, parametric EOS can be chosen, for which the underlying physics are not specified. In this case, a user can freely specify S₀, L, and γ, which will uniquely specify Y_p,β . For the EOS in our sample, we find that this approach is able to accurately extrapolate from β-equilibrium to arbitrary proton fraction, introducing errors of ≲10% for densities of interest (above 0.5 n_sat), and errors of ≲3% at high densities.

We now turn to the thermal energy, which was first defined in Equation 12(c) as

$$\begin{eqnarray*}\begin{array}{rcl}{E}_{\mathrm{th}}(n,{Y}_{p},T) & = & {E}_{\mathrm{nucl},\mathrm{th}}(n,{Y}_{p}=1/2,T)\\ & & +\,{E}_{\mathrm{sym},\mathrm{th}}(n,T)\times {\left(1-2{Y}_{p}\right)}^{2}\\ & & +\,{E}_{\mathrm{lepton},\mathrm{th}}(n,{Y}_{p},T).\end{array}\end{eqnarray*}$$

It is useful to further divide the thermal energy into density regimes, over which the matter displays distinct behaviors. At the lowest densities, the contribution from relativistic leptons and photons dominates. At intermediate densities, an ideal-fluid description suffices. However, at high densities, matter can remain partially degenerate, even at intermediate-to-high temperatures. In the high-density regime, some of the available energy goes into lifting the degeneracy of the particles rather than adding thermal support and, accordingly, the thermal pressure can dip well below the prediction for an ideal fluid. (See Figure 6 for the markedly different behaviors in thermal pressure across these three regimes.)

It is, therefore, convenient to write the thermal energy as

$$\begin{eqnarray}&&\begin{array}{l}{E}_{\mathrm{th}}(n,{Y}_{p},T)\,=\\ \left\{\begin{array}{ll}{E}_{\mathrm{rel}}(n,T), & n\lt {n}_{1}\\ {E}_{\mathrm{ideal}}(T), & {n}_{1}\lt n\lt {n}_{2}\\ {E}_{\mathrm{th},\deg .}(n,{Y}_{p}=1/2,T) & \\ \quad +\,{E}_{\mathrm{sym},\mathrm{th}}(n,T){\left(1-2{Y}_{p}\right)}^{2}, & n\gt {n}_{2}\end{array}\right.\end{array}\end{eqnarray} \tag{ 23 }$$

where the relativistic component,

$$\begin{eqnarray}&&{E}_{\mathrm{rel}}(n,T)=\displaystyle \frac{4\sigma }{c}\displaystyle \frac{{f}_{s}}{n}{T}^{4},\end{eqnarray} \tag{ 24 }$$

and the ideal component,

$$\begin{eqnarray}&&{E}_{\mathrm{ideal}}(T)=\displaystyle \frac{3}{2}{k}_{{\rm{B}}}T\end{eqnarray} \tag{ 25 }$$

are given as in Equations (2) and (3). Here, ${E}_{\mathrm{th},\deg .}(n,{Y}_{p}\,=1/2,T)$ is the degenerate thermal energy of symmetric matter, which we introduce below. We note that, because the ideal-fluid and relativistic terms do not depend on the proton fraction, the symmetry-energy correction is only relevant in the degenerate regime. Finally, we define the first transition density, n₁, as the density at which the relativistic and ideal-fluid energies are equal. The second transition density, n₂, is the density at which the ideal-fluid energy is equal to the degenerate thermal energy, for a given temperature and proton fraction.

This piecewise expression of the thermal energy is convenient for later calculations of the thermal pressure and sound speed. However, the discontinuities at the transition densities are artificial and will create problems in numerical simulations, potentially leading to undesired reflections of matter waves at density boundaries. Thus, whenever we actually implement the thermal energy or pressure, we use a smoothed version instead. This smoothed version is of the form

$$\begin{eqnarray}&&\begin{array}{l}{E}_{\mathrm{th}}(n,{Y}_{p},T)\approx {E}_{\mathrm{rel}}(n,T)\\ \,+\,\left[{E}_{\mathrm{ideal}}{\left(T\right)}^{-1}\right.{\left.+{E}_{\mathrm{th},\deg .}{\left(n,{Y}_{p},T\right)}^{-1}\right]}^{-1},\end{array}\end{eqnarray} \tag{ 26 }$$

where we have added the latter two terms inversely to ensure that the ideal term dominates at intermediate densities and the degenerate term dominates at the highest densities. The smoothed approximation is also more computationally efficient than the piecewise version because it does not require the calculation of transition densities, which will vary with the temperature and proton fraction.

To calculate the thermal energy in the degenerate regime, we consider the nucleons as a free Fermi gas. In that limit, the leading-order thermal energy of degenerate matter is given by

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{th},q}^{\deg }(n,{Y}_{q},T) & = & a({Y}_{q}n,{M}^{* })\left(\displaystyle \frac{{N}_{q}}{{N}_{p}+{N}_{n}}\right){T}^{2},\\ & = & a({Y}_{q}n,{M}^{* }){Y}_{q}{T}^{2}\end{array}\end{eqnarray} \tag{ 27 }$$

for a single-species system of particle q. For simplicity, we have introduced the level-density parameter a, which is defined as

$$\begin{eqnarray}&&a({n}_{q},{M}^{* })\equiv \displaystyle \frac{{\pi }^{2}{k}_{{\rm{B}}}^{2}}{2}\displaystyle \frac{\sqrt{{\left(3{\pi }^{2}{n}_{q}\right)}^{2/3}{\left({\hslash }c\right)}^{2}+{M}^{* }{\left({n}_{q}\right)}^{2}}}{{\left(3{\pi }^{2}{n}_{q}\right)}^{2/3}{\left({\hslash }c\right)}^{2}},\end{eqnarray} \tag{ 28 }$$

where ${M}^{* }({n}_{q})$ is the Dirac effective mass of the relevant species at a specific density. (For a complete derivation at next-to-leading order in temperature, see Constantinou et al. 2015.)

As an example, the thermal nuclear energy for symmetric matter would be

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{th},\mathrm{nucl}}^{\deg }(n,T) & = & \left[\displaystyle \frac{a({n}_{p},{M}_{p,\mathrm{SM}}^{* }){N}_{p}+a({n}_{n},{M}_{n,\mathrm{SM}}^{* }){N}_{n}}{{N}_{p}+{N}_{n}}\right]{T}^{2}\\ & = & a(0.5n,0.5{M}_{\mathrm{SM}}^{* }){T}^{2},\end{array}\end{eqnarray} \tag{ 29 }$$

where the subscript SM stands for symmetric matter and, in the second line, we have used the fact that ${n}_{n}={n}_{p}=0.5n$ in symmetric matter. We have further made the approximation that the effective masses of neutrons and protons are comparable in symmetric matter and that the average of these two effective masses gives the overall effective mass of symmetric matter; i.e., ${M}_{n,\mathrm{SM}}^{* }\approx {M}_{p,\mathrm{SM}}^{* }\approx 1/2{M}_{\mathrm{SM}}^{* }$.

By likewise defining the thermal energy per baryon for pure neutron matter, we can calculate the thermal contribution to the symmetry energy, as

$$\begin{eqnarray}&&\begin{array}{l}{E}_{\mathrm{sym},\mathrm{th}}(n,T)\,=\\ \left\{\begin{array}{ll}0, & n\lt {n}_{2}\\ \left[a(n,{M}_{\mathrm{PNM}}^{* })-a(0.5n,0.5{M}_{\mathrm{SM}}^{* })\right]{T}^{2}, & n\gt {n}_{2},\end{array}\right.\end{array}\end{eqnarray} \tag{ 30 }$$

where the low-density limit of ${E}_{\mathrm{sym},\mathrm{th}}$ arises from the fact that both pure neutron matter and symmetric matter behave identically as ideal or relativistic fluids at n < n₂.

In principle, this symmetry energy term extrapolates the thermal energy of symmetric nuclear matter to arbitrary proton fraction. However, we find that including this term has a negligible effect on the results. In particular, making the approximation ${E}_{\mathrm{th},\mathrm{nucl}}(n,{Y}_{p},T)\approx {E}_{\mathrm{th},\mathrm{nucl}}(n,{Y}_{p}=1/2,T)$ introduces an average error of ≲1% in the total pressure across the density range of interest. We thus neglect the thermal correction to the symmetry energy for the rest of this paper.

For leptons, the degenerate thermal pressure is even simpler. The effective mass of electrons is approximately constant, due to their small cross-sections of interaction. Hence, ${M}_{e}^{* }\approx {m}_{e}$. This allows us to write Equation (27) simply as

$$\begin{eqnarray}&&{E}_{\mathrm{th},{e}^{-}}^{\deg }(n,{Y}_{p},T)=a({Y}_{p}n,{m}_{e}){Y}_{p}{T}^{2},\end{eqnarray} \tag{ 31 }$$

where we have required that the electron fraction balance the proton fraction to satisfy the requirement of charge neutrality and we have used Equation (6) to substitute Y_p . We note that in the presence of a significant population of positrons, the proton fraction in Equation (31) should be replaced by the net lepton fraction.

With expressions for the degenerate and ideal-fluid thermal terms in hand, we can now write a complete version of Equation 12(c) for E_th as follows:

$$\begin{eqnarray}&&\begin{array}{l}{E}_{\mathrm{th}}(n,{Y}_{p},T)\,=\\ \left\{\begin{array}{ll}4\sigma {f}_{s}{T}^{4}/({cn}), & n\lt {n}_{1}\\ (3/2){k}_{{\rm{B}}}T, & {n}_{1}\lt n\lt {n}_{2}\\ \left[a(0.5n,0.5{M}_{\mathrm{SM}}^{* })+a({Y}_{p}n,{m}_{e}){Y}_{p}\right]{T}^{2}, & n\gt {n}_{2}\end{array}\right.\end{array}\end{eqnarray} \tag{ 32 }$$

where we have neglected the thermal contribution to the symmetry energy, as discussed above.

We thus have a complete expression for the thermal energy of matter as a function of only the density, temperature, proton fraction, and the effective mass of the nucleons in symmetric matter.

5.1.  M*-approximation

A full calculation of E_th using Equation (32) requires knowledge of the Dirac effective masses in symmetric matter, and hence the scalar meson interactions and particle potentials of a particular EOS. We instead choose to express the Dirac effective mass with a physically motivated yet computationally simple approximation. At low densities, the effective mass must approach the dominant nucleon mass, while at higher densities, M* must decrease as particle interactions become important. We represent this behavior by introducing a power-law expression,

$$\begin{eqnarray}&&{M}^{* }({n}_{q})={\left\{{\left({{mc}}^{2}\right)}^{-b}+{\left[{{mc}}^{2}{\left(\displaystyle \frac{{n}_{q}}{{n}_{0}}\right)}^{-\alpha }\right]}^{-b}\right\}}^{-1/b},\end{eqnarray} \tag{ 33 }$$

where m is the nucleon mass (which we take to be the neutron mass, ${{mc}}^{2}=939.57$ MeV) ⁵ and n₀ is the transition density above which M* starts to decrease. The exponent b determines the sharpness of the transition and α specifies the power-law slope at high densities. We find that b = 2 works well to represent the curvature connecting the low- and high-density regimes, and we thus fix it to this value in the following analysis, leaving just two free parameters to describe the effective mass, ${M}^{* }={M}^{* }({n}_{0},\alpha )$.

We fit the effective masses together at ${k}_{{\rm{B}}}T=1,10,$ and 47.9 MeV for nine of the EOS in our sample, using a standard least-squares method across the entire density range provided. We exclude the LS and SLY4-RG models here because the effective masses for these EOS are not currently published (but see Section 5.3 for a separate comparison with these models). The results of these fits for symmetric matter are given in Table 2. For completeness, we also include the fits for pure neutron matter in Table 2, which can be used to calculate ${E}_{\mathrm{sym},\mathrm{th}}(n,T)$ in Equation (30).

Table 2. Parameters Characterizing M*, Fit Together at k_B T = 1, 10, and 47.9 MeV, for either Pure Neutron Matter (PNM) or Symmetric Matter (SM)

	PNM (Yp  = 0.01)		SM (Yp  = 0.5)
EOS	n0 (fm−3)	α	n0 (fm−3)	α
TM1	0.11	0.73	0.12	0.86
TMA	0.11	0.65	0.13	0.77
NL3	0.10	0.90	0.11	1.08
FSUGold	0.10	0.61	0.11	0.72
IUFSU	0.11	0.72	0.12	0.85
DD2	0.08	0.68	0.10	0.84
STOS	0.11	0.76	0.12	0.90
SFHo	0.21	0.82	0.22	0.89
SFHx	0.16	0.77	0.17	0.88
Range	0.08–0.21	0.61–0.90	0.10–0.22	0.72–1.08
Mean	0.12	0.74	0.13	0.87

Note. We fix b = 2 and m = m_n in all fits.

Download table as:  ASCIITypeset image

We show the performance of the fit for NL3 in Figure 4. In this fit, we use the NL3 tables calculated at k_B T = 1, 10, and 47.9 MeV (shown in purple, orange, and blue, respectively) with a proton fraction of Y_p  = 0.01, to emulate pure neutron matter. We show our approximation for M* as a black solid line. We find that the M*-approximation accurately captures the behavior predicted by the full EOS, with fit parameters n₀ = 0.10 fm⁻³ and α = 0.90. Figure 5 shows the M* predictions and our approximation thereof for Y_p  = 0.5. The M*-approximation is similarly able to capture the behavior of symmetric matter, with slightly adjusted parameters n₀ = 0.11 fm⁻³ and α = 1.08.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a brief aside, we note a discontinuity in the first derivative of M* at approximately half the nuclear saturation density for large Y_p and low temperatures (seen most clearly in the purple stars in Figure 5, at n_sat/2 ≈ 0.08 fm⁻³). This discontinuity is an artifact of the treatment of the first-order phase transition to uniform nuclear matter at these densities in the original EOS calculations.

There is an easily understood origin of this artifact. Lattimer & Swesty (1991), Shen et al. (1998), and Hempel & Schaffner-Bielich (2010) all use a Maxwell construction to calculate the phase transition at approximately half the nuclear saturation density. At low proton fractions, where matter is approximately made up of a single species, the Maxwell construction works well to represent the phase transition. However, the Maxwell construction is invalid for multi-component species. When a system has more than one significant component, the Gibbs construction must instead be used (Glendenning 1992, 2000). Because all of the EOS that are included in this section use the Maxwell construction, they all suffer from artifacts due to this choice at roughly half the saturation density, where the transition to uniform nuclear matter occurs.

Correcting these artifacts would require us to re-calculate all of the EOS with a different formalism, which is beyond the scope of this paper. However, we note that at high temperatures (k_B T ≳ 15 MeV), the non-uniform phase of matter disappears (see the discussion around Figure 5 in Shen et al. 1998). Thus, we can avoid the issue altogether by performing our fit to M* at only the highest temperatures, when Y_p is large. In practice, we find that whether we fit only the k_B T = 47.9 MeV curve for M* or if we fit the curves for all the temperatures together, the difference in the resulting parameters is small. Therefore, we choose to perform the fits to three temperatures (k_B T = 1, 10 and 47.9 MeV) together and use the same method for both low and high proton fractions.

Returning to our discussion of the M* model, we note that the errors introduced by using our M*-approximation are comparable to those shown in Figures 4 and 5 for the full set of nine EOS in this section. We thus conclude that our M*-approximation reasonably captures the density-dependence of the Dirac effective mass, while greatly simplifying subsequent calculations.

Moreover, we find that the range of inferred fit parameters is relatively narrow. In particular, for a wide range of temperatures and EOS, we find that the transition density lies in the range n₀ ∈ (0.08, 0.22) fm⁻³, with an average value of ∼0.13 fm⁻³ for both pure neutron matter and symmetric matter. The power-law index characterizing the decay of M* is similarly well constrained, with α ∈ (0.61–0.90), with an average value of 0.74 for pure neutron matter; and α ∈ (0.72–1.08), with a slightly higher average value of 0.87 for symmetric matter. We find only a weak dependence of n₀ and α on the temperature, thus suggesting that these parameters could be treated as constants for use in numerical simulations.

5.2. Performance of the M*-approximation of Thermal Effects at Fixed Y_p

We now turn to a comparison between the M*-approximation of the thermal effects and the nine EOS listed in Table 2. As in Section 4, we make the comparison in terms of the pressure rather than the energy and we use Equation (13) to convert between the two. The expressions for P_th(n, Y_p , T) are given in Box II in Section 6. In particular, all of the results shown here use the smoothed approximation of the thermal pressure, as defined in Equation (39).

In order to focus specifically on the thermal pressure, we calculate the thermal contribution to the pressure from each realistic EOS in our sample by subtracting the cold component at the same Y_p .

In general, we find excellent agreement between the M*-approximation and the thermal pressures calculated from the full EOS. We show an example in Figure 6 for NL3. We find that our approximation of P_th closely recreates the full calculation for NL3 for nearly all of the densities and temperatures explored here. For comparison, we also include in Figure 6 the hybrid approximation with Γ_th = 1.67 as dashed lines. ⁶ The full thermal pressure only agrees with the hybrid approximation at intermediate densities. At the lowest densities, this value of Γ_th overestimates the contribution from relativistic species. At higher densities that are relevant for forming and merging neutron stars, particle interactions become important and the ideal-fluid approximation grossly overestimates the thermal pressure, remaining several orders of magnitude above the true thermal pressure.

Zoom In Zoom Out Reset image size

In order to gain an intuitive understanding of the behavior of P_th, we also explore an extreme range of the M* parameters. Specifically, in Figure 7, we zoom in on P_th at k_B T = 10 MeV and Y_p  = 0.5 and show the effect of varying the parameters n₀ and α for symmetric matter. We intentionally take extreme values for the parameters, well beyond the ranges found in Table 2, to emphasize that the variations between more realistic parameter choices will be small. Even for these unreasonable choices of values for n₀ and α, we find that P_th approximates the full thermal pressure reasonably well and is in all cases better than the ideal-fluid approximation. Analyzing the specific dependencies more closely, we see in Figure 7 that the parameter n₀ controls the density at which the rise in the thermal pressure starts to slow. This corresponds to the density at which particle interactions become significant and where degenerate thermal effects can no longer be ignored. The parameter α, which controls the power-law slope of M*, directly controls the height of the dip in P_th. This makes intuitive sense: if particle interactions are stronger, then M* decreases more rapidly, α will be larger, and the thermal pressure will deviate even more drastically from the ideal-fluid approximation as part of the free energy is taken up by these interactions.

Zoom In Zoom Out Reset image size

Finally, we compare the M*-approximation of the thermal pressure against the full sample of EOS listed in Table 2. We show the corresponding residuals at three temperatures in Figure 8 and find that the residuals are typically ≲30% at densities above 0.5 n_sat. For comparison, Figure 8 also shows a sample set of residuals between the full thermal pressure from NL3 and the hybrid approximation (Γ_th = 1.67) as a black dashed line. We find that the M*-approximation produces residuals that are up to three orders of magnitude smaller than the ideal-fluid approximation used in hybrid EOS, with only two additional parameters that are easy to specify.

Zoom In Zoom Out Reset image size

5.3.  M*-approximation for Non-RMF Models

We have so far only calculated the thermal pressures using the sub-sample of EOS for which there exist published tables of the effective masses. While this allowed us to directly test the performance of the M*-approximation, this set of EOS happens to also be calculated exclusively with RMF models. In this section, we compare the M*-approximation to the LS and SLY4-RG models, which are calculated using non-relativistic Skyrme energy functionals (see Section 2). We also include here Zhang & Prakash's (2016) two-loop exchange model, which is an extension of mean field theory. We note that the pressures of the Zhang & Prakash (2016) EOS are reported only at Y_p  = 0 and 0.5, which is why this EOS is not included in our full sample. As a result of these and other limitations in the publicly available values for this EOS, all of the comparisons in this section are made at Y_p  = 0.5 and T = 20 MeV. We also fix n₀ and α to the mean values for symmetric matter from Table 2 for all three EOS.

Figure 9 shows the residuals between the M*-approximation of the thermal pressure and the true EOS for these three models. For comparison, this figure also shows the corresponding residuals between the hybrid approximation and the true EOS (dashed lines). In general, we find that the M*-approximation of the thermal pressure results in larger residuals for these EOS compared to the RMF models but that it still offers a significant improvement over the hybrid approximation at densities above ∼n_sat.

Zoom In Zoom Out Reset image size

We also compared the residuals at T = 50 MeV and found that the M*-approximation performed comparably to the hybrid approximation at this temperature. In fact, for densities between n_sat and 0.7 fm⁻³, the hybrid approximation produces slightly smaller residuals in the thermal pressure for these non-RMF models. In this regime, the hybrid approximation tends to overestimate the thermal pressure for these models, while the M*-approximation tends to under-estimate P_th by a similar degree. However, even in this case, the M*-approximation still offers an appreciable improvement over the hybrid approximation at the highest densities, above ∼0.7 fm⁻³.

We now summarize the equations and approximations that we have developed so far to represent the total energy per particle in Box I.

Box I.  Total Energy Expressions for Finite-temperature Dense Gas.

The energy per particle of n–p–e matter is given by

$$\begin{eqnarray*}&&\begin{array}{l}E(n,{Y}_{p},T)=\left(\mathrm{Cold}\ \mathrm{EOS}\ \mathrm{in}\ \beta \mbox{-} \mathrm{equilibrium}\right)+\,3K\left({Y}_{p}^{4/3}-{Y}_{p,\beta }^{4/3}\right){n}^{1/3}\\ \,\,\,+\,{E}_{\mathrm{sym}}(n,T=0)\left[{\left(1-2{Y}_{p}\right)}^{2}-{\left(1-2{Y}_{p,\beta }\right)}^{2}\right]+\,\left\{\begin{array}{ll}4\sigma {f}_{s}{T}^{4}/({cn}), & n\lt {n}_{1}\\ (3/2){k}_{{\rm{B}}}T, & {n}_{1}\lt n\lt {n}_{2}\\ \left[a(0.5n,0.5{M}_{\mathrm{SM}}^{* })+a({Y}_{p}n,{m}_{e}){Y}_{p}\right]{T}^{2}, & \,n\gt {n}_{2},\end{array}\right.\end{array}\end{eqnarray*}$$

where the symmetry energy is approximated as

$$\begin{eqnarray*}&&{E}_{\mathrm{sym}}(n,T=0)=\eta {E}_{\mathrm{sym}}^{\mathrm{kin}}(n)+\left[{S}_{0}-\eta {E}_{\mathrm{sym}}^{\mathrm{kin}}({n}_{\mathrm{sat}})\right]{\left(\displaystyle \frac{n}{{n}_{\mathrm{sat}}}\right)}^{\gamma },\end{eqnarray*}$$

$$\begin{eqnarray*}&&{E}_{\mathrm{sym}}^{\mathrm{kin}}(n)=\displaystyle \frac{3}{5}\left(1-{2}^{1/3}\right){E}_{{\rm{f}}}(n),\end{eqnarray*}$$

$$\begin{eqnarray*}&&\eta =\displaystyle \frac{5}{9}\left[\displaystyle \frac{L-3{S}_{0}\gamma }{\left(1-{2}^{1/3}\right)\left(2/3-\gamma \right){E}_{{\rm{f}}}({n}_{\mathrm{sat}})}\right],\end{eqnarray*}$$

and the terms of the M*-approximation are given by

$$\begin{eqnarray*}&&\begin{array}{l}a({n}_{q},{M}_{q}^{* })\equiv \displaystyle \frac{{\pi }^{2}{k}_{{\rm{B}}}^{2}}{2}\displaystyle \frac{\sqrt{{M}_{q}^{* }{\left({n}_{q}\right)}^{2}+{\left(3{\pi }^{2}{n}_{q}\right)}^{2/3}{\left({\hslash }c\right)}^{2}}}{{\left(3{\pi }^{2}{n}_{q}\right)}^{2/3}{\left({\hslash }c\right)}^{2}}\,\,\,\mathrm{and}\,\,\,{M}^{* }({n}_{q})={\left\{{\left({{mc}}^{2}\right)}^{-b}+{\left[{{mc}}^{2}{\left(\displaystyle \frac{{n}_{q}}{{n}_{0}}\right)}^{-\alpha }\right]}^{-b}\right\}}^{-1/b}.\end{array}\end{eqnarray*}$$

• The parameters S₀, L, and γ ∈ (0.2–1.2) are freely specified, which will uniquely specify Y_p,β .
• Alternatively, S₀, L, and Y_p,β may be specified and the proton fraction may be fit for γ. We provide fits to γ for eleven EOS in Table 1.
• We find that for ${M}_{\mathrm{SM}}^{* }$, n₀ ∼ 0.13 fm⁻³ and α ∼ 0.9 provide reasonable fits to most EOS.

Using the expressions for the energy from Box I, we can derive the pressure via the standard thermodynamic relations of Equation (13), where the derivatives are evaluated at constant Y_p , Y_p,β , and S. The total entropy of the relativistic, ideal-fluid, and degenerate terms is given by

$$\begin{eqnarray}S(n,{N}_{p},{N}_{n},{N}_{e},T)=\left\{\begin{array}{ll}{S}_{\mathrm{rel}}, & n\lt {n}_{1}\\ {S}_{\mathrm{ideal}}, & {n}_{1}\lt n\lt {n}_{2}\\ {S}_{\deg }, & n\gt {n}_{2},\end{array}\right.\end{eqnarray} \tag{ 34 }$$

where n₁ and n₂ are the thermal energy transition densities, as defined in Section 5.

The entropy of a gas of relativistic leptons and photons is given by

$$\begin{eqnarray}&&{S}_{\mathrm{rel}}=\displaystyle \frac{16\sigma {f}_{s}}{3c}\left(\displaystyle \frac{{N}_{p}+{N}_{n}}{n}\right){T}^{3}.\end{eqnarray} \tag{ 35 }$$

The entropy of a monatomic ideal fluid is given by the Sackur–Tetrode equation,

$$\begin{eqnarray}&&\begin{array}{l}{S}_{\mathrm{ideal}}=\left({N}_{p}+{N}_{n}+{N}_{e}\right){k}_{{\rm{B}}}\\ \times \,\left\{\mathrm{ln}\left[\left(\displaystyle \frac{{N}_{p}+{N}_{n}}{{N}_{p}+{N}_{n}+{N}_{e}}\right){n}^{-1}{\left(\displaystyle \frac{{{mk}}_{{\rm{B}}}T}{2\pi {{\hslash }}^{2}}\right)}^{3/2}\right]+\displaystyle \frac{5}{2}\right\}.\end{array}\end{eqnarray} \tag{ 36 }$$

Finally, the entropy of a degenerate Fermi gas in our framework is given by

$$\begin{eqnarray}&&{S}_{q}=2{a}_{q}{N}_{q}T\end{eqnarray} \tag{ 37 }$$

for a particle q, so that the total entropy for the degenerate terms is

$$\begin{eqnarray}&&{S}_{\deg }=2\left\{a(0.5n,0.5{M}_{\mathrm{SM}}^{* })[{N}_{n}+{N}_{p}]+a({Y}_{p}n,{m}_{e}){N}_{e}\right\}T.\end{eqnarray} \tag{ 38 }$$

We summarize the resulting pressure equations in Box II.

Box II.  Pressure Expressions for Finite-temperature Dense Gas.

The pressure of n–p–e matter is given by

$$\begin{eqnarray*}&&\begin{array}{l}P(n,{Y}_{p},T)=\left(\mathrm{Cold}\ \mathrm{EOS}\ \mathrm{in}\ \beta \mbox{-} \mathrm{equilibrium}\right)+\,K\left({Y}_{p}^{4/3}-{Y}_{p,\beta }^{4/3}\right){n}^{4/3}\\ \,+\,{P}_{\mathrm{sym}}(n,T=0)\left[{\left(1-2{Y}_{p}\right)}^{2}-{\left(1-2{Y}_{p,\beta }\right)}^{2}\right]+\,\left\{\begin{array}{lr}4\sigma {f}_{s}{T}^{4}/(3c), & n\lt {n}_{1}\\ {{nk}}_{{\rm{B}}}T, & {n}_{1}\lt n\lt {n}_{2}\\ -\,\left[\displaystyle \frac{\partial a(0.5n,0.5{M}_{\mathrm{SM}}^{* })}{\partial n}+\displaystyle \frac{\partial a({Y}_{p}n,{m}_{e})}{\partial n}{Y}_{p}\right]{n}^{2}{T}^{2},\, & \,n\gt {n}_{2},\end{array}\right.\end{array}\end{eqnarray*}$$

where n₁ and n₂ are the thermal energy transition densities for a particular temperature and proton fraction.

The symmetry pressure, corresponding to our model of the symmetry energy, is

$$\begin{eqnarray*}&&{P}_{\mathrm{sym}}(n,T=0)=\displaystyle \frac{2\eta }{3}{{nE}}_{\mathrm{sym}}^{\mathrm{kin}}(n)+\left[{S}_{0}-\eta {E}_{\mathrm{sym}}^{\mathrm{kin}}({n}_{\mathrm{sat}})\right]{\left(\displaystyle \frac{n}{{n}_{\mathrm{sat}}}\right)}^{\gamma }\gamma n.\end{eqnarray*}$$

The full analytic expression for Y_p,β is given in Equation (20) and is derived in Appendix A.

The M*-approximation derivatives are given by

$$\begin{eqnarray*}&&\begin{array}{l}{\left.\displaystyle \frac{\partial a({n}_{q},{M}^{* })}{\partial n}\right|}_{{Y}_{q}}=-\displaystyle \frac{2a({n}_{q},{M}^{* })}{3n}\left\{1-\displaystyle \frac{1}{2}\left[\displaystyle \frac{{M}^{* }{\left({n}_{q}\right)}^{2}}{{M}^{* }{\left({n}_{q}\right)}^{2}+{\left(3{\pi }^{2}{n}_{q}\right)}^{2/3}{\left({\hslash }c\right)}^{2}}\right]\right.\left.\left.{\left(\displaystyle \frac{{\left(3{\pi }^{2}{n}_{q}\right)}^{2/3}{\left({\hslash }c\right)}^{2}}{{M}^{* }{\left({n}_{q}\right)}^{2}}+3\displaystyle \frac{\partial \mathrm{ln}[{M}^{* }({n}_{q})]}{\partial \mathrm{ln}n}\right|}_{{Y}_{q}}\right)\right\},\end{array}\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{\left.\displaystyle \frac{\partial \mathrm{ln}[{M}^{* }({n}_{q})]}{\partial \mathrm{ln}n}\right|}_{{Y}_{q}}=-\alpha \left[1-{\left(\displaystyle \frac{{M}^{* }({n}_{q})}{{Y}_{q}{{mc}}^{2}}\right)}^{2}\right],\end{eqnarray*}$$

where, for symmetric matter, we replace ${M}^{* }({n}_{q})\to 0.5{M}_{\mathrm{SM}}^{* }(0.5n)$ and for the electrons, ${M}^{* }({n}_{q})\to {m}_{e}$.

• As in Box I, there are five free parameters: ${S}_{0},L,\gamma ,{n}_{0},\alpha$.
• A user may freely specify S₀, L, and Y_p,β and fit for γ. Alternatively, a user may specify S₀, L, and γ, which will uniquely specify Y_p,β .
• We provide fits for γ, n₀, and α for the EOS in our sample in Tables 1 and 2.

The piecewise definitions of the thermal energy and pressure are mathematically convenient, but the sharp transitions are themselves unphysical, as discussed in Section 5. Therefore, we instead implement the thermal pressure using a smoothed approximation of the form

$$\begin{eqnarray}&&{P}_{\mathrm{th}}(n,{Y}_{p},T)\approx {P}_{\mathrm{rel}}+{\left({P}_{\mathrm{ideal}}^{-1}+{P}_{\deg }^{-1}\right)}^{-1}.\end{eqnarray} \tag{ 39 }$$

This smoothed approximation of the thermal pressure is used for the figures throughout this paper. We note that we use this separate smoothing for both the thermal pressure and the thermal energy (as in Equation (26)) to keep the problem tractable. However, this is not mathematically exact because, formally, the energy is the proper thermodynamic function and the pressure should, ideally, be derived from the smoothed energy. Nevertheless, the errors introduced by the separate smoothing approximations will be limited to the regions close to the transition points. Physically, the mismatch between the approximate thermal energy and pressure will correspond to a small error in the sound speed in these regions, which we neglect for the present purposes.

Finally, we note that our model allows significant freedom in creating a new finite-temperature EOS. We have provided a set of parameters that correspond to physically motivated EOS. If one wishes to vary these parameters significantly, it will be useful to check that the resulting EOS is still physical. One requirement of a realistic EOS is that the sound speed remain subluminal at all densities and temperatures of interest. For this reason, we include in Appendix B a calculation of the sound speed for astrophysical merger scenarios.

In Section 4, we found that our model is able to extrapolate from β-equilibrium to an arbitrary proton fraction with resulting errors of ≲10% at densities above 0.5 n_sat. Similarly, in Section 5.2, we showed that the M*-approximation is able to reproduce the thermal pressure of realistic EOS, at fixed Y_p , to within ∼30% for a variety of EOS based on RMF theory. In this section, we quantify the performance of our complete model: starting with a cold EOS in β-equilibrium, and extrapolating to arbitrary temperature and proton fraction.

Figure 10 shows an example of a complete model for NL3 and DD2 at three different temperatures. For our approximation, we start with the relevant cold EOS in β-equilibrium and add the corrections outlined in Box II to extrapolate the pressure to Y_p  = 0.1 and the three indicated temperatures. We take the values for n₀, α, and γ listed in Tables 1 and 2 for each EOS. We show the results as the solid lines in Figure 10, while the predictions of the full EOS are shown as the diamonds. We find close agreement between our approximation and the full pressures predicted by NL3 and DD2, especially at densities above ∼0.5 n_sat.

Zoom In Zoom Out Reset image size

Figure 11 shows the corresponding residuals between our approximation and the full EOS for NL3 and DD2, as well as the rest of our sample of EOS. For each EOS in this figure, we use the values for n₀, α and γ listed in Tables 1 and 2, where possible. For LS and SLY4-RG, for which we do not have fit values for n₀ and α, we use the average parameter values for symmetric matter in Table 2. We find that our approximation works comparably well to recreate any of the EOS in our sample. Moreover, we find that for n ≳ n_sat, the residuals are ≲20% at all three temperatures.

Zoom In Zoom Out Reset image size

For all the EOS in our sample, the error introduced by our model increases in the vicinity of ∼0.5n_sat. This is a result of the break-down in the E_sym approximation at low densities. Our derivation of E_sym in Section 4 assumed uniform n–p–e matter, but at densities below ∼0.5 n_sat, the matter becomes inhomogeneous. Nevertheless, with the exception of LS, the errors at these densities are still typically ≲50%.

We have thus verified that our model is able to recreate realistic EOS at relevant densities, with a simple set of parameters. The implications of this result are two-fold. First, this approximation can be used in lieu of more complicated calculations to analytically represent the EOS that are commonly used in the literature with reasonable accuracy. Second, it implies that our approximation can be reliably used to create new finite-temperature EOS for n–p–e matter that probe different physics through the choice of n₀, α, γ, S₀, and L. Our model allows further freedom to create a new finite-temperature EOS through the choice of the cold, β-equilibrium EOS. We thus find that this model can span a broad range of possible physics, with parameters that are directly tied to the underlying physics and which can be integrated with minimal computational cost to a large array of numerical calculations.

In this paper, we have developed a general framework to calculate the pressure of neutron-star matter at arbitrary proton fraction and finite temperature. Our model is designed so that the corrections that we have developed here can be added to any cold n–p–e EOS in neutrinoless β-equilibrium. The model is based on a set of five physically motivated parameters: S₀, L, γ, n₀, and α. The first three, S₀, L, and γ, characterize the symmetry energy and can be chosen to match a particular EOS or set of priors from laboratory experiments. The parameters n₀ and α are introduced through our M*-approximation, where n₀ represents the density at which particle interactions become important and α characterizes the strength of those interactions. We find that the effective masses of nine realistic EOS can be well characterized by our M*-approximation with a relatively narrow range of these parameters, with average values of n₀ ∼ 0.13 fm⁻³ and α ∼ 0.9.

The complete model is able to extrapolate from cold matter in β-equilibrium to arbitrary proton fraction and temperature. We find that our model is able to recreate a sample of eleven realistic EOS with resulting errors of ≲20% at a variety of temperatures and proton fractions, above n_sat. In particular, by including the effects of degenerate matter, our M*-approximation reproduces the thermal pressure of realistic EOS with residuals that are several orders of magnitude smaller than the hybrid EOS that are commonly used in the literature.

In addition to providing a 1–3 orders-of-magnitude improvement over the ideal-fluid approximation of the thermal pressure, this model also includes the effects of changing the proton fraction, which is particularly relevant in simulating the formation and cooling of proto-neutron stars.

The complete model can thus be used to accurately recreate the realistic EOS that are currently in use in the literature with a set of simple, analytic functions. Furthermore, the model can be used to calculate new finite-temperature EOS that span a wide range of underlying physics, following one of two possible paths. One possibility is to choose a physically motivated cold EOS, which will provide predictions for the β-equilibrium proton fraction and symmetry energy parameters. These can then be used to fit for the free parameter γ and then extrapolate to an arbitrary proton fraction. Alternatively, one can use a cold, parametric EOS that does not specify the microphysics. In this case, there is freedom to choose the symmetry energy parameters to probe entirely new physics. In either case, one can freely choose the interaction parameters to control the relative importance of thermal effects. All together, these possibilities will allow a new and wide range of physics to be robustly probed in studies of dynamical neutron star phenomena.

We thank Vasileios Paschalidis for useful discussions and comments on this work. C.R. is supported by the NSF Graduate Research Fellowship Program Grant DGE-1143953. F.O. and D.P. acknowledge support from NASA grant NNX16AC56G.

In this appendix, we derive the relationship shown in Equation (20), which asserts that the cold β-equilibrium proton fraction is uniquely specified by the symmetry energy.

At zero-temperature for n–p–e matter, the total energy per baryon is given by

$$\begin{eqnarray}&&{E}_{\mathrm{tot}}(n,{Y}_{p})={E}_{n}(n,{Y}_{n})+{E}_{p}(n,{Y}_{p})+{E}_{e}(n,{Y}_{e}),\end{eqnarray} \tag{ 40 }$$

where E_n is the energy per baryon of neutrons, E_p is the energy per baryon of protons, and E_e the energy per baryon of electrons. Here, we also introduce the neutron fraction, Y_n  ≡ (1 − Y_p ), and the electron fraction, Y_e  = Y_p , where the latter equality holds in a charge-neutral system.

To find the minimum of the total energy, we differentiate with respect to Y_p and get

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{\mathrm{tot}}(n,{Y}_{p})}{\partial {Y}_{p}}=\displaystyle \frac{\partial {E}_{n}}{\partial {Y}_{n}}\displaystyle \frac{\partial {Y}_{n}}{\partial {Y}_{p}}+\displaystyle \frac{\partial {E}_{p}}{\partial {Y}_{p}}+\displaystyle \frac{\partial {E}_{e}}{\partial {Y}_{e}}\displaystyle \frac{\partial {Y}_{e}}{\partial {Y}_{p}},\end{eqnarray} \tag{ 41 }$$

where all the partial derivatives here and throughout this appendix are evaluated at constant entropy and baryon density and we have suppressed the notation for clarity. By substituting in the chemical potential of a species i, given by ${\mu }_{i}\equiv \partial {E}_{i}/\partial {Y}_{i}{| }_{S,n}$, Equation (41) simplifies to

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{\mathrm{tot}}(n,{Y}_{p})}{\partial {Y}_{p}}=-{\mu }_{n}+{\mu }_{p}+{\mu }_{e},\end{eqnarray} \tag{ 42 }$$

which is zero in β-equilibrium.

Alternatively, we can write the total energy as an expansion about nuclear symmetric matter with electrons added; i.e.,

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{tot}}(n,{Y}_{p}) & = & {E}_{\mathrm{nucl}}(n,1/2)\\ & & +\,{E}_{\mathrm{sym}}(n){\left(1-2{Y}_{p}\right)}^{2}+{E}_{e}(n,{Y}_{e}).\end{array}\end{eqnarray} \tag{ 43 }$$

This results in

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{\mathrm{tot}}(n,{Y}_{p})}{\partial {Y}_{p}}=-4(1-2{Y}_{p}){E}_{\mathrm{sym}}(n)+\displaystyle \frac{\partial {E}_{e}}{\partial {Y}_{e}}\displaystyle \frac{\partial {Y}_{e}}{\partial {Y}_{p}}.\end{eqnarray} \tag{ 44 }$$

By charge neutrality and the definition of the chemical potential, the second term is simply ${\mu }_{e}$. Combining Equations (42) and (44) in β-equilibrium gives

$$\begin{eqnarray}&&{\mu }_{e}=4(1-2{Y}_{p,\beta }){E}_{\mathrm{sym}}(n).\end{eqnarray} \tag{ 45 }$$

For relativistic electrons,

$$\begin{eqnarray}&&{\mu }_{e}=\sqrt{{p}_{{\rm{f}}}^{2}{c}^{2}+{m}_{e}^{2}{c}^{4}}\approx {p}_{{\rm{f}}}c,\end{eqnarray} \tag{ 46 }$$

where ${p}_{{\rm{f}}}c={(3{\pi }^{2}{Y}_{e}n)}^{1/3}{\hslash }c$ is the Fermi momentum of the electrons. Combining this expression for μ with Equation (45) yields

$$\begin{eqnarray}&&{\left(3{\pi }^{2}{Y}_{p,\beta }n\right)}^{1/3}{\hslash }c=4(1-2{Y}_{p,\beta }){E}_{\mathrm{sym}}(n),\end{eqnarray} \tag{ 47 }$$

or, rearranged to match Equation (20),

$$\begin{eqnarray}&&\displaystyle \frac{{Y}_{p,\beta }}{{\left(1-2{Y}_{p,\beta }\right)}^{3}}=\displaystyle \frac{64{E}_{\mathrm{sym}}{\left(n\right)}^{3}}{3{\pi }^{2}n{\left({\hslash }c\right)}^{3}}.\end{eqnarray} \tag{ 48 }$$

Solved for Y_p,β , this gives

$$\begin{eqnarray}&&{Y}_{p,\beta }=\displaystyle \frac{1}{2}+\displaystyle \frac{{\left(2{\pi }^{2}\right)}^{1/3}}{32}\displaystyle \frac{n}{\xi }\left\{{\left(2{\pi }^{2}\right)}^{1/3}-\displaystyle \frac{{\xi }^{2}}{n}{\left[\displaystyle \frac{{\hslash }c}{{E}_{\mathrm{sym}}(n,T=0)}\right]}^{3}\right\},\end{eqnarray} \tag{ 49 }$$

as in Equation (21), and where ξ is defined as

$$\begin{eqnarray}\begin{array}{rcl}\xi & \equiv & {\left[\displaystyle \frac{{E}_{\mathrm{sym}}(n,T=0)}{{\hslash }c}\right]}^{2}\\ & & \times \,{\left\{24n\left[1+\sqrt{1+\displaystyle \frac{{\pi }^{2}n}{288}{\left(\displaystyle \frac{{\hslash }c}{{E}_{\mathrm{sym}}(n,T=0)}\right)}^{3}}\right]\right\}}^{1/3},\end{array}\end{eqnarray} \tag{ 50 }$$

as in Equation (22).

Thus, if the form of E_sym(n) is known, then this will uniquely specify Y_p,β . Alternatively, if the β-equilibrium proton fraction is known from the cold EOS, then it can be used to fit for the parameters of the particular model of E_sym. In the context of this paper, specifying Y_p,β , S₀, and L can be used to fit for the parameter γ; or, specifying S₀, L, and γ can be used to calculate Y_p,β . The latter option is particularly useful because it allows our framework to be applied to parametric EOS that may not calculate Y_p,β directly.

Finally, we provide the derivative of ${Y}_{p,\beta }$, which is required to calculate the sound speed in Appendix B. The derivative at constant entropy is given by

$$\begin{eqnarray}\begin{array}{rcl}{\left.\displaystyle \frac{\partial {Y}_{p,\beta }}{\partial n}\right|}_{S} & = & \displaystyle \frac{1}{16}{\left(\displaystyle \frac{\pi }{2}\right)}^{2/3}\left\{{\left(2{\pi }^{2}\right)}^{1/3}\left(\displaystyle \frac{1-n\phi }{\xi }\right)\right.\\ & & +\,\left.\xi {\left[\displaystyle \frac{{\hslash }c}{{E}_{\mathrm{sym}}(n,T=0)}\right]}^{3}\left(3x-\phi \right)\right\}\end{array}\end{eqnarray} \tag{ 51 }$$

where for simplicity we have introduced the quantities

$$\begin{eqnarray}\begin{array}{rcl}\phi & = & \displaystyle \frac{1}{2n}\left(3{xn}+1\right)\\ & & +\,\displaystyle \frac{1}{6n}{\left\{1+\displaystyle \frac{{\pi }^{2}n}{288}{\left[\displaystyle \frac{{\hslash }c}{{E}_{\mathrm{sym}}(n,T=0)}\right]}^{3}\right\}}^{-1/2}(3{xn}-1)\end{array}\end{eqnarray} \tag{ 52 }$$

and

$$\begin{eqnarray}&&x=\displaystyle \frac{1}{{n}^{2}}\displaystyle \frac{{P}_{\mathrm{sym}}(n,T=0)}{{E}_{\mathrm{sym}}(n,T=0)}.\end{eqnarray} \tag{ 53 }$$

In this paper, we have provided the complete set of expressions necessary to extend any cold EOS to non-equilibrium conditions and arbitrary temperature. These expressions can be used to create a new finite-temperature EOS, by varying either the cold, underlying EOS or any of the five parameters of our model. In creating a new EOS, it is useful to always to check that the choice of parameters results in a model that is causal at all densities and temperatures of interest. To that end, we here provide a sample calculation of the adiabatic sound speed for our model.

The sound speed will need to be calculated differently depending on the relevant timescales for the astrophysical system at hand. If the sound-crossing timescale is longer than the time for weak interactions, then matter will remain in β-equilibrium as the system evolves and the proton fraction will change accordingly. This scenario may correspond to the early phases of a neutron star merger or the cooling of proto-neutron stars. Alternatively, if the dynamical timescale is shorter than the timescale required to maintain β-equilibrium, as in the late stages of a merger, the proton fraction will remain approximately constant.

In this appendix, we will calculate the sound speed for the latter case: of a system with a constant proton fraction. For such a system, the adiabatic sound speed, c_s , is defined as

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{c}_{s}}{c}\right)}^{2}\equiv \displaystyle \frac{\partial P(n,T)}{\partial \epsilon }{| }_{S,{Y}_{p}}\end{eqnarray} \tag{ 54 }$$

where $\epsilon \equiv E(n,T)n\,+\,{{mc}}^{2}n$ is the relativistic energy density, consisting of the classical internal energy density and the rest mass density. Here, we have suppressed the proton-fraction dependence of the pressure and energy models because we are considering a system that maintains its initial proton fraction, Y_p,β (n). We can expand this derivative as follows

$$\begin{eqnarray}&&\begin{array}{l}{\left(\displaystyle \frac{{c}_{s}}{c}\right)}^{2}\,=\,{\left.\displaystyle \frac{\partial P(n,T)}{\partial n}\right|}_{S,{Y}_{p}}{\left(\displaystyle \frac{\partial \epsilon }{\partial n}{| }_{S,{Y}_{p}}\right)}^{-1}\\ \,=\,{\left.\displaystyle \frac{\partial P(n,T)}{\partial n}\right|}_{S,{Y}_{p}}\times \left[\displaystyle \frac{n}{E(n,T)n+P(n,T)+{{mc}}^{2}n}\right].\end{array}\end{eqnarray} \tag{ 55 }$$

For each term in the pressure expressions of Box II, we calculate and provide the derivatives at constant entropy below.

For a tabular, cold EOS in β-equilibrium, the cold pressure derivative must be calculated numerically. For a polytropic cold EOS, however, the derivative is simply

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial {P}_{\mathrm{cold}}(n,T=0)}{\partial n}\right|}_{S,{Y}_{p}}=\displaystyle \frac{{\rm{\Gamma }}}{n}{P}_{\mathrm{cold}}(n,T=0),\end{eqnarray} \tag{ 56 }$$

where Γ is the polytropic index. In the following expression for the complete derivative, we assume the cold EOS can be represented as a polytrope. However, if this is not the case, then the first term should simply be replaced by the numerical derivative of the cold EOS. The total derivative for the case of a constant proton fraction is then

$$\begin{eqnarray}\begin{array}{rcl}{\left.\displaystyle \frac{\partial P(n,T)}{\partial n}\right|}_{S,{Y}_{p}} & = & \displaystyle \frac{{\rm{\Gamma }}}{n}{P}_{\mathrm{cold}}(n,T=0)\\ & & +\,\left\{\begin{array}{ll}16{f}_{s}\sigma {T}^{4}/(9{cn}), & n\lt {n}_{1}\\ 5T/3 & {n}_{1}\lt n\lt {n}_{2}\\ \displaystyle \frac{\partial {P}_{\mathrm{th},\deg }(n,T)}{\partial n}{| }_{S,{Y}_{p}} & n\gt {n}_{2}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 57 }$$

The degenerate thermal pressure of nucleons and electrons, ${P}_{\mathrm{th},\deg }$, is given for symmetric matter by

$$\begin{eqnarray}&&{P}_{\mathrm{th},\deg }(n,T)=-\left[{a}_{\mathrm{SM}}^{\prime} +{a}_{e}^{\prime} {Y}_{p,\beta }\right]{n}^{2}{T}^{2},\end{eqnarray} \tag{ 58 }$$

as in Box II. We assume that this is approximately equal to the β-equilibrium expression because the thermal symmetry-energy correction is small, as discussed in Section 5. This assumption will likely introduce a small error into the final sound speed, which we neglect for the present purposes.

In this appendix, for brevity, we will use the following notation: ${a}_{\mathrm{SM}}\equiv a(0.5n,0.5{M}_{\mathrm{SM}}^{* })$ and ${a}_{e}\equiv a({Y}_{p,\beta }n,{m}_{e})$. Additionally,

$$\begin{eqnarray}&&\begin{array}{l}{\left.{a}_{\mathrm{SM}}^{\prime} \equiv \displaystyle \frac{\partial a(0.5n,0.5{M}_{\mathrm{SM}}^{* })}{\partial n}\right|}_{{Y}_{q}}\\ =-\displaystyle \frac{2{a}_{\mathrm{SM}}}{3n}\left\{1-\displaystyle \frac{B}{2}\left(C+3A\right)\right\},\end{array}\end{eqnarray} \tag{ 59a }$$

$$\begin{eqnarray}&&{\left.{a}_{e}^{\prime} \equiv \displaystyle \frac{\partial a({Y}_{p,\beta }n,{m}_{e})}{\partial n}\right|}_{{Y}_{q}}=-\displaystyle \frac{2{a}_{e}}{3n}\left\{1-\displaystyle \frac{{BC}}{2}\right\},\end{eqnarray} \tag{ 59b }$$

where we have introduced

$$\begin{eqnarray}&&{\left.A\equiv A({n}_{q},{M}^{* })=\displaystyle \frac{\partial \mathrm{ln}{M}^{* }({n}_{q})}{\partial \mathrm{ln}n}\right|}_{{Y}_{q}}=-\alpha \left[1-{\left(\displaystyle \frac{{M}^{* }({n}_{q})}{{{mc}}^{2}}\right)}^{2}\right]\end{eqnarray} \tag{ 60a }$$

$$\begin{eqnarray}&&B\equiv B({n}_{q},{M}^{* })=\displaystyle \frac{{M}^{* }{\left({n}_{q}\right)}^{2}}{{M}^{* }{\left({n}_{q}\right)}^{2}+{\left(3{\pi }^{2}{n}_{q}\right)}^{2/3}{\left({\hslash }c\right)}^{2}}\end{eqnarray} \tag{ 60b }$$

$$\begin{eqnarray}&&C\equiv C({n}_{q},{M}^{* })=\displaystyle \frac{{\left(3{\pi }^{2}{n}_{q}\right)}^{2/3}{\left({\hslash }c\right)}^{2}}{{M}^{* }{\left({n}_{q}\right)}^{2}}.\end{eqnarray} \tag{ 60c }$$

For the symmetric nuclear terms, ${n}_{q}\to 0.5n$ and ${M}^{* }({n}_{q})\to 0.5{M}_{\mathrm{SM}}^{* }(0.5n)$. For the lepton term, ${n}_{q}\to {Y}_{p,\beta }n$ and the effective mass is simply the electron mass.

Using the entropy expressions of Section 6, we can then write the derivative of the degenerate thermal pressure as

$$\begin{eqnarray}&&\begin{array}{l}{\left.\displaystyle \frac{\partial {P}_{\mathrm{th},\deg }(n,T)}{\partial n}\right|}_{S,{Y}_{p}}\\ \,=\,-{\left[\displaystyle \frac{\partial ({a}_{\mathrm{SM}}^{\prime} )}{\partial n}\right|}_{S,{Y}_{p}}+\left.{\left.\displaystyle \frac{\partial ({a}_{e}^{\prime} )}{\partial n}\right|}_{S,{Y}_{p}}{Y}_{p,\beta }\right]{n}^{2}{T}^{2}\\ \quad +\,2{P}_{\mathrm{th},\deg }(n,T)\left[\displaystyle \frac{1}{n}-\displaystyle \frac{{a}_{\mathrm{SM}}^{\prime} +{a}_{e}^{\prime} \cdot {Y}_{p,\beta }}{{a}_{\mathrm{SM}}+{a}_{e}{Y}_{p,\beta }}\right].\end{array}\end{eqnarray} \tag{ 61 }$$

The second derivative of a_SM is given by

$$\begin{eqnarray}&&\begin{array}{l}{\left.\displaystyle \frac{\partial ({a}_{\mathrm{SM}}^{\prime} )}{\partial n}\right|}_{S,{Y}_{p}}={a}_{\mathrm{SM}}^{\prime} \left(\displaystyle \frac{{a}_{\mathrm{SM}}^{\prime} }{{a}_{\mathrm{SM}}}-\displaystyle \frac{1}{n}\right)+\displaystyle \frac{2{a}_{\mathrm{SM}}}{3{n}^{2}}B\\ \quad \times \,\left[3{A}^{2}-\displaystyle \frac{1}{3}B{\left(3A+C\right)}^{2}+\displaystyle \frac{1}{3}C+\displaystyle \frac{3n}{2}\displaystyle \frac{\partial A}{\partial n}\right].\end{array}\end{eqnarray} \tag{ 62 }$$

The second derivative of the electron term, for which $A\to 0$ due to its constant effective mass, is simply

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial ({a}_{e}^{\prime} )}{\partial n}\right|}_{S,{Y}_{p}}={a}_{e}^{\prime} \left(\displaystyle \frac{{a}_{e}^{\prime} }{{a}_{e}}-\displaystyle \frac{1}{n}\right)+\displaystyle \frac{2{a}_{e}}{9{n}^{2}}{BC}(1-{BC}).\end{eqnarray} \tag{ 63 }$$

Finally, the second derivative of the M* term is given by

$$\begin{eqnarray}&&\displaystyle \frac{\partial A}{\partial n}=\displaystyle \frac{\partial }{\partial n}\left[\displaystyle \frac{\partial \mathrm{ln}{M}^{* }({n}_{q})}{\partial n}\right]=\displaystyle \frac{2\alpha }{n}{\left[\displaystyle \frac{{M}^{* }({n}_{q})}{{Y}_{q}{{mc}}^{2}}\right]}^{2}\left[\displaystyle \frac{\partial \mathrm{ln}{M}^{* }({n}_{q})}{\partial n}\right],\end{eqnarray} \tag{ 64 }$$

where we have assumed that M* is defined here for symmetric matter, so that ${M}^{* }({n}_{q})\to 0.5{M}_{\mathrm{SM}}^{* }(0.5n)$.