A New Equation of State for Astrophysical Simulations
Abstract
We generate a new complete equation of state (EOS) of nuclear matter for a wide range of temperatures, densities, and proton fractions ready for use in astrophysical simulations of supernovae and neutron star mergers. Our previous two papers tabulated the EOS at over 180,000 grid points in the temperature range = 0 to 80 MeV, the density range = 10 to 1.6 fm, and the proton fraction range = 0 to 0.56. In this paper we combine these data points using a suitable interpolation scheme to generate a single equation of state table on a finer grid. This table is thermodynamically consistent and conserves entropy during adiabatic compression tests. We present various thermodynamic quantities and the composition of matter in the new EOS, along with several comparisons with existing EOS tables. Our equation of state table is available for download.
pacs:
21.65.Mn,26.50.+x,26.60.Kp,21.60.JzI Introduction
The equation of state (EOS), pressure as a function of density and temperature, for warm dense nuclear matter is an important ingredient in the theory of core collapse supernovae MPAreview (). The EOS influences the shock formation and evolution, and determines the compactness of the nascent protoneutron star. The EOS also impacts the position of the neutrinosphere, which is the site of last scattering for neutrinos in supernovae, and the emitted neutrino spectra light07 (). This helps determine neutrino flavor oscillations Duan10 (), and ultimately the spectra of supernova neutrinos in terrestrial detectors.
In a previous paper SHT10a () we used a relativistic mean field (RMF) model to selfconsistently calculate nonuniform matter at intermediate density and uniform matter at high density. In a following paper SHT10b (), we used a Virial expansion for a nonideal gas of nucleons and nuclei to obtain the EOS at low densities. Our Virial approach Virial1 (); Virial2 (); Virial3 () uses elastic scattering phase shifts and nuclear masses as input. We also include Coulomb corrections that can be important for neutrino interactions ioncorrelations (). Alternative approaches to the equation of state at low densities, such as typel (), can include model dependent strong interactions. Altogether our relativistic mean field and Virial EOS models cover the large range of temperatures, densities, and proton fractions shown in Table 1.
Parameter  minimum  maximum  number of points 

