Finite-temperature extension for cold neutron star equations of state
Carolyn A. Raithel Feryal Ozel Dimitrios Psaltis






Observations of isolated neutron stars place constraints on the equation of state (EOS) of cold, neutron-rich matter, while nuclear physics experiments probe the EOS of hot, symmetric matter. Many dynamical phenomena, such as core-collapse supernovae, the formation and cooling of proto-neutron stars, and neutron star mergers, lie between these two regimes and depend on the EOS at finite temperatures for matter with varying proton fractions. In this paper, we introduce a new framework to accurately calculate the thermal pressure of neutron-proton-electron matter at arbitrary density, temperature, and proton fraction. This framework can be expressed using a set of five physically-motivated parameters that span a narrow range of values for realistic EOS and are able to capture the leading-order effects of degenerate matter on the thermal pressure.
We base two of these parameters on a new approximation of the Dirac effective mass, with which we reproduce the thermal pressure to within ≲30% for a variety of realistic EOS at densities of interest. Three additional parameters, based on the behavior of the symmetry energy near the nuclear saturation density, allow for the extrapolation of any cold EOS in β-equilibrium to arbitrary proton fractions.  Our model thus allows a user to extend any cold nucleonic EOS, including piecewise-polytropes, to arbitrary temperature and proton fraction, for use in calculations and numerical simulations of astrophysical phenomena. We find that our formalism is able to reproduce realistic finite-temperature EOS with errors of ≲20% and offers a 1-3 orders-of-magnitude improvement over existing ideal-fluid models. 





§ INTRODUCTION
 

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 suffice, to temperatures of up to 10-100 MeV in neutron star mergers (e.g., ). Second, the composition may range from nearly pure neutron matter to symmetric matter, with some dynamical timescales 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 Fig. <ref>. For recent reviews, see e.g., .) Such extrapolations to arbitrary proton fraction and temperature add further uncertainty to the EOS and complicate numerical simulations of these phenomena. 

In the zero-temperature limit, a large number of EOS have been calculated, ranging from purely nucleonic models (e.g., ) to models incorporating quark degrees of freedom using state-of-the-art results from perturbative QCD <cit.>. 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, as were first introduced by <cit.> and <cit.>. 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 <cit.>; as well as the EOS of <cit.>, which was calculated using relativistic mean field (RMF) theory with a Thomas-Fermi approximation. More recently, the statistical model developed in <cit.> has been applied to an additional ∼10 combinations of RMF models and nuclear mass tables.



Just as parametrizations of the cold EOS have proven 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 <cit.>. The ideal-fluid thermal component is parametrized in terms of a simple adiabatic index as P_ th = ϵ_ th (Γ_ th-1), where P_ th and ϵ_ th are the thermal pressure and energy density and Γ_ th is the adiabatic index, which is assumed to be constant. Such an approach is computationally simple, but 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 <cit.>. <cit.> 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. <cit.> later expanded this work and showed that the formalism can be used to recreate even models beyond mean field theory, such as the two-loop exchange model of <cit.>. 

The strong dependence of thermal properties on the effective mass can also be seen in the behavior Γ_ th. For example, <cit.>  compared two EOS with similar zero-temperature properties but with different single-particle potentials, and hence different density-dependences 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. <cit.> 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 has important consequences for dynamical simulations as well. For example, <cit.> 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 the <cit.> and <cit.> EOS to hybrid EOS that were constructed from the zero-temperature versions of these same EOS with either Γ_ 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 what is 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 and the post-collapse accretion disk mass around the resulting black hole can differ by up to 30% when the simplified thermal effects are used <cit.>. These results all suggest that it is indeed important to account for the effect of degeneracy on the thermal pressure when simulating neutron star mergers.

The Sommerfeld expansion results of <cit.> 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 for incorporating 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. The framework we present in this paper is specific to neutron-proton-electron () matter, but 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  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 <ref> with a brief review of existing finite-temperature EOS and a discussion of the regimes in which thermal effects become important. In <ref>, we outline our model. We provide the symmetry-energy dependent extrapolation to arbitrary proton fraction in <ref>. In <ref>, we introduce our M^*-approximation of the thermal effects. We summarize the model in <ref>, in which all of the relevant equations can be found in Boxes I and II. Finally, we quantify the performance of our model in <ref>. 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. 



§ OVERVIEW OF FINITE-TEMPERATURE EOS


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