[MeV]  0, 10  10  36 
log(n) [fm]  8.0  0.2  83 
Y  0, 0.05  0.56  53 
Note that pure neutron matter (Y=0) is also included in the table. The equation of state at zero temperature is obtained by extrapolation from the two lowest temperatures. There are 73,840 data points from the Virial calculation at low densities, 17,021 data points from the nonuniform Hartree calculation, and 90,478 data points from uniform matter calculations. The overall calculations took 7,000 CPU days in Indiana University’s supercomputer clusters.
It is most efficient for astrophysical simulations of supernovae and neutron star mergers to interpolate within an existing EOS table to update the properties of stellar matter with three independent thermodynamics parameters (, and ). Therefore in this paper we generate, from these original results, a sufficiently large EOS table that can easily be interpolated in a way that preserves thermodynamic consistency.
There exist only two realistic EOS tables that are in widespread use for astrophysical simulations, the LattimerSwesty (LS) equation of state LS (), that uses a compressible liquid drop model with a Skyrme force, and the H. Shen, Toki, Oyamatsu and Sumiyoshi (SS) equation of state Shen98a (); Shen98 (), that uses the ThomasFermi and variational approximations with an RMF model. The two EOS show clear differences for example in the nuclear composition during the corecollapse phase and in the adiabatic index around and above the phase transition to homogeneous nuclear matter MPAreview (). In this paper, we will discuss several comparisons between our new EOS and these two existing ones.
The paper is organized as follows: In section II we present the numerical details for the interpolation scheme we use to generate the full EOS table. Section III shows results for our EOS, including various thermodynamic properties and the composition. The = 0 beta equilibrium EOS is also presented and used to calculate the neutron star massradius relation. We also give some comparisons between our EOS and the two existing EOS tables. Section IV presents a summary of our results and gives an outlook for future work. In the appendix V the contents of our EOS table are explained in more detail and the public access to our EOS table is given.
This paper presents results based on the relatively stiff RMF interaction of Lalazissis, König, and Ring NL3 () that has a high pressure at high densities. In later work we will present additional EOS tables for other softer interactions. This will allow one to more easily run astrophysical simulations for different EOS tables and identify quantities that are sensitive to the EOS. Note that a variety of uncertainties in the EOS and some laboratory measurements and astronomical observations related to the EOS have been discussed in our previous papers SHT10a (); SHT10b ().
Ii Numerical details
In this section we describe how we interpolate our free energy results and calculate other quantities from derivatives of the free energy in a thermodynamically consistent manner.
ii.1 Equation of state at zero temperature
At subMeV temperatures, the EOS is not very sensitive to variation in the temperature. In our calculations, the = 0 EOS is obtained by quadratic extrapolation from the results at the two lowest temperatures = 0.158 and 0.251 MeV Lattimer09 (). Specifically, thermodynamic quantities (such as the free energy) have approximate temperature dependence,
(1) 
and this equation is used to find .
We start with the zero temperature part of the EOS that constitutes the most important contribution to the adiabatic index and the speed of sound, particularly at low temperature. The latter two quantities mostly rely on the derivative , which involves a second derivative of free energy with respect to density and this is sensitive to noise in the bicubic interpolation scheme that we use. We start from the original Table 1, which is used to generate the pressure at zero T (we also include the electron pressure in the following smoothing procedure). The thermodynamic pressure can be obtained numerically from the free energy per baryon ,
(2) 
We smooth the pressure in a 10 points per decade table by removing points that differ significantly from the geometric mean of neighboring points and replace them with interpolated values. Then we use monotonic cubic Hermite interpolation COtt () on the smoothed table of pressure to get a finer table with 40 points pre decade in the density axis. Finally we obtain the (free) energy at zero temperature by integrating this pressure with respect to density.
ii.2 The entropy at finite temperature
Next, we calculate the entropy at finite temperatures MeV. From the free energy points at temperature , density , and proton fraction as in Table 1, the entropy per baryon can be obtained numerically,
(3) 
Finally, the energy per baryon is
(4) 
Note that for higher temperatures 12.5 MeV, matter is uniform for all proton fractions and densities. For uniform matter, all thermodynamic quantities can be obtained directly from relativistic mean field calculations, with good thermodynamic consistency.
The bicubic interpolation scheme we use does not guarantee the monotonicity of the second derivative, e.g.
(5) 
As a result, the internal energy per baryon is not guaranteed to be monotonic either,
(6) 
The bicubic interpolation with slope limiter can guarantee that the free energy is monotonic and that the first law of thermodynamics will be satisfied. But it is not sufficient to preserve monotonicities of the entropy and internal energy, and it can not guarantee smoothness of as well.
The following scheme is used to generate a big table of EOS that is thermodynamically consistent, keeping monotonicities of entropy and energy, and perserving the smoothness of . We start with the entropy in the () plane as in Table 1 with ten points per decade of density and temperature. Then we sweep through the raw table of entropy, discarding numerically noisy points that violate the following constraints:
(7) 
We replace those points by interpolation from neighboring points. Therefore we have a table of entropy that satisfies conditions (7). (The electron and photon parts can be handled more accurately, by monotonic interpolation on individual quantities like energy, pressure and entropy.) Then we further smooth the entropy as a function of density while keeping (7), by further replacing noisy points with interpolated values. This procedure ensures the smoothness of .
Now we perform monotonic cubic Hermite interpolation on this small table of entropy with ten points per decade, to generate a larger EOS table with 40 points per decade in both the temperature and density directions, as indicated in Table 2. Then we integrate this smoothed entropy table as a function of temperature to get values of the free energy (adding the energy at zero temperature to make the full free energy). Thus we obtain the free energy on this finer 40 by 40 points per decade grid in a thermodynamically consistent fashion.
Parameter  minimum  maximum  number of points 

T [MeV]  0, 10  10  110 
log(n) [fm]  8.0  0.2  329 
Y  0, 0.05  0.56  1(Y=0)+52 
ii.3 Bicubic Interpolation of F/A
The final step is to carry out bicubic interpolation of the previous free energy values (as in Table 2) to generate the entropy and pressure by thermodynamic derivatives Eqs. (2,3). This prescription guarantees the monotonicity of entropy and pressure in the final table, and conserves the first law of thermodynamics in adiabatic compression tests.
We first apply bicubic interpolation numericalrecipe () for the free energy. The first derivatives on the grid points are generated from monotonic cubic Hermite interpolation COtt (). The second derivative  the cross derivative , on the grid points is generated as in Ref. numericalrecipe (). The bicubic interpolation can then fit free energies with cubic functions in temperature and density coordinates, and provides the first and second (cross) derivatives. Then the entropy and pressure are otained from Eqs. (2,3). Finally, the boundary points at the highest temperature or density are discarded to avoid boundary artifacts in the interpolation.
Iii Results for the equation of state
In this section we discuss various thermodynamic quantities and the composition in the EOS. First, we discuss the free energy per baryon from Virial gas, nonuniform Hartree mean field, and uniform matter calculations. Next we map out the phase boundaries of our EOS. Then we show various thermodynamic quantities such as the pressure, entropy per baryon, chemical potentials of neutrons, protons, and electron neutrinos, and the adiabatic index. We also show information on the composition, i.e., the average mass number and proton number of heavy nuclei in nonuniform matter, and the mass fractions of different species. The zero temperature EOS in chemical equilibrium is presented next, and this is used to calculate neutron star structure. Finally, we make some comparisons between our EOS and those of LattimerSwesty LS () and H. Shen et al.Shen98 ().
iii.1 Free energy and phase boundaries
In Fig. 1, the free energy per baryon for matter at = 1, 3.16, 6.31 and 10 MeV with different proton fractions is shown as a function of density . The free energy is obtained from Virial gas, nonuniform Hartree mean field, and uniform matter calculations, respectively. In most cases the transition (as the density grows) is found at the density where Hartree calculations or uniform matter calculations give the lowest free energy. For matter at very low temperature (not greater than 1 MeV) and low proton fraction (not greater than 0.1), some matching points are obtained at densities where Virial gas and Hartree calculations give the closest free energies (the difference is less than hundreds of KeV, and note that that mass table FRDM () used in our Virial EOS SHT10b () has itself 600 KeV rms deviations in nuclear binding energies compared to data). Detailed information has been reported in our previous paper SHT10b (). We reproduce the figure here for completeness.
In Fig. 2, we show the phase boundaries of nuclear matter at different proton fractions = 0.05, 0.1, 0.2 and 0.4. The mass fraction of heavy nuclei with mass number is . The boundary at low densities indicates when is greater or less than 10. The boundary at high densities indicates the transition between nonuniform matter and uniform matter. At very low densities, the matter is dominated by free nucleons and alpha particles. As the density rises, heavy nuclei persist to higher temperatures. Finally uniform matter takes over at sufficiently high density. Fig. 2 also shows that as proton fraction rises, the temperature regime with appreciable heavy nuclei grows and the transition density to uniform matter increases. The density for the nonuniform matter to uniform matter transition has a weak temperature dependence. Fig. 3 shows the phase diagram of nuclear matter at different temperatures = 1, 3.16, 6.31, and 10 MeV. Similarly as in Fig. 2, the left boundary at low density indicates where is greater or less than 10. The right boundary at high density indicates the transition between nonuniform matter and uniform matter. As the temperature grows, the dependence of the phase boundaries becomes larger, and the density regime with appreciable heavy nuclei rapidly shrinks. Our Figs. 2 and 3 are very simillar to the corresponding Figs. 2 and 3 of ref. Shen98 () for the H Shen et al. EOS.
iii.2 Pressure
In Fig. 4, the pressure of nuclear matter with different proton fractions is shown as a function of density, for four different temperatures (1, 3.16, 6.31, and 10 MeV). The pressure is obtained via the bicubic Hermite interpolation, in Eq. (2). For matter at low density in the Virial gas phase, the pressure scales nicely with density. When the density increases, formation of additional nuclei decreases the pressure. The Coulomb correction to the lattice energy reduces the pressure even further. Note in part of Virial gas and Hartree WignerSeitz cell regions, the pressure is negative. Figure 4 only shows hadronic contribution to the pressure. In addition one needs to take into account the contribution from electrons and photons, which makes the total pressure (not shown here) positive. We also note that with inclusion of electrons and photons in the final table, the condition for thermodynamic stability against density fluctuations is always preserved to the numerical precision of the table,
(8) 
Finally after the Hartree to uniform matter transition, the pressure rises rapidly with density because of the large incompressibility of nuclear mater.
iii.3 Entropy
In the upper left panel of Fig. 5, the entropy per baryon is shown as a function of density for nuclear matter with = 1 MeV at different proton fractions. The entropy is calculated from Eq. (3). Near a density of a few times 10 fm, the entropy drops rapidly with density, which indicates formation of nuclei, first mostly alpha particles and then heavier nuclei at higher density. Later we will illustrate this by plotting the mass fractions of different species. For nuclear matter with a large proton fraction, larger mass fractions of heavy nuclei are present so that the entropy drops more rapidly with density. When nuclei compose more than 90% of the mass fraction, the entropy remains almost constant as the density increases. The other panels in Fig. 5 show the entropy per baryon in nuclear matter at higher temperatures 3.16, 6.31, and 10 MeV. They have similar characteristics as that for T = 1 MeV, except that the range of densities where the entropy drops rapidly, due to formation of nuclei, becomes smaller because it is harder to form heavy nuclei at higher temperatures.
iii.4 Chemical potentials of , , and .
In Fig. 6, the neutron (left subpanels, a, c, e, g) and proton (right subpanels, b, d, f, h) chemical potentials are shown for nuclear matter with = 0.1 (black curves) and 0.5 (red curves) at different temperatures. For nuclear matter at low density and high temperature, the chemical potentials scale with density. The formation of nuclei breaks the scaling behavior, i.e., the bound nucleons have much smaller chemical potentials.
In Fig. 7, the equilibrium electron neutrino chemical potential, defined as , is shown versus density for nuclear matter at various temperatures and proton fractions. Let us focus on the top left panel where nuclear matter has = 1 MeV. For symmetric nuclear matter, grows rapidly with density, since almost cancels and . For proton rich nuclear matter (), will grow even faster than for symmetric nuclear matter. For neutron rich nuclear matter, has a maximum value at some density. This is due to the competition between the electron chemical potential and the difference of the proton and neutron chemical potentials. The latter difference between proton and neutron is negative when the proton fraction is less than 0.5. The density where has a maximum decreases with decreasing proton fraction. It is 0.2 fm when = 0.3, and 0.05 fm when = 0.05. In the limit of pure neutron matter there is no maximum for , which is just negative by construction.
iii.5 Adiabatic index .
The left panel of Fig. 8 shows the temperature of the adiabat with entropy = 1, for nuclear matter with proton fraction = 0.1, 0.2, 0.3, and 0.4 including baryon, electron, positron, and photon contributions to . At low density, matter is dominated by free nucleons in the gas phase. There are lots of kinetic degrees of freedom and one needs only a small temperature for the matter to have = 1. When the density rises, the formation of heavy nuclei suppresses the entropy and it takes higher temperatures for the matter to reach = 1. The entropy, in this case, comes from a fraction of free nucleons in a gas and the partial occupation of levels in the heavy nuclei. When uniform nuclear matter forms at high density, it takes even a higher temperature to reach the same entropy = 1 because the uniform matter is highly degenerate and the entropy comes only from excitations above the Fermi surfaces.
The adiabatic index ,
(9) 
describes the stiffness of the EOS at constant entropy. In the right panel of Fig. 8, the adiabatic index is shown versus density for nuclear matter with constant entropy = 1 and = 0.1, 0.2, 0.3, and 0.4. At subnuclear density, remains almost constant with small variation. It then rises rapidly at the transition from nonuniform to uniform matter. This characterizes the stiffening of the EOS due to the large nuclear incompressibility at high density.
iii.6 Average and of heavy nuclei.
Figure 9 shows the average mass number of heavy nuclei (with ) versus baryon density at different temperatures and proton fractions. Note that alpha particles are not counted as heavy nuclei. Let us first look at the upper left panel, where = 1 MeV. Nuclear shell effects give rise to several approximate plateaus in vs density for each , for example = 12, 50, 80 and 100. Usually is larger in matter with smaller . There are oscillations in in the Hartree mean field regime. This is due to both nuclear shell effects and small errors in the free energy minimization due to our using a finite step in the Wigner Seitz cell size SHT10a (). The average mass can be as large as 3,000 at high density before the final transition to uniform matter. At higher temperatures, as shown in the other panels of Fig. 9, grows more rapidly to several thousand in a narrower range of density. As a result, the width of the plateaus in with density become much shorter and finally vanish.
The following Fig. 10 shows the average proton number of heavy nuclei versus baryon density at different temperatures and proton fractions. The average proton number has very similar characteristics as for the average mass number, except that the differences in between different s are smaller.
iii.7 Mass fractions of nucleons and different nuclei.
In Fig. 11, the mass fractions of free neutrons, free protons, alpha particles and heavy nuclei are shown versus baryon density at four different temperatures for a proton fraction = 0.05. The matter is extremely neutron rich, so the mass fraction of free neutrons is always appreciable at any density. Note that in the Hartree mean field results there is only a single nucleus associated with each WignerSeitz cell. We define a nucleon level to be free when it has positive energy. The upper left panel is for = 1 MeV. At very low densities  less than a few times fm, the system is basically a free neutron and proton gas. The alpha fraction rises rapidly to 10% around a density of fm. Some of the free neutrons are used to form alphas, so the fraction of free neutrons drops to 90%. The free neutron fraction drops further beyond a density of fm, where heavy nuclei appear with a mass fraction of about 15%. The fraction of heavy nuclei grows faster after the Virial gas  Hartree mean field transition around 3 fm. More and more free neutrons are bound to form heavier nuclei, as can be seen in Fig. 9. The fraction of heavy nuclei reaches a maximum value of 70% before the final transition to uniform matter. The other panels give the mass fractions of different species at higher temperatures. In these cases, the appearance of alpha particles and heavy nuclei occur at higher densities and with smaller mass fractions compared to the = 1 MeV case. When = 6.31 MeV, there is only a narrow range of densities where alpha particles have a mass fraction of a couple of percent and the amount of heavy nuclei is very small.
Figs. 12 and 13 show the mass fractions as were shown in Fig. 11, but for matter with higher proton fractions = 0.2 and 0.5, respectively. Alpha particles appear at similar densities as for = 0.05, but with much larger mass fractions. For example, the alpha fraction can almost reach 100% around fm when = 1 MeV and = 0.5. Heavy nuclei also appear at similar densities as in = 0.05, and with much larger mass fractions. As the proton fraction rises, heavy nuclei can survive to higher temperatures, even though only at very high density. For example, when = 10 MeV and = 0.5, heavy nuclei have a large mass fraction for densities around a few times 10 fm.
iii.8 Beta equilibrium at zero temperature.
Matter in beta equilibrium at zero temperature is particularly interesting and determines neutron star structure, including the massradius relation and the neutron star crust. The = 0 EOS is included in our table 2, which can be used to find the beta equilibrium EOS by imposing the constraint . The resulting energy per baryon, pressure, and proton fraction are shown as function of density in Fig. 14. The energy per baryon starts near 8.5 MeV, a typical value for stable nuclei, at 10 fm, becomes positive around 10 fm, and finally rises rapidly after the transition to bulk uniform matter around 0.05 fm. The pressure vs density is polytropic like below 10 fm, and then it develops small plateaus. After the transition to uniform matter, the pressure rises more rapidly which indicates a stiffening of the EOS. At 10 fm matter is made up of stable nuclei with similar numbers of neutrons and protons. That’s why the proton fraction at 10 fm is as high as 0.45. The proton fraction drops as the density increases, which means matter becomes more and more neutron rich due to the large electron Fermi energy. The proton fraction reaches a minimum value of 0.007 around 0.025 fm. However, it increases again due to the large symmetry energy at high density for the NL3 NL3 () interaction, which favors more nearly equal neutron and proton fractions.
iii.9 Neutron star structure
The mass and radius of a neutron star are obtained by solving the TolmanOppenheimerVolkoff (TOV) equations oppenheimer39 (); tolman39 ().
By using the EOS in beta equilibrium at zero temperature from the previous section in the core and using BaymPethickSutherland EOS in the crust BPS (), we obtain the neutron star massradius and masscentral density relations, which are plotted in Fig. 15. The neutron star’s maximum gravitational mass is about 2.77 solar mass with a radius of 13.3 km. The corresponding central density is about 10 g/cm. The large value for the maximum neutron star mass is not too surprising since NL3 provides a very stiff EOS as we have emphasized several times.
iii.10 Comparisons with LattimerSwesty’s EOS and H. Shen et al.’s EOS
It is instructive to compare our EOS with the LattimerSwesty equation of state, that uses a simple liquid drop model with a nonrelativistic Skyrme interaction, and the H. Shen et al.equation of state, that uses the Thomas Fermi and variational approximations to a relativistic mean field model. The LS EOS quoted in this work corresponds to the one with incompressibility coefficient =180 MeV. In Fig. 16, the pressure for matter at = 1, or 6.31 MeV and = 0.05 or 0.4 is shown for our EOS, LattimerSwesty’s (LS) and H. Shen et al.’s (SS) EOSs. Here the pressure includes contributions from electrons, positrons and photons. Below nuclear saturation density 0.16 fm, ours and SS EOS agree very well. The LS EOS gives a slightly lower pressure at densities above 10 fm when = 1 MeV and = 0.05. Above saturation density, our EOS gives the largest pressure.
Fig. 17 compares the mass fraction of heavy nuclei from our EOS, LS EOS and SS EOS, for matter at = 1 or 6.31 MeV, and = 0.05 or 0.4. Although the pressure agrees well among the three EOSs, the mass fraction of heavy nuclei and alpha particles can be different and this contributes to the differences in entropy shown in Fig. 18. In addition our EOS may have a higher entropy in Fig.18 (b) because we may have smaller nuclei for less than fm compared to SS and because we can have a distribution of several heavy nuclei, while the LS and SS EOSs use only a single average species. The average and of heavy nuclei is shown in Fig. 19. This figure clearly shows the effects of shell structure, that is included in our EOS but is neglected in both LS and SS EOSs. As a result, and for our EOS have a serious of steps while and for LS and SS EOSs increase smoothly with density.
For = 1 MeV and = 0.05, the three EOSs agree very well on entropy, except that at low densities the LS EOS gives a slightly smaller value. It is comprehensible since the matter is dominated by free neutrons until very high density when heavy nuclei dominate (upper left panel of Fig. 17). But heavy nuclei have low entropy at low temperatures so the difference in entropy between the EOSs is small. For = 1 MeV and = 0.4, the differences in entropy between the three EOSs are larger. This is partially due to the smaller scale used. It is also because the different EOS predict different mass fractions of heavy nuclei, whose formation decreases the entropy. The LS EOS has 80% heavy nuclei at 10 fm (upper right panel in Fig. 17) compared to 20% in our EOS and the SS EOS, therefore the LS EOS has the lowest entropy. At the higher temperature = 6.31 MeV, the three EOSs agree well on entropy for different values of , because free nucleons dominate at most densities (see lower panels in Fig. 17).
The differences in abundances and average and between EOSs may arise because of differences in symmetry energy and approximations made. For example LS uses a very simple liquid drop model while the SS EOS uses variational and Thomas Fermi approximations that are avoided in our EOS. These differences could be important for neutrino interactions and should be studied further.
iii.11 Adabatic compression tests
One important test of thermodynamic consistency for our EOS table is to check that entropy is conserved as matter undergoes adiabatic compression. This is very closely related to the conservation of energy in the first law of thermodynamics. Let us consider the following adiabatic compression procudure.