Two of the most widely-used finite-temperature EOS are the models of <cit.>, which is based on a finite-temperature liquid drop model with a Skyrme nuclear force, and <cit.> which is an RMF model that is extended with the Thomas-Fermi approximation. An additional eight EOS have been calculated with the framework of <cit.>, 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 an RMF model and additionally includes excluded volume effects. Of the RMF models that have been used with the HS method, six are nucleonic: TMA <cit.>, TM1 <cit.>, NL3 <cit.>, FSUGold <cit.>, IUFSU <cit.>, DD2 <cit.>; while the models BHBΛϕ and BHBΛ include hyperons with and without the repulsive hyperon-hyperon interaction mediated by the ϕ meson, respectively <cit.>. Additionally, <cit.> created a set of two finite-temperature EOS, SFHo/x, that also used the statistical method of HS, but with 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 models FSUGold <cit.> and NL3 <cit.>. Tables of these various EOS can be found on the website of M. Hempel,[https://astro.physik.unibas.ch/people/matthias-hempel/equations-of-state.html] , and/or the  database.[https://compose.obspm.fr/home/]

More recently, several new finite-temperature EOS have been added to the  database. These include the SLY4-RG model, which is calculated in nuclear statistical equilibrium using a Skyrme energy functional <cit.>, as well as chiral mean field theory models, which include hyperons as additional degrees of freedom <cit.>, generalized relativistic density functional models <cit.>, and models calculated using a variational method applied to two- and three-body nuclear potentials <cit.>.

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 as well as the eight nucleonic EOS 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 <cit.> and have been used in many subsequent works <cit.>. In these hybrid EOS, the thermal pressure is written as

    P_ th, hybrid(n,T) = n E_ th, hybrid(n,T)  ( Γ_ th - 1),

where E_ th,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 ≤Γ_ th≤ 2. In the hybrid approximation, Γ_ th is assumed to be constant.

Following <cit.>, the hybrid temperature-dependence of E_ th,hybrid is included as an ideal fluid plus a contribution from relativistic particles, i.e.,

    E_ th,hybrid(n,T) = 3/2 k_B T + 4 σ/cf_s/ n T^4,

where k_B is the Boltzmann constant, T is the temperature, and σ≡π^2 k_B^4 / [60 ħ^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_BT≪ 2 m_e c^2, where m_e is the mass of an electron, photons will dominate and f_S=1. For k_BT≫ 2 m_e c^2, electrons and positrons become relativistic as well and yield f_S = 1 + 2× (7/8) = 11/4. Finally, for k_BT ≳ 10 MeV, thermal neutrinos and anti-neutrinos appear, rendering f_s = 11/4 +  3 × (7/8) = 43/8. If right-handed neutrinos were to exist, this would become f_s = 11/4+ 3×2×(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 in wherever we use STOS in this paper. For the STOS thermal lepton and photon contribution, we use eq. (<ref>) with the appropriate lepton density. For the cold lepton energy, we add the contribution for a degenerate gas of relativistic electrons. Because all the EOS neglect neutrinos, we will also neglect neutrinos in our comparisons and thus we will calculate f_S only as 

    f_S = 
    
    1,       k_B T < 1 MeV , 
    
     11/4,      k_B T ≥ 1 MeV.


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, 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 Fig. <ref>.  In this plot, we show various regions calculated for the EOS STOS, all 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.[We note that, throughout this paper, we use the coldest HS calculation, performed at k_BT=0.1 MeV, as an approximation of the zero-temperature EOS. Even though the STOS EOS is calculated at T=0 MeV, we use the k_BT=0.1 MeV table as our cold component for this EOS as well, in order to maintain consistency with the HS set of EOS.] Finally, the thermal contribution, P_ th, is defined as P_ total -P_ cold for the same proton fraction.

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_ th, ideal=nK_BT), 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, we also show in Fig. <ref> 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  neutron stars, as simulated using the EOS STOS <cit.>. 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 <cit.>. 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 such 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 towards Y_p=0.5, the white degeneracy region in Fig. <ref> 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. 




§ GENERIC MODEL OF A FINITE TEMPERATURE EOS


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 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 <ref>.

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

    E_ nucl(n, Y_p, T) = E_ nucl(n, Y_p=12, T) + E_ sym(n, T)(1-2 Y_p)^2,

where E_ nucl(n, Y_p = 12, T) represents the energy of symmetric nuclear matter and

    E_ sym(n,T) ≡1/2∂^2 E_ nucl(n, Y_p, T)/∂ (1-2 Y_p)^2|_Y_p=1/2

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

    Y_p = n_p/n = N_p/N_n + N_p,

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

    n_e = Y_p n.

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

    n_n = (1-Y_p) n.


We can further expand eq. (<ref>) by separating the energy of cold, symmetric matter from its thermal contribution, i.e., 

    E_ nucl(n, Y_p, T) =     E_ nucl(n, Y_p=12, T=0)  
    
    	 +     E_ nucl,th(n, Y_p=12,T) 
    
    	 +     E_ sym(n, T)(1-2 Y_p)^2.

Here and throughout the 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 eq. (<ref>) to yield

    E_ nucl(n, Y_p, T) =     E_ nucl(n, Y_p,β, T=0)  
    
    	 +     E_ nucl,th(n, Y_p=12,T) 
    
    	 +     E_ sym(n, T)(1-2 Y_p)^2 
    
    	 -     E_ sym(n, T=0)(1-2 Y_p,β)^2,

where Y_p,β 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,β =Y_p,β(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

    E_ lepton(n, Y_p, T=0) =  3 K Y_p (Y_p n)^1/3,

where the extra factor of Y_p comes from our definition of E as the energy per baryon, combined with eqs. (<ref>) and (<ref>). Here, K≡ (3 π^2)^1/3 ( ħ c/4). Additionally, there will also be a thermal contribution, E_ lepton,th(n,Y_p,T), which we derive in <ref>.

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



    E(n, Y_p, T)   = E(n, Y_p, T=0) + E_ th(n, Y_p, T)
    
    
        E(n, Y_p, T=0)  = E(n, Y_p,β, T=0)  
           + E_ sym(n, T=0)[(1-2 Y_p)^2 - (1-2 Y_p,β)^2] 
          + 3 K  ( Y_p^4/3 - Y_p,β^4/3) n^1/3
        E_ th(n, Y_p, T)  = E_ nucl,th(n, Y_p=12,T) 
            + E_ lepton,  th(n, Y_p, T) 
            + E_ sym,th(n, T)(1-2 Y_p)^2 .




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

    P ≡ -∂ U/∂ V|_N_q,S = n^2 [ ∂ E(n,T=0)/∂ n] |_Y_p, S

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 eq. (<ref>), 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 <ref>.

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, 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. 



§ DERIVATION OF THE COLD SYMMETRY ENERGY IN THE FERMI GAS LIMIT


We turn first to the symmetry energy correction term, E_ sym(n, T) of eq. (<ref>). 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 eq. (<ref>), we perform the expansion with respect to the proton fraction and, in the following section, 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_ sym(n,T) at zero-temperature. For the thermal contribution to the symmetry energy, which turns out to be negligible, see <ref>.

It is particularly useful to parameterize the symmetry energy in terms of its separate kinetic and potential components at zero-temperature <cit.>, 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 <cit.>. In this framework, we parameterize the symmetry energy of eq. (<ref>) as

    E_ sym(n, T=0)= η E_ sym^ kin(n)  
    	+ [ S_0 - η E_ sym^ kin(n_ sat) ] (n/n_ sat)^γ,

as in <cit.>. Here, E_ sym^ 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^-3 is the nuclear saturation density,[We note that n_ sat does vary slightly among the EOS in our sample, but we fix the value to n_ sat=0.16 fm^-3 in order to more easily compare the various EOS. We find that this does not significantly affect the results.]  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 ≡ E_ sym(n_ 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

    ε_k,q/n = 3/5 E_f(n_q)

where ε_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,

    E_f(n_q) = ħ^2/2 m( 3π^2 n_q )^2/3,

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

    E_ sym^ kin(n)     = 3/5[ 2E_f(n_p = n_n =1/2 n) - E_f(n_n = n)  ]  
       = 3/5(2^1/3 -1 )  E_f(n).


We can also eliminate the parameter η in eq. (<ref>) by introducing the constant L, which is related to the overall slope at the saturation density via,

    L ≡  3  n_ sat[ ∂ E_ sym(n, T=0)/∂ n]|_n_ sat.

Combining eqs. (<ref>) and (<ref>), we can solve for η in terms of the quantities S_0 and L, which are constrained by nuclear physics experiments for matter near Y_p=1/2 <cit.>. We find

    η = 5/9[  L-3 S_0 γ/(2^1/3-1)(2/3 - γ) E_f(n_ sat)]  ,

thereby leaving one free parameter, γ, which is constrained by nuclear experiments to lie in the range ∼ 0.2 to 1.2 (see, e.g., Fig. 2 of ; ).

llll

0.48 
 Symmetry energy parameters characterizing each EOS at k_BT=0.1 MeV.  


 EOS      
S_0 (MeV)      
L (MeV)      
γ   

TM1      36.95     110.99      0.75 

          TMA      30.66     90.14      0.66 

          NL3      37.39     118.49      0.62 

          FSG      32.56     60.43      1.11 

          IUF      31.29     47.20      0.52 

          DD2      31.67     55.03      0.91 

        STOS      36.95     110.99      0.77 

         SFHo      31.57     47.10      0.41 

         SFHx      28.67     23.18      -0.04[The inferred value for γ for SFHx is highly sensitive to the density range that is included in the fit; see the discussion in the text for details.] 

	 LS      29.3      74.0	     1.05 

 	SLY4-RG    	 32.04	   	46.00    	0.35 

S_0 and L are fixed to the values predicted for each EOS, while γ is a fit parameter. All fits are performed for densities above n≥0.01 fm^-3 and n_ sat=0.16 fm^-3.


We thus have a complete expression for the symmetry energy that depends only on the three parameters γ, S_0, 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  matter in neutrinoless β-equilibrium,



    Y_p,β/(1-2 Y_p,β)^3 = 64/3 π^2 n[ E_ sym(n, T=0)/ħ c]^3

<cit.>.  When solved for Y_p,β, this becomes

    Y_p,β = 1/2 + (2 π^2)^1/3/32n/ξ{ (2π^2)^1/3 - ξ^2/n[ħ c/E_ sym(n,T=0)]^3 },

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

    ξ≡[ E_ sym(n,T=0)/ħ c]^2  ×
    { 24 n [ 1+ √( 1 +  π^2 n/288(ħ c/E_ sym(n,T=0) )^3)]  }^1/3.


For each of the EOS in our sample, we stitch together a complete cold EOS at β-equilibrium from the publically-available tables at fixed Y_p, by requiring that μ_e + μ_p - μ_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 eqs. (<ref>)-(<ref>) and keeping S_0 and L fixed for each EOS. We perform the fits using a standard least-squares method and limit the density range to n ≥ 10^-2 fm^-3. In principle, eqs. (<ref>)-(<ref>) apply only to  matter, which will be uniform only above 0.5 n_ sat. However, in practice, we find a very small difference in the fits for γ whether we include densities above 0.5n_ sat=0.08 fm^-3 or whether we start the fits at a slightly lower but still astrophysically relevant cutoff of n=10^-2 fm^-3. We show the resulting fit values in Table <ref>. 

We note that the range of EOS provided in Table <ref> is intentionally broad. While the symmetry energy parameters of some of these EOS disagree with the combined set of experimental constraints (see  for a recent review), or are in disagreement with certain theoretical considerations such as chiral effective field theory results for pure neutron matter <cit.>, they are all consistent with at least some experimental constraints on S_0 and L.

We find that γ spans roughly the range of experimentally-allowed values, between 0.15 and 1.0, as expected, with the exception of SFHx. SFHx has an extremely low value of L, which 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^-2 fm^-3 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 eq. (<ref>) to fit for γ from the β-equilibrium proton fractions of realistic EOS. We wish to also emphasize that eq. (<ref>) can, of course, be used to calculate Y_p,β, given a choice of S_0, L, and γ. Once these three parameters are specified, eqs. (<ref>)-(<ref>) can be used to calculate Y_p,β for any EOS. As a result, 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 Fig. <ref> for the EOS NL3  <cit.> and DD2 <cit.>. 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 <ref>). The top panel of Fig. <ref> 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 eqs. (<ref>)-(<ref>), and correcting for the leptons, all according to eq. (<ref>). For these models, we take the values of S_0, L, and γ for each EOS from Table <ref>. We note that we are plotting pressures, but could have similarly shown the energy. We use eq. (<ref>) to convert the equations of this section to pressures; for the complete set of pressure expressions, see <ref> and Box II. 


The bottom panel of Fig. <ref> 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 the residuals between our model and NL3 and DD2 are comparable to those shown in Fig. <ref>. We, therefore, 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_0, L, and the narrowly-constrained parameter γ. There are two possible routes for creating 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_0, L, and Y_p,β. In this case, eq. (<ref>) can be used to fit for γ. We have provided such fits for the EOS in our sample in Table <ref>. 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_0, 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.



§ THERMAL CONTRIBUTION TO THE ENERGY


We now turn to the thermal energy, which was first defined in eq. (<ref>) as

    E_ th(n, Y_p, T)     = E_ nucl, th(n, Y_p=12,T) 
        + E_ sym, th(n,T) × (1-2 Y_p)^2 
        + E_ lepton, th(n, Y_p, T) .


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 Fig. <ref> for the markedly different behaviors in thermal pressure across these three regimes.) 

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

    E_ th   (n, Y_p, T)  =  
       
    E_ rel(n,T),       n_1< n < n_1 
    
    E_ ideal(T),     n_1 < n <n_2 
    
    E_ th, deg.(n,Y_p=1/2,T) 
      + E_ sym, th(n, T)  (1-2 Y_p)^2,    n_1< n>n_2

where the relativistic component,

    E_ rel(n,T) = 4 σ/cf_s/ n T^4,

and the ideal component,

    E_ ideal(T) = 3/2 k_BT

are given as in eqs. (<ref>) and (<ref>). Here, E_ 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_1, as the density at which the relativistic and ideal-fluid energies are equal. The second transition density, n_2, 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 the 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

    E_ th   (n, Y_p, T) ≈ E_ rel(n,T) 
       + [ E_ ideal(T)^-1 +  E_ th, deg.(n,Y_p,T)^-1]^-1,
 
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, as it does not require the calculation of transition densities, which will vary with the temperature and proton fraction.

In order 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

    E_ th, q^ deg(n, Y_q, T)    =  a(Y_qn, M^*) (N_q/N_p + N_n)T^2,
        = a(Y_qn, M^*)  Y_q T^2

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

    a(n_q, M^*) ≡π^2 k_B^2/2√(( 3 π^2 n_q )^2/3  (ħ c)^2 + M^*(n_q)^2 )/( 3 π^2 n_q )^2/3 (ħ c)^2 ,

where M^*(n_q) is the Dirac effective mass of the relevant species at a specific density. <cit.>.

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

    E_ th, nucl^ deg(n, T)    =  [ a(n_p, M^*_ p, SM)N_p + a( n_n, M^*_ n, SM)N_n /N_p + N_n]  T^2 
       = a(0.5 n, 0.5 M^*_ SM) T^2,

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, SM≈ M^*_ p, SM≈12 M^*_ SM.

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

    E_ sym, th(n, T) =  
       
     0,        n < n_2 
    [ a(n, M^*_ PNM) - a(0.5 n, 0.5 M^*_ SM) ] T^2,      n >  n_2,

 where the low-density limit of E_ sym, th arises from the fact that both pure neutron matter and symmetric matter behave identically as ideal or relativistic fluids at n < n_2. 

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_ th, nucl(n, Y_p, T)≈ E_ th, nucl(n, Y_p=12, 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 remainder of the 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 ≈ m_e. This allows us to write eq. (<ref>) simply as

    E_th, e^-^ deg(n, Y_p, T) = a(Y_p n, m_e) Y_p T^2,
 
where we have required that the electron fraction balance the proton fraction in order to satisfy the requirement of charge neutrality and we have used eq. (<ref>) to substitute Y_p. We note that in the presence of a significant population of positrons, the proton fraction in eq. (<ref>) 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 eq. (<ref>) for E_ th  as follows:

    E_ th(n, Y_p, T)  =  
       
    4 σ f_s  T^4/ (c n),       n_1< n < n_1 
    
    (3/2) k_B T,     n_1< n <n_2 
    [ a(0.5 n, 0.5 M^*_ SM) + a(Y_p n, m_e) Y_p   ] T^2 ,    n_1< n>n_2

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 only of the density, temperature, proton fraction, and the effective mass of the nucleons in symmetric matter. 




 §.§ M^*-approximation


A full calculation of E_ th using eq. (<ref>) 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,

    M^*(n_q) = { (m c^2)^-b + [ mc^2 (  n_q /n_0 )^-α]^-b}^-1/b,

where m is the nucleon mass (which we take to be the neutron mass, mc^2=939.57 MeV)[The EOS in our sample vary in their low-density limit of M^* from 938-939.57 MeV. This parameter can easily be adjusted to any low-density value for M^*. For simplicity, however, we take it to simply be the neutron mass. We find that this simplification has a negligible effect on our results. ] and n_0 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 thus fix it to this value in the following analysis, leaving just two free parameters to describe the effective mass, M^* = M^*(n_0,α). 

We fit the effective masses together at k_BT=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 models LS and SLY4-RG here because the effective masses for these EOS are not currently published (but see <ref> for a separate comparison with these models). The results of these fits are given in Table <ref> for symmetric matter. For completeness, we also include in Table <ref> the fits for pure neutron matter, which can be used to calculate E_ sym,th(n,T) in eq. (<ref>).

lcccccc

0.48  
 Parameters characterizing M^*, fit together at k_BT=1, 10, and 47.9 Mev, for either pure neutron matter (PNM) or symmetric matter (SM).  


    
PNM (Y_p = 0.01)   
    
    
SM (Y_p= 0.5)   
     
 


 EOS      
n_0 (fm^-3)      
α      
n_0 (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
We fix b=2 and m=m_n in all fits.


We show the performance of the fit for NL3 in Fig. <ref>. In this fit, we use the NL3 tables calculated at k_BT=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 (top panel), and Y_p=0.5, to represent symmetric matter (bottom panel). We show our approximation for M^* as the black solid line. We find that the M^*-approximation accurately captures the behavior predicted by the full EOS, with fit parameters n_0=0.10 fm^-3 and α=0.90 for Y_p=0.01 and fit parameters n_0=0.11 fm^-3 and α=1.08 for Y_p=0.5. 








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 the bottom panel of Fig. <ref>, at n_ sat/2 ≈ 0.08fm^-3). 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. <cit.>, <cit.>, and <cit.> 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 <cit.>. Because all 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 re-calculating all EOS with a different formalism and is beyond the scope of this paper. However, we note that at high temperatures (k_BT≳15 MeV), the non-uniform phase of matter disappears (see discussion around Fig. 5 in ). 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_BT=47.9 MeV curve for M^* or we fit the curves for all the temperatures together, the difference in the resulting parameters is small. We, therefore, choose to perform the fits to three temperatures (k_BT=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 Fig. <ref> 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 ∈ (0.08,0.22) fm^-3, with an average value of ∼0.13 fm^-3 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_0 and α on the temperature, thus suggesting that these parameters could be treated as constants for use in numerical simulations. 



 §.§ 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 <ref>. As in <ref>, we make the comparison in terms of the pressure, rather than the energy, and use eq. (<ref>) to convert between the two. The expressions for P_ th(n, Y_p, T) are given in Box II in <ref>. In particular, all results shown here use the smoothed approximation of the thermal pressure, as defined in eq. (<ref>).

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 Fig. <ref> for NL3. We find that our approximation of P_ th closely recreates the full calculation for NL3 for nearly all densities and temperatures explored here. For comparison, we also include in Fig. <ref> the hybrid approximation with Γ_ th=1.67 as dashed lines.[We choose the relatively low value of Γ_ th=1.67 in order to minimize the residuals of the hybrid model. This value of Γ_ th ensures the hybrid EOS matches an ideal fluid at intermediate densities. Larger values, as are more commonly used in numerical simulations, would cause the hybrid P_ th to overestimate even the ideal regime.]  The full thermal pressure agrees with the hybrid approximation only 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.


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 Fig. <ref>, we zoom in on P_ th at k_BT=10 MeV and Y_p=0.5 and show the effect of varying the parameters n_0 and α for symmetric matter. We intentionally take extreme values for the parameters, well beyond the ranges found in Table <ref>, in order to emphasize that the variations between more realistic parameter choices will be small. Even for these unreasonable choices of values for n_0 and α, we find that P_ th approximates the full thermal pressure reasonably well and, in all cases, better than the ideal fluid approximation. Analyzing the specific dependences more closely, we see in Fig. <ref> that the parameter n_0 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 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, 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 those interactions.

Finally, we compare the M^*-approximation of the thermal pressure against the full sample of EOS listed in Table <ref>. We show the corresponding residuals at three temperatures in Fig. <ref> and find that the residuals are typically ≲30% at densities above 0.5 n_ sat. For comparison, Fig. <ref> also shows a sample set of residuals between the full thermal pressure from NL3 and the hybrid approximation (Γ_ th=1.67) as the 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.



 §.§ 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  <ref>). We also include here the two-loop exchange model of <cit.>, which is an extension of mean field theory. We note that the pressures of the <cit.> 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 comparisons in this section are made at Y_p=0.5 and T=20 MeV. We also fix n_0 and α to the mean values for symmetric matter from Table <ref> for all three EOS.




Figure <ref> 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.

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^-3, the hybrid approximation produces slightly smaller residuals in the thermal pressure for these non-RMF models. In this regime, the hybrid approximation tends to over-estimate 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^-3.




§ PUTTING IT ALL TOGETHER


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




 0.9
 Box I: Total Energy Expressions for Finite-Temperature Dense Gas. 


The energy per particle of  matter is given by

    E(n, Y_p, T)      = ( Cold EOS in β-equilibrium) 
    		+ 3 K  ( Y_p^4/3 - Y_p,β^4/3) n^1/3
        + E_ sym(n, T=0)[ (1-2 Y_p)^2 - (1-2 Y_p,β)^2] 
     	 +	
    
    4 σ f_s  T^4/ (c n),       n_1< n < n_1 
    
    (3/2) k_B T,     n_1< n <n_2 
    [ a(0.5 n, 0.5 M^*_ SM) + a(Y_p n, m_e) Y_p   ] T^2 ,    n_1< n>n_2,



where the symmetry energy is approximated as



    E_ sym(n, T=0) =  η E_ sym^ kin(n) + [ S_0 - η E_ sym^ kin(n_ sat)](n/n_ sat)^γ,


    E_ sym^ kin(n) = 3 /5( 1 -2^1/3) E_f(n),


    η = 5/9[  L-3 S_0 γ/(1-2^1/3)(2/3 - γ) E_f(n_ sat)],



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



    a(n_q, M_q^*) ≡π^2 k_B^2/2√( M_q^*(n_q)^2 + (3 π^2 n_q)^2/3(ħ c)^2 )/(3 π^2 n_q)^2/3(ħ c)^2   and  
    M^*(n_q) = { (m c^2)^-b + [ mc^2 (  n_q /n_0 )^-α]^-b}^-1/b.



  * The parameters S_0, L, and γ∈ (0.2-1.2) are freely specified; this will uniquely specify Y_p,β. 

  * Alternatively, S_0, L, and Y_p,β may be specified and the proton fraction may be fit for γ. We provide fits to γ for eleven EOS in Table <ref>.

  * We find that for M^*_ SM, n_0 ∼ 0.13 fm^-3 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 eq. (<ref>), 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

    S(n,N_p, N_n,N_e, T) =
     
    S_ rel,       n_1< n < n_1 
    
    S_ ideal ,     n_1< n <n_2 
    
    S_ deg ,    n_1< n>n_2,

where n_1 and n_2 are the thermal energy transition densities, as defined in <ref>.


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

    S_ rel =  16 σ f_s/3 c( N_p + N_n/n) T^3.


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

    S_ ideal =  ( N_p + N_n + N_e ) k_B 
    ×{ln[ ( N_p+N_n/N_p + N_n + N_e) n^-1( m k_B T/ 2πħ^2)^3/2] + 5/2}.


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

    S_q = 2 a_q N_q T

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

    S_ deg =2 { a(0.5 n, 0.5 M^*_ SM)[N_n + N_p] + a(Y_p n, m_e) N_e } T.
 


We summarize the resulting pressure equations in Box II.




 0.9
 Box II: Pressure Expressions for Finite-Temperature Dense Gas. 


The pressure of  matter is given by

    P(n, Y_p, T)     = ( Cold EOS in β-equilibrium) 
    	+ K  ( Y_p^4/3 - Y_p,β^4/3) n^4/3
        + P_ sym(n, T=0)[ (1-2 Y_p)^2 - (1-2 Y_p,β)^2 ] 
     +
     
    4 σ f_s  T^4/ (3 c),       n_1< n < n_1 
    
    n k_B T,     n_1< n <n_2 
    
    -[ ∂ a(0.5 n, 0.5 M^*_ SM)/∂ n + ∂  a(Y_p n, m_e)/∂ n Y_p    ] n^2 T^2 ,    n_1< n>n_2,



where n_1 and n_2 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 



    P_ sym(n, T=0)     = 2η/3  n E_ sym^ kin(n) + 
    	[ S_0 - η E_ sym^ kin(n_ sat)] ( n/n_ sat)^γγ n  .



The full analytic expression for Y_p,β is given in eq. (<ref>) and derived in Appendix A.


The M^*-approximation derivatives are given by



    ∂ a(n_q, M^*)/∂ n|_Y_q = 
    -2 a(n_q, M^*)/3 n{ 1 -1/2[ M^*(n_q)^2/ M^*(n_q)^2+(3π^2 n_q)^2/3(ħ c)^2] ( (3π^2 n_q)^2/3(ħ c)^2/M^*(n_q)^2+ 3 ∂ln [ M^*(n_q)]/∂ln n|_Y_q) },

  and

    ∂ln[ M^*(n_q)]/∂lnn|_Y_q = -  α[ 1- (M^*(n_q)/Y_q mc^2)^2 ],

where, for symmetric matter, we replace M^*(n_q)→ 0.5 M^*_ SM(0.5 n) and for the electrons, M^*(n_q) → m_e.


  * As in Box I., there are five free parameters:  S_0, L, γ, n_0, α. 

  * A user may freely specify S_0, L, and Y_p,β and fit for γ. Alternatively, a user may specify  S_0, L, and γ, which will uniquely specify Y_p,β.

  * We provide fits for γ, n_0, and α for the EOS in our sample in Tables <ref> and <ref>.









The piecewise definitions of the thermal energy and pressure are mathematically convenient, but the sharp transitions are themselves unphysical, as discussed in <ref>. We, therefore, instead implement the thermal pressure using a smoothed approximation of the form

    P_ th(n, Y_p, T) ≈ P_ rel + (  P_ ideal^-1 +  P_ deg^-1)^-1.

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 eq. <ref>) in order to keep the problem tractable. However, this is not mathematically exact since, 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, but 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 sub-luminal 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.



§ COMPLETE MODEL:  COMPARISON OF REALISTIC EOS AT ARBITRARY Y_P AND T


In <ref>, 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 <ref>, 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 <ref> 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_0, α, and γ listed in Tables <ref> and <ref> for each EOS. We show the results as the solid lines in Fig. <ref>, 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. 

Figure <ref> 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_0, α and γ listed in Tables <ref> and <ref>, where possible. For LS and SLY4-RG, for which we do not have fit values for n_0 and α, we use the average parameter values for symmetric matter in Table <ref>. 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.

For all the EOS in our sample, the error introduced by our model increases in the vicinity of ∼ 0.5 n_ sat. This is a result of the break-down in the E_ sym approximation at low densities. Our derivation of E_ sym in <ref> assumed uniform  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  matter that probe different physics through the choice of n_0, α, γ, S_0, and L. Our model allows further freedom in creating 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 that can be integrated with minimal computational cost to a large array of numerical calculations.




§ CONCLUSIONS


In this paper, we have developed a general framework for calculating the pressure of neutron-star matter at arbitrary proton fraction and finite temperature. Our model is designed so that the corrections we have developed here can be added to any cold  EOS in neutrinoless β-equilibrium. The model is based on a set of five physically-motivated parameters: S_0, L, γ, n_0, and α. The first three, S_0, 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_0 and α are introduced through our M^*-approximation, where n_0 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∼ 0.13 fm^-3 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 used to 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.

Acknowledgements. We thank Vasileios Paschalidis for useful discussions and comments on this work. CR is supported by the NSF Graduate Research Fellowship Program Grant DGE-1143953. FO and DP acknowledge support from NASA grant NNX16AC56G.



apj






§ A. RELATIONSHIP BETWEEN THE Β-EQUILIBRIUM PROTON FRACTION AND THE SYMMETRY ENERGY FOR  MATTER 


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

At zero-temperature for  matter, the total energy per baryon is given by

    E_ tot(n, Y_p) = E_n(n, Y_n) + E_p(n, Y_p) + E_e(n, Y_e),

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.

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

    ∂ E_ tot(n, Y_p) /∂ Y_p  = 
     ∂ E_n/∂ Y_n∂ Y_n/∂ Y_p +
     ∂ E_p/∂ Y_p  +
     ∂ E_e/∂ Y_e∂ Y_e/∂ Y_p,

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 μ_i ≡∂ E_i/ ∂ Y_i |_S,n, eq. (<ref>) simplifies to

    ∂ E_ tot(n, Y_p) /∂ Y_p = -μ_n + μ_p + μ_e,

which is zero in β-equilibrium. 

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

    E_ tot(n, Y_p) = E_ nucl(n, 1/2) + E_ sym(n) (1-2 Y_p)^2 + E_e(n, Y_e).

This results in

    ∂ E_ tot(n, Y_p) /∂ Y_p = -4 (1-2 Y_p)  E_ sym(n) +  ∂ E_e/∂ Y_e∂ Y_e/∂ Y_p.

By charge neutrality and the definition of the chemical potential, the second term is simply μ_e. 
Combining eqs. (<ref>) and (<ref>) in β-equilibrium gives

    μ_e = 4(1-2 Y_p,β) E_ sym(n).

For relativistic electrons, 

    μ_e = √(p_f^2c^2 + m_e^2 c^4)≈ p_fc,

where p_fc = (3 π^2 Y_e n)^1/3ħ c is the Fermi momentum of the electrons. Combining this expression for μ with eq. (<ref>) yields

    (3 π^2 Y_p,β n)^1/3ħ c = 4 (1-2 Y_p,β)E_ sym(n),

or, rearranged to match eq. (<ref>),

    Y_p,β/(1-2 Y_p,β)^3 = 64 E_ sym(n)^3/ 3 π^2 n (ħ c)^3.

Solved for Y_p,β, this gives

    Y_p,β = 1/2 + (2 π^2)^1/3/32n/ξ{ (2π^2)^1/3 - ξ^2/n[ħ c/E_ sym(n,T=0)]^3 },

as in eq. (<ref>), and where ξ is defined as

    ξ≡[ E_ sym(n,T=0)/ħ c]^2  
    { 24 n [ 1+ √( 1 +  π^2 n/288(ħ c/E_ sym(n,T=0) )^3)]  }^1/3,

as in eq. (<ref>).

Thus, if the form of E_ sym(n) is known, this will uniquely specify Y_p,β. Alternatively, if the β-equilibrium proton fraction is known from the cold EOS, 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_0, and L can be used to fit for the parameter γ; or, specifying S_0, L, and γ can be used to calculate Y_p,β. The latter option is particularly useful as 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,β, which is required to calculate the sound speed in Appendix B. The derivative at constant entropy is given by

    ∂ Y_p,β/∂ n|_S = 1/16(π/2)^2/3{ (2π^2)^1/3( 1 - n ϕ/ξ) + ξ[ ħ c/E_ sym(n,T=0)]^3 ( 3x -ϕ) }

where for simplicity we have introduced the quantities

    ϕ = 1/2 n( 3 x n + 1 ) +
     1/6n{ 1+π^2 n/288[ ħ c/E_ sym(n,T=0)]^3 }^-1/2 (3xn - 1)

and

    x = 1/n^2P_ sym(n,T=0)/E_ sym(n,T=0).





§ B. CALCULATION OF THE SOUND SPEED





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

    ( c_s/c)^2 ≡∂ P(n, T)/∂ϵ|_S,Y_p

where ϵ≡ E(n,T)n + m c^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

    ( c_s/c)^2  = ∂ P(n, T)/∂ n|_S,Y_p( ∂ϵ/∂ n|_S,Y_p)^-1 =  ∂ P(n,T)/∂ n|_S,Y_p[ n/E(n, T) n + P(n, T) + mc^2n],

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

    ∂ P_ cold(n, T=0)/∂ n|_S,Y_p = Γ/n P_ cold(n, T=0),

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, 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


    ∂ P(n, T)/∂ n|_S,Y_p =
        Γ/n P_ cold(n, T=0) 
     + 
    
     16 f_s σ T^4/(9c n) , 						    n_1< n < n_1 
    
     5T/3 								    n_1 < n < n_2 
    ∂ P_ th, deg(n,T)/∂ n|_S,Y_p   n_1<  n > n_2 .


The degenerate thermal pressure of nucleons and electrons, P_ th, deg, is given for symmetric matter by

    P_ th, deg(n, T)  = - [ a_ SM' + a_e'  Y_p,β] n^2 T^2,

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 <ref>. 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_ SM≡ a(0.5 n, 0.5 M^*_ SM) and a_ e≡ a(Y_p,β n, m_e). Additionally,


    a_ SM'   ≡∂ a(0.5 n, 0.5 M^*_ SM)/∂ n|_Y_q = 
    -2 a_ SM/3 n{ 1 -B/2(C+ 3 A ) } ,
    
    
    a_ e'   ≡∂ a(Y_p,β n, m_e)/∂ n|_Y_q = 
    -2 a_ e/3 n{ 1 -BC/2} ,


where we have introduced


    A    ≡  A(n_q, M^*) = ∂ln M^*(n_q)  /∂lnn|_Y_q = -  α[ 1- ( M^*(n_q)/mc^2)^2 ] 
    
    
    B    ≡  B(n_q, M^*) =  M^*(n_q)^2/M^*(n_q)^2 + (3π^2 n_q)^2/3(ħ c)^2
    
    C    ≡  C(n_q, M^*) =(3π^2 n_q)^2/3(ħ c)^2/M^*(n_q)^2.


For the symmetric nuclear terms, n_q → 0.5 n and M^*(n_q) → 0.5 M^*_ SM(0.5 n). For the lepton term, n_q → Y_p,βn and the effective mass is simply the electron mass.

Using the entropy expressions of <ref>, we can then write the derivative of the degenerate thermal pressure as

    ∂ P_ th, deg(n, T)/∂ n|_S,Y_p =    -[ ∂ (a_ SM')/∂ n|_S,Y_p + ∂ (a_ e') /∂ n|_S,Y_p Y_p,β] n^2T^2 
    + 2 P_ th, deg(n,T) [ 1/n -  a_ SM'  + a_ e'· Y_p,β/a_ SM +a_ e Y_p,β]


The second derivative of a_ SM is given by












    ∂ (a_ SM') /∂ n|_S,Y_p = a_ SM'( a_ SM'/a_ SM - 1/n) + 2a_ SM/3n^2 B [ 3 A^2  - 1/3 B(3A+C)^2 + 1/3 C + 3n/2∂ A/∂ n]

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


    ∂ (a_ e')/∂ n|_S,Y_p = a_ e'( a_ e'/a_ e -1/n)
    + 2a_ e/9n^2 
    BC(1-BC)

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

    ∂ A/∂ n 
    	= ∂/∂ n[ ∂ln M^*(n_q)/∂ n] 
    	= 2 α/n[ M^*(n_q)/Y_qmc^2]^2[ ∂ln M^*(n_q)/∂ n],

where we have assumed that M^* is defined here for symmetric matter, so that M^*(n_q) → 0.5 M^*_ SM(0.5 n).