Start from an energy per baryon , a density and a temperature . Here stands for a sequence of times starting at . The pressure is then given by the EOS.

Now compress the matter by increasing the density by an amount so that with .

Next ensure that the first law of thermodynamics is satisfied by using a backward difference equation and requiring
(10) by solving for the new pressure and temperature .

Finally the updated entropy follows from our EOS at the new density and temperature .
The entropy should be conserved in this adiabatic compression so that is independent of . In Fig. 20, the temperature and entropy versus density during adiabatic compression are shown for nuclear matter with fixed proton fraction 0.15. The initial density is 5.16738 10 fm and initial temperature is 0.5 MeV. The solid curve is test result for our EOS table with lepton and photon contributions included in addition to contributions from baryons. The dashed curve is a test result for our EOS table with baryons alone. The variation in entropy for both cases is less than 1%. For a comparison the test result for H. Shen EOS is also shown (as dotdashed curve) for the similar initial condition. Note that it is important to use an accurate interpolation scheme with the EOS table to ensure that the first law is satisfied so that entropy is conserved. The adiabatic compression test result for H. Shen EOS was obtained using the routine developed in Ref. Ott09 ().
Iv Summary and Outlook
In this paper we generate a new complete equation of state for use in supernova and neutron star merger simulations. We use a density dependent RMF model for nuclear matter at intermediate and high density with a spherical WignerSeitz approximation for nonuniform matter, which incorporates nuclear shell effects SHT10a (). For nuclear matter at low density SHT10b (), we use a Virial expansion for a nonideal gas consisting of neutrons, protons, alpha particles and thousands of heavy nuclei from the finite range droplet model (FRDM) mass table FRDM (). We include second order virial corrections for light elements 4, nuclear partition functions for heavy nuclei, and Coulomb corrections between electrons and heavy nuclei. As the density decreases, the mean field results match smoothly to the Virial gas. As the density goes down, the Virial expansion reduces to nuclear statistical equilibrium and the Virial expansion is exact in the low density limit. The computed EOS covers 180,000 grid points in the temperature range = 0 80 MeV, the density range = 10 1.6 fm, and the proton fraction range = 0 0.56.
We use a hybrid interpolation scheme to generate a full EOS table on a fine grid that is thermodynamic consistent. This ensures that the first law of thermodynamics is satisfied and that entropy is conserved during adiabatic compression. Our EOS is an improvement over the existing LattimerSwesty LS () and H. Shen et al. Shen98a (); Shen98 (), equations of state because our EOS includes thousands of heavy nuclei and is exact in the low density limit.
We discuss the thermodynamic properties of our EOS in detail. The free energy per baryon, pressure, entropy per baryon and chemical potentials for neutrons, protons, and electron neutrinos are shown for matter at high and low temperatures and for large and small proton fractions. We also perform adiabatic compression tests for our EOS table and show that it preserves entropy to better than 1%.
Finally, we show some comparisons between our EOS and two existing EOS tables. At high density, our RMF model reduces to the normal NL3 set. It produces a stiffer EOS at high density than the SS EOS, which is already stiffer than the LS EOS. One can find this difference, for example, in Fig. 16 for the pressure at high density. As a result, neutron stars obtained from our EOS have a large maximum mass of around 2.77 solar mass and radius of 13.3 km. At low density, our model is essentially the virial expansion of a nonideal gas, which includes many nuclei. Our model usually has a broad distribution of nuclei as shown in our previous paper SHT10b (). This spread in the distribution of nuclei makes the composition of nuclear matter in our EOS different from that in the LS and SS EOS tables. This composition difference also introduces some differences in the entropy for nuclear matter at low temperature and large proton fraction, see for example the upper right panel in Fig. 18 for matter with = 1 MeV and = 0.4. More importantly, the composition of nuclear matter, especially the wide spreads in distribution of nuclei, could influence the position of the neutrinosphere and the spectra of radiated neutrinos and antineutrinos.
The complete equation of state table is not only useful for core collapse supernova simulations, but also for other astrophysical applications, including the long term evolution of protoneutron stars, and neutron star mergers. Our procedure, for matching mean field and virial calculations to produce a thermodynamically consistent EOS, can be used for a variety of effective interactions. In this paper we have presented an EOS based on the NL3 interaction that is relatively stiff at high densities. In future work we will present additional EOS tables based on softer interactions that have lower pressures at high densities.
Finally, our full EOS tables, both with and without lepton and photon contributions, are available for download as described in Appendix V.
We thank Lorenz Hedepohl, Thomas Janka, Andreas Marek, Evan O’Connor, and Christian Ott for important help running astrophysical simulations to debug our equation of state. This work was supported in part by DOE grant DEFG0287ER40365.
V Appendix: format of EOS tables.
Here we describe the physical quantities provided in the tables for the equation of state, NL3eos1.03.dat and NL3eosb1.03.dat and where they can be downloaded. One should download the gzip compressed files NL3eos1.03.dat.gz and or NL3eosb1.03.dat.gz (that are about 100 MB each) and use gunzip to decompress them. The grid structures of these tables are indicated in Table 2 and contain approximately 517 MB of data each. The tables, a sample FORTRAN computer program, and a readme file are available for download both at our website http://cecelia.physics.indiana.edu/gang_shen_eos/ and, from the Electronic Physics Auxiliary Publication Service  EPAPS web site, EPAPS (). Please check our web site for any updated information regarding these EOS tables.
The tables contain a number of quantities. For a single triplet of , , and , there are 16 items that make up a row of the table. These items are described below. In the table NL3eosb1.03.dat, only the contribution from baryons is taken into account for items 4,5,6. In the table NL3eos1.03.dat, the contributions from electrons, positrons, and photons are also included. The electron mass is 0.511 MeV.
v.1 Items in EOS tables.

Temperature [MeV]. The range of temperature is first =0 then from MeV, see Table 2.

Proton fraction . The range of proton fraction is first 0, and then from, 0.05 0.56. The step in proton fraction is 0.01, see Table 2.

Baryon number density [fm]. The range of density is from 10 to 10 fm. See Table 2.

Free energy per baryon [MeV] which has subtracted the free nucleon mass 939 MeV.

Pressure [MeV/fm].

Entropy per baryon [k].

Chemical potential for neutrons [MeV]. The tabulated value is relative to the nucleon mass 939 MeV.

Chemical potential for protons [MeV]. The tabulated value is relative to the nucleon mass.

Chemical potential for electrons [MeV]. The tabulated value includes the electron mass.

Average mass number of heavy nuclei with , which exclude alpha particles.

Average proton number of heavy nuclei with , which exclude alpha particles.

Mass fraction of free neutrons.

Mass fraction of free protons.

Mass fraction of alpha particles.

Mass fraction of heavy nuclei with , which exclude alpha particles.

Effective nucleon mass [MeV]. In uniform matter it is obtained from RMF theory. For virial gas and nonuniform matter, it is chosen to be the free nucleon mass MeV.
v.2 Sample FORTRAN computer program readeos.f
The Fortran program readeos.f includes a very short main program that calls the subroutine load_table, to read NL3eos1.03.dat or NL3eosb1.03.dat, and then calls the subroutine readeos with inputs (in MeV), proton fraction , and density (in fm). The subroutine readeos uses triliner interpolation (in , , and ) to return the above 16 values plus the internal energy per baryon (in MeV) and the chemical potential for electron neutrinos in chemical equilibrium (in MeV). Note that one needs to call load_table only once and then one can call readeos many times. For further details please see the comments in readeos.f.
Footnotes
 email: gshen@indiana.edu
 email: horowit@indiana.edu
 email: steige@indiana.edu
References
 H.Th. Janka, K. Langanke, A. Marek, G. MartínezPinedo, and B. Müller, Phys. Rep. 442, 38 (2007).
 A. Arcones, G. MartínezPinedo, E. O’Connor, A. Schwenk, H.T. Janka, C. J. Horowitz, and K. Langanke, Phys. Rev. C 78, 015806 (2008).
 See for a recent review, H. Duan, G. M. Fuller, and Y.Z. Qian, Ann. Rev. Nucl. Part. Sci. 60, 569 (2010).
 G. Shen, C. J. Horowitz and S. Teige, Phys. Rev. C 82, 015806 (2010).
 G. Shen, C. J. Horowitz and S. Teige, Phys. Rev. C 82, 045802 (2010).
 C.J. Horowitz, A. Schwenk, Nucl. Phys. A776, 55 (2006).
 C.J. Horowitz, A. Schwenk, Phys. Lett. B638, 153 (2006).
 C.J. Horowitz, A. Schwenk, Phys. Lett. B642, 326 (2006).
 C.J. Horowitz, Phys. Rev. D 55, 4577 (1997).
 S. Typel, G. Röpke, T. Klähn, D. Blaschke, and H. H. Wolter, Phys. Rev. C 81, 015803 (2010).
 J. M. Lattimer and F. D. Swesty, Nucl. Phys. A 535, 331 (1991).
 H. Shen, H. Toki, K. Oyamatsu and K. Sumiyoshi, Nucl. Phys. A 637, 435 (1998).
 H. Shen, H. Toki, K. Oyamatsu and K. Sumiyoshi, Prog. Theo. Phys. 100, 1013 (1998).
 G. A. Lalazissis, J. König, and P. Ring, Phys. Rev. C 55, 540 (1997).
 We thank C. Ott for providing us the subroutine for 1D cubic Hermite interpolation with slope limiter.
 W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes, pp 95101, Cambridge University Press 1989.
 Jim Lattimer (private communication, 2009).
 P. Möller, J. R. Nix, and W. J. Swiatecki, Atom. Data Nucl. Data Tables, 59, 185 (1995).
 J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374 (1939).
 R. C. Tolman, Phys. Rev. 55, 364 (1939).
 G. Baym, C. Pethick, and P. Sutherland, Ap. J. 170, 299 (1971).
 E. O’Connor and C. D. Ott, Class. Quant. Grav. 27, 114103 (2010).
 See supplemental material at [http://link.aps.org/supplemental/10.1103/PhysRevC.83.035802].