Equation of State of nuclear matter in a Virial expansion of nucleons and nuclei
Abstract
We study the equation of state (EOS) of nuclear matter at subnuclear density in a Virial expansion for a nonideal gas. The gas consists of neutrons, protons, alpha particles, and 8980 species of nuclei with and masses from the finiterange droplet model (FRDM). At very low density, the Virial expansion reduces to nuclear statistical equilibrium. At higher density, the Virial results match smoothly to the relativistic mean field results discussed in our previous paper. We tabulate the resulting EOS at over 73,000 grid points in the temperature range = 0.158 to 15.8 MeV, the density range = 10 to 0.1 fm, and the proton fraction range = 0.05 to 0.56. In the future we plan to match these low density results to our earlier high density mean field results, and generate a full EOS table for use in supernova and neutron star merger simulations. This Virial EOS is exact in the low density limit.
pacs:
21.65.Mn,26.50.+x,26.60.Kp,51.30.+i,97.60.BwI Introduction
One of the main ingredients in simulations simulation1 ; simulation2 of core collpase supernovae and neutron star mergers is the Equation of State (EOS) for hot dense nuclear matter. The EOS, along with detailed information on the composition of nuclear matter, may play an important role in the neutrinomatter dynamics Burrows04 and supernova explosion mechanism. In a previous paper SHT10a we used a relativistic mean field (RMF) model to calculate the freeenergy of nonuniform matter at intermediate density and uniform matter at high density.
Simulations of the core collapse and explosion of supernovae depend heavily on the EOS at subnuclear density, especially the detailed composition. In this work, we study subnuclear density nuclear matter in a Virial expansion for a nonideal gas, consisting of neutrons, protons, alpha particles, and 8980 species of heavy nuclei () with masses from the finiterange droplet model (FRDM) FRDM . The Virial results will cover the density range = 10 to 0.1 fm, the temperature range = 0.158 to 15.8 MeV, and the proton fraction range = 0.05 to 0.56. (For temperature higher than 15.8 MeV, matter is uniform and fully described in the RMF model SHT10a .) The distribution of nuclei for given conditions is obtained in this approach, while the existing EOS tables of LattimerSwesty (LS) LS and H. Shen, Toki, Oyamatsu and Sumiyoshi (SS) Shen98a ; Shen98 , use a single heavy nucleus approximation. Our Virial EOS will be matched, using a thermodynamically consistent interpolation scheme, to the RMF EOS obtained in our previous paper SHT10a to generate a full EOS table for supernova.
Detailed information on the distributions of nuclei in the EOS table is important for neutrinomatter dynamics. Neutrinos radiate 99% of the energy released in supernovae. Besides gravitational wave signals, neutrinos are the only messenger through which one can directly probe the EOS inside supernovae. The neutrinosphere is the surface of last scattering before or escape. The neutrinosphere is expected to occur at a density around g/cm and this is consistent with the available information from a handful events in SN1987a SN1987a . The composition of matter at subnuclear density constrains the position of the neutrinosphere and influences the spectra of emitted neutrinos and antineutrinos. For example, in a recent study light07 , light nuclei with mass 2, 3 and 4 were found to have an important influence on the spectra of antielectron neutrinos.
For matter at subnuclear density and low entropy, heavy nuclei tend to form. Nuclear statistical equilibrium (NSE) models treat low density nuclear matter as a system of noninteracting nuclei in statistical equilibrium, taking into account the binding properties of heavy nuclei. This has been widely used in nuclear astrophysics NSE . Recently, there have been several NSE based studies of the supernova EOS, see for example Refs. NSE09 and Hempel09 . These studies use modern mass tables with thousands of nuclei and include excluded volume effects Hempel09 . The NSE models have the advantage of generating thermodynamically consistent EOS tables. However they also have several disadvantages. First, NSE models themselves can not be used to describe well the nonuniform matter at nearly nuclear density, when exotic pasta phases may appear. To generate a complete EOS table, NSE models need to be matched to uniform matter models at high density. Thus low and high density matter are usually described with different models. Second, NSE models do not provide a systematic way to include the strong interactions between nuclei. These interactions could be important as the density increases. Moreover, neutron matter at low density is closely analogous to a unitary gas, whose properties can not be described satisfactorily in either NSE models or normal mean field approaches (see for example unitarygas ). The Virial expansion, on the other hand, has been successfully used to describe the properties of unitary gases Ho04 as well as neutron matter at low density HS05b .
In a previous work HS05 on low density nuclear matter consisting of neutrons, protons and alpha () particles, the Virial expansion has been used to systematically incorporate second virial corrections from nucleonnucleon (NN), N and elastic scattering. These light nuclei components are expected to dominate in matter with high entropy. The effects of second virial corrections are found to be important for the composition and for the equation of state. In this work we use the Virial expansion for a nonideal gas consisting of neutrons, protons, alpha particles, and thousands of heavy nuclei. The heavy nuclei included in this work are from FRDM mass tables FRDM , which have 8979 nuclei with 16. We also include C in the mass table. As mentioned above, the mass 2 and 3 nuclei, and other nuclei asside from He may be important for the spectra of antielectron neutrinos. The inclusion of these light elements will require information on second virial corrections among them and nucleons, which are under current investigation mass3 . In this work the light elements are represented by alpha particles, which typically have the largest abundance among the light elements in statistical equilibrium. At very low density, the Virial expansion agrees with nuclear statistical equilibrium models. As the density rises, we find that the Virial gas matches smoothly to previous mean field results.
In this work we include the second order virial corrections among nucleons and alphas as in Ref. HS05 . Partition functions for heavy nuclei are included using the recipes of Fowler et al.Fowler78 . In addition, some calculations are presented using partition functions based on the recipe of Rauscher et al.Rauscher97 for comparison. Equivalently, the partition functions can be considered as the sum of successive high orders of the Virial expansion for heavy nuclei. There have been many studies on the level density and partition functions of hot nuclei in astrophysical environments. For large scale astrophysical applications, it is necessary to find both reliable and computationally practical methods for the level density. Most of these studies Gilbert65 ; Fowler78 ; Mazurek79 ; Engelbrecht91 ; Rauscher97 followed the original noninteracting Fermi gas model of Hans Bethe Bethe36 . For astrophysical nuclear reactions with temperature below a few times 10 Rauscher97 ; Rauscher00 , this phenomenological approach gave excellent agreement with more sophisticated Monte Carlo shell model calculations mc , as well as combinatorial approaches Engelbrecht91 ; Paar97 . This justifies the application of the Fermigas description at and above the neutron separation energy. For temperatures higher than a few times 10, there are big ambiguities in the values of the partition functions. However, as suggested by some authors NSE09 and supported by our own calculations, these uncertainties have only a small impact on the thermodynamics of dense matter.
The effects of Coulomb interactions in the plasma can be estimated by the plasma parameter , where is the atomic number of the nucleus, is the temperature and is the spacing between nuclei. For matter at low density, is smaller than one and the effect of Coulomb corrections is small. However for matter at higher density and when the dominant species carry large charges, can be much greater than one and the effect of Coulomb interactions should be taken into account. The Coulomb correction to the plasma has been studied analytically up to high by the cluster expansion coulomb . Generally the correction due to electronion interactions will reduce the free energy of the plasma and eventually crystalize the matter at high density. For simplicity, in this work the Coulomb interactions between nuclei and electrons are included via a WignerSeitz approximation with effective ion spheres for each species of nuclei, wherein local electrical neutrality is maintained. This WignerSeitz approximation for the Coulomb correction will be compared with a more rigorous cluster expansion method.
Based on the above Virial expansion, we generate an equation of state table which covers the range of temperatures, densities and proton fractions shown in Tab. 1.
Parameter  minimum  maximum  number of grids 

T [MeV]  10  10  20 
log(n) [fm]  8.0  1.0  71 
Y  0.05  0.56  52 
There are 7 points in temperature (0.16, 0.26, 0.40, 0.63, 0.71, 0.79 and 0.89 MeV) for temperature below 1 MeV. For higher temperatures we use a spacing of 0.1 in . This is a total of 20 different temperatures from = 0.16 to 15.8 MeV. We use a spacing of 0.1 in , giving a total of 71 points in density for 0.1 fm. We use a linear step 0.01 in proton fraction, giving a total of 52 points in proton fraction for Y = 0.05 to 0.56. There is a total of 73,840 data points in the Virial gas calculation. This took 1,000 CPU days in Indiana University’s supercomputer clusters.
The paper is organized as follows: in section II our Virial expansion for a nonideal gas is explained in detail. In section III we present the recipes for the nuclear partition function used in the Virial expansion. Section IV discusses the effect of different mass tables on the EOS, and shows several examples for the free energy and the distribution of nuclei in the Virial EOS. Section V presents a summary of our results and gives an outlook for future work.
Ii Formalism
We now describe our Virial expansion formalism for a gas consisting of neutrons (n), protons (p), alpha particles and heavy nuclei. The grand partition function for a gas at pressure and volume is expanded to second order in the neutron fugacity , the proton fugacity and the alpha particle fugacity as follows,
(3)  
where is the partition function for nuclei and and are the second Virial coefficients as defined in Ref. HS05 . The sum on runs over different heavy nuclei, for which we use the FRDM mass table FRDM for A 16 and include C. The thermal wavelength for species is ,
(4) 
From now on are used for sums over heavy nuclei and are used for sums over all species.
There exist several recipes for the nuclear partition function . We will use that of Fowler et al.Fowler78 in this work. We also consider the choice of Rauscher et al.Rauscher97 as an option. Different choices of partition functions change the matching densities to our mean field results slightly, but the influence on the thermodynamics is negligible. This is also discussed in Ref. NSE09 .
Chemical equilibrium between nucleons and a heavy nucleus with protons and neutrons insures,
(5) 
where are chemical potentials of the heavy nucleus, protons and neutrons, respectively. Therefore the fugacity of a heavy nucleus is readily obtained,
(6) 
where is the binding energy of heavy nucleus .
We consider the Coulomb interaction between the electron background and nucleus following Baym et al.BBP , but generalized to multiple species of nuclei. The total Coulomb energy of a nucleus in an electron background is,
(7) 
where fm is the nuclear radius (in accordance with that in FRDM FRDM ), and is the average ion sphere radius defined through
(8) 
We emphasize that the sum over runs over heavy nuclei. This definition ensures local charge neutrality inside each ion sphere. The first term in Eq. (7) comes from the protonproton Coulomb energy in the nucleus and is also included in the binding energy of the nucleus. When approaches the total Coulomb energy in a heavy nucleus vanishes as expected. So the Coulomb correction to the binding energy of nucleus is
(9) 
Adding this correction to the binding energy of a nucleus, Eq. (6) becomes
(10) 
The densities of each species can be obtained from
(11) 
This gives,
(12)  
(13)  
(14)  
(15) 
The mass fraction of species is defined as,
(16) 
The two conditions used to determine the fugacities of the neutrons and protons are that the total baryon density is conserved,
(17) 
and that the proton fraction is given by,
(18) 
Since the Coulomb correction is included, one extra loop is needed to selfconsistently determine the values of the ion sphere radii for each species.
The entropy density of the Virial gas is obtained from,
(22)  
where the prime indicates a derivative with respect to temperature. Note that the last term is the Coulomb correction to the entropy, which is nontrivial to evaluate directly. However the free energy of the Virial gas can be obtained directly (see discussion below) so that we can evaluate the entropy from the derivative of the free energy.
The energy density can be calculated from the entropy density,
(23) 
where the sum on runs over neutrons, protons, alphas and heavy nuclei.
The free energy density is given by
(24)  
(26)  
From the above equation we define an effective Coulomb correction to free energy per nucleon ,
(27) 
The free energy per nucleon is
(28) 
Iii Numerical details
In Subsection III.1 we describe some details of the evaluation of the partition functions and then in Subsection III.2 we briefly describe the parallel computations of the EOS for many different temperatures, densities, and proton fractions.
iii.1 Recipes for partition function
Fowler et al.Fowler78 proposed an efficient approximation for the partition function of hot nuclei, which has a closed form,
(31) 
where is the contribution from known discrete states and is the continuum subtraction. One could use where is the ground state spin. Inaccuracies from this approximation will become progressively unimportant beyond 10 K. The continuum subtraction, on the other hand, only becomes important beyond 10 K (see also Nadyozhin04 ), by this temperature uniform matter will be more stable in most regions of phase space, as will be shown in the following discussions. So in the temperature range we are interested, one can discard the latter term. Therefore the nuclear partition function becomes,
(32) 
and its derivative versus temperature is,
(33) 
A widely used expression for the level density is the backshifted Fermi gas formula Bethe36 ,
(34) 
where is the level spacing, and with a backshift parameter related to pairing. Various prescriptions of these two parameters are available, which reproduce more rigorous results for the level density from Monte Carlo or combinatorial calculations. We will use the prescription from Fowler et al Fowler78 in our calculation and include that of Rauscher et al.Rauscher97 as an option in our code.
The integral limits are determined by the following equations Fowler78 :
(35)  
(36) 
and are the single neutron and single proton separation energies, which could be obtained from the mass table. is the zeropoint energy, with . The Coulomb barrier is . For unstable nuclei ( or 0), the partition function is set to the ground state value . There is still a large ambiguity regarding the value of the integral upper limit , which is related to the contribution of the continuum. This ambiguity will introduce big changes in the partition function only when the temperature is above several MeV (few times 10). We compared results with different choices of and found the ambiguity changes the transition density to our relativistic mean field EOS SHT10a slightly, but the influence on the thermodynamics is negligible. This was also discussed in a previous work NSE09 .
When is larger than , should be used. In this case Eq. (34) is inapplicable due to the pole in the lower integration limit. A constant temperature formula,
(37) 
is used instead. The constants and are obtained by the continuity of the level density and its first derivative when matched to Eq. (34).
iii.1.1 Choices of a and by Fowler et al.. Fowler78
In our calculation we use the prescription of Fowler et al.Fowler78 for the level spacing parameter (in MeV) and pairing parameter (in MeV):
(38) 
where . When , is set to , and Eq. (37) is used below an energy 2. By continuity of the level density and its derivative, the constants and are obtained,
(39)  
(40) 
In our calculation is chosen to be .
iii.1.2 Choices of a and by Rauscher et al.. Rauscher97
In the parameterization of Rauscher et al., the level spacing depends on energy,
(41) 
where describes properties of a nucleus differing from the spherical macroscopic energy, which is already given in the FRDM mass table. Based on the FRDM mass formula, the best fitted values for , and are obtained by comparing the level densities with experimental analysis: = 0.1337, = 0.06571, = 0.04884 Rauscher97 .
The pairing parameter is given by,
(42) 
where the neutron pairing energy is
(43) 
and is the binding energy of nucleus . The proton pairing energy is calculated in a similar way. The constants and are obtained as following,
(44)  
(45) 
where
(46) 
iii.2 Computational methodology
We calculate, in parallel, the points in parameter space covered by the Virial EOS, analogous to the way used in our previous paper to perform relativistic mean field calculationsSHT10a . There are a total of 73,840 points in the 3dimensional parameter space of temperature , proton fraction , and density . Each point takes up to 20 minutes to calculate. The overall Virial EOS took about 1,000 CPU days in Indiana University’s supercomputer cluster.
Each point in the parameter space was mapped to a unique integer that we refer to as the job index. A file, runlist, was prepared with a list of job indicies for the whole parameter space, and a single character (A=available, R=running, r=Rerunning, C=complete, T=timelimited and F=failed) that gives the status of calculations for that job index. A Message Passing Interface (MPI) parallel wrapper code manages the running of the many requested tasks. Typically, one parallel job requests a set of compute cores (usually 256). Each MPI rank, using a single CPU core, is assigned one job index corresponding to one point in the parameter space and evaluates the required quantities.
Initially, rank zero of the MPI job

locks the job listing file runlist,

reads runlist until a list of available tasks is filled,

closes runlist and releases the lock,

passes a job index to each MPI rank and begins the calculation for that job index.
When the calculation completes (or timelimits or fails) for a given MPI rank, the status character for the job index in runlist is modified appropriately. The now available MPI rank will search runlist for the next available task and the calculation restarts for the new job index. Since completion occurs asynchronously file locking is not used for this part of the process.
A simple batch job runs through the points in parameter space. A wall clock limit (48 hours) is used. Each rank of the MPI job can run a series of points using the above procedure, efficiently using each available core for the requested wall clock period. One job per core is running when the wall clock limit is reached. These jobs are identified by being left in the ”R” state after the batch job completes. This procedure allows us to calculate 99 of the points in the runlist file.
Iv Results
In this section we discuss the Virial EOS in detail. First we show the effect of different mass tables on the EOS and its influence on the matching to RMF results. We also discuss the coulomb correction and compare our results with some analytical cluster expansion analysis. Second we discuss the matching between the Virial EOS and RMF EOS for several choices of temperature and proton fraction. Finally we show some examples of the distribution of nuclei obtained in the Virial EOS. We also compare some examples of mass fractions of nuclear matter between Virial EOS and existing EOS tables.
iv.1 Effect of different mass tables
We use both the FRDM FRDM and HFB14 HFB14 mass tables in the Virial expansion. Here we compare the free energy and average charge number (average mass number = /) from the Virial expansion with the two different mass tables. We also match them to the relativistic mean field results from our previous paper SHT10a .
iv.1.1 Free energy
In Figure 1, free energies are shown as a function of density for nuclear matter at = 1 and 3.16 MeV and = 0.3. Virial gas results are obtained from Eq. (28) for the FRDM (black solid curve) and HFB14 (red dashed curve) mass tables. Nonuniform (or Hartree) mean field calculations (circles) and uniform matter calculations (dotted dashed curve) are also shown (See Eqs. (10) and (18) in Ref. SHT10a ). These calculations give lower free energies at high densities. The two different mass tables in the Virial expansion give very similar free energies for nuclear matter. The difference between the black curve (FRDM) and red curve (HFB14) is very small until densities where the Hartree mean field results are more stable than either mass formula. For these two cases, the transition from Virial to Hartree EOS occurs at a density of about 3.9810 fm.
iv.1.2 Average and
Similar to Fig. 1, the lower panels in Fig. 2 give the corresponding average charge number from the Virial EOS with either FRDM or HFB14 mass tables, and from Hartree mean field results. Again the Virial EOS with the two mass tables predict very similar values for . Moreover, the Virial EOS with either mass table gives very similar to that from the Hartree mean field results at the transition density 3.9810 fm (blue dotted line). The fluctuation of in Hartree results below the transition density (at = 3.16 MeV) is probably due to finite step error in search of cell size. In conclusion, the composition, free energy, and transition density to Hartree mean field results, depend little on the choice of mass tables.
In the two upper panels in Fig. 2, we compare the Coulomb energy correction in our Virial expansion, Eq. (27), with an analytic cluster expansion for the onecomponent plasma, Eq. (22) in coulomb . The overall agreement for densities below the VirialHartree transition is good, though at higher temperature the differences become larger. The reason is probably that we calculate the multicomponent contribution to the Coulomb correction in the Virial gas while we only used the average charge number in the analytic formula for the onecomponent plasma.
iv.2 Free energy and phase boundaries
We show in Fig. 3 the transition densities between the Virial and RMF EOSs. We also show the free energy per nucleon as a function of density for = 1, 3.16, 6.31 and 10 MeV. At low densities, is obtained from Eq. (28) in the Virial expansion. The free energy per nucleon is also shown for Hartree mean field calculations at intermediate densities, and from uniform matter at high densities. The latter two have been obtained in our first paper SHT10a .
In most cases the transition (as density grows) is found at the density when the Hartree or uniform matter calculation gives a lower free energy compared to Virial results. 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 the density when the Virial gas and Hartree calculation give the closest free energy. The difference is below hundreds of KeV and is comparable to the mean deviation of binding energy for nuclei ( 600 KeV) in the FRDM mass table FRDM .
In each panel in Fig. 3, the vertical red dashed curves give the Virial gas  Hartree WignerSeitz cell transition densities, and the red solid curves give the transition densities to uniform matter. As temperature increases, the second transiton may happen at lower density than that of first transiton, which means the Hartree WignerSeitz cell is energetically unfavorable for all densities. This indicates the critical temperature for the nucleon liquid  Virial gas transition. We note here that the Virial gas may still contain a small amount of alpha particles and heavier nuclei even above this transition temperature. As seen from the figures, for matter at any temperature and proton fraction the Virial  mean field transition densities are larger than a few times 10 fm, which is about the neutrinosphere density. So in almost all regions around and below the neutrinosphere density, our EOS is represented by the Virial results, which include multiple nuclei.
For matter at low density in the Virial gas phase, the free energy scales nicely with density. The reason is as follows. Low density matter is dominated by neutrons and protons, and the interactions among them are not important due to the large particle spacing, so the kinetic energy dominates and scales as the logarithmic value of the density. As the temperature rises, the scaling behavior extends to higher densities until heavy nuclei appear, and electromagnetic and strong interaction corrections become important. Formation of nuclei greatly decreases the free energy, both in the Virial gas and Hartree mean field regimes.
iv.3 Mass fractions of species
The Virial expansion gives the distribution of heavy nuclei (refer to Eq.(15, 16)), where 8980 species of nuclei are in thermal and chemical equilibrium with free neutrons, free protons and alpha particles. This is an improvement over the LS EOS and SS EOS that both used a singlenucleus representation.
Figure 4 shows mass fractions of different nuclei () for matter with = 1 MeV, = fm, and = 0.2. In upper panel different colors indicate the mass fraction using a scale. In lower panel different lines indicate the contours of mass fraction in scale from 1 to 7. The total mass fraction of heavy nuclei is 56.5% (the rest is free neutrons). The majority of the mass fractions are centered around = 25 30, with a wide range of other nuclei with smaller abundances. In Figure 5 for a higher =0.4, the total mass fraction of heavy nuclei is 99.1% (free neutrons and protons have 0.8% and 0.1% respectively). The majority mass fractions are centered around = 30 35.
In Fig. 6 we show mass fractions of different nuclei for matter with = 1 MeV, = 0.4, and = fm. The total mass fraction of heavy nuclei is near one and the distribution is centered around = 35 and 50. The mass distribution of heavy nuclei in this high case is a doublepeaked Gaussian distribution, as shown in Fig. 8, where is sum of the abundances of heavy nuclei with same proton number . At a higher = 3.16 MeV shown in Fig. 7, the total mass fraction of heavy nuclei is 22.3% (free neutron, proton and alpha particles have 24.3%. 6.3% and 47.1% respectively). Here the distribution of heavy nuclei is centered around a lower Z region and has multiple peaks in as plotted in Fig. 8. It is also interesting to note that in this case there is oddeven effect in the value of for the abundances .
It is instructive to compare the composition of matter in the Virial EOS with the existing EOS tables of LattimerSwesty and H Shen et al.. The location of the neutrinosphere in a supernova is sensitive to the composition of matter and is important for the emitted neutrino spectra. Studies of collective flavor oscillations of neutrinos during their streaming outside of the neutrinosphere have already indicated a sensitivity of neutrino oscillations to the emitted neutrino spectra from the neutrinosphere Duan10 . Below we will compare some examples for the composition of matter around the neutrinosphere from the Virial EOS, the LS EOS, and SS EOS.
Fig. 9 shows the mass fractions of neutrons, protons, alpha particles, and nuclei in matter at densities from 10 to 10 fm. The matter has a temperature of 3.16 MeV and a proton fraction of 0.1 (a) or 0.3 (b), respectively. In panel (a) for a proton fraction of 0.1, free neutrons and protons dominate until the density reaches 10 fm in all three EoSs. Free nucleons dominate at low densities because of the high entropy. Above 10 fm, alpha particles appear. The SS EOS is close to our Virial results at densities below roughly 10 fm. The LS EOS significantly underestimates X and this may be due to an error in the alpha particle binding energy. Alpha particles have larger abundances and exist up to higher densities in our Virial EOS than in the other EOSs. This is partly because the attractive interactions between neutrons and alpha particles in the Virial expansion favors more alpha particles HS05 . Heavy nuclei begin to appear around 410 fm in the LS EOS, and at higher densities in the SS EOS and our Virial EOS. Moreover, the LS EOS predicts the largest abundance for heavy nuclei, while ours predicts the smallest abundance. Free neutrons have the largest abundance in our Virial EOS. This is due to the strong attractive interaction between neutrons in the Virial expansion which lowers the energy and enhances the abundance of neutrons. Note in the = 0.1 case the VirialHartree transition happens at 0.0158 fm. The right panel of Fig. 9, for a different proton fraction of 0.3, has similar characteristics. However here, alpha particles and heavy nuclei have much larger abundances than for the = 0.1 case, since a higher proton fraction favors formation of nuclei. In this = 0.3 case, the transition density from Virial gas to Hartree mean field calculations occurs at 6.3 10 fm as indicated by the dotted line in the figure.
In future work EOSIII we will carefully interpolate the free energy results of this paper in a thermodynamically consistent way in order to accurately determine derivatives of the free energy with respect to temperature or density. From these derivatives we will calculate a number of additional quantities such as the pressure or entropy.
V Summary and Outlook
In this paper we present our model for nuclear matter at subnuclear density, the Virial expansion for a nonideal gas consisting of neutrons, protons, alpha particles and thousands of heavy nuclei. We include second order virial corrections for light elements 4, nuclear partition functions for heavy nuclei, and Coulomb corrections. At very low density, the Virial expansion reduces to nuclear statistical equilibrium. We calculate the free energy, and tabulate the resulting EOS at over 73,000 grid points in the temperature range = 0.158 to 15.8 MeV, the density range = 10 to 0.1 fm, and the proton fraction range = 0.05 to 0.56. These calculations took over 1000 CPU days. The above parameter space is complementary to that of the relativistic mean field results in our previous paper SHT10a .
The treatment of Coulomb corrections in WignerSeitz approximation agrees reasonably with an analytical cluster expansion. Our results do not appear to be very sensitive to the mass table employed, or to the form of the partition function for heavy nuclei. However, results at higher densities for mean field calculations are sensitive to the interaction employed. In the future, we intend to match the Virial results of this paper to mean field calculations using a number of different interactions, see for example IUFSU .
Our EOS includes broad distributions of nuclei, that are not included in two commonly used EOS tables. These distributions of nuclei make the composition of nuclear matter in our EOS different from LS and SS EOS tables. These differences may be important for neutrino interactions.
This paper provides the second part of our results for a complete EOS table that will cover a broad range of temperatures, densities, and proton fractions for use in supernova and neutron star merger simulations. In the future we will use a thermodynamically consistent interpolation scheme to match the Virial EOS in this paper and the RMF EOS from our previous paper, or similar RMF models, and generate a complete EOS table EOSIII . Our Virial EOS is exact in the low density limit. Finally the Virial expansion also provides a suitable framework for future work that includes other modern mass tables, such as HFB14 HFB14 or KTUY05 KTUY05 .
Vi Acknowledgement
We thank Lorenz Hedepohl, Thomas Janka, Jim Lattimer, Andreas Marek, Evan O’Connor, and Christian Ott for helpful discussions. This work was supported in part by DOE grant DEFG0287ER40365. This material is based upon work supported by the National Science Foundation under Grants No. ACI0338618l, OCI0451237, OCI0535258, OCI0504075 and CNS0521433. This research was supported in part by Shared University Research grants from IBM, Inc. to Indiana University and by the Indiana METACyt Initiative. The Indiana METACyt Initiative of Indiana University is supported in part by Lilly Endowment.
References
 (1) M. Liebendörfer, M. Rampp, H.Th. Janka, and A. Mezzacappa, Astrophys. J. 620, 840 (2005).
 (2) R. Walder, A. Burrows, C. D. Ott, E. Livne, I. Lichtenstadt, and M. Jarrah, Astrophys. J. 626, 317 (2005).
 (3) A. Burrows, S. Reddy, and T. A. Thompson, Nucl. Phys. A 777, 356 (2006).
 (4) G. Shen, C. J. Horowitz, and S. Teige, Phys. Rev. C 82, 015806 (2010).
 (5) P. Möller, J. R. Nix, and W. J. Swiatecki, Atom. Data Nucl. Data Tables, 59, 185 (1995).
 (6) J. M. Lattimer and F. D. Swesty, Nucl. Phys. A 535, 331 (1991).
 (7) H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi, Nucl. Phys. A 637, 435 (1998).
 (8) H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi, Prog. Theo. Phys. 100, 1013 (1998).
 (9) M. Costantinit, A. Ianni, and F. Visanni, Phys. Rev. D 70, 043006 (2004); C. Lunardini and A.Y. Smirnov, Astropart. Phys. 21, 703 (2004).
 (10) 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).
 (11) B.S. Meyer, Ann. Rev. Astron. Astrophys. 32, 153 (1994).
 (12) S. I. Blinnikov, I. V. Panov, M. A. Rudzsky, and K. Sumiyoshi, arXiv: 0904.3849.
 (13) Matthias Hempel and Jürgen SchaffnerBielich, arXiv:0911.4073v1.
 (14) See for example, J. Carlson, S.Y. Chang, V. R. Pandharipande, and K. E. Schmidt, Phys. Rev. Lett. 91, 050401 (2003).
 (15) T.L. Ho and E. J. Mueller, Phys. Rev. Lett. 92, 160404 (2004).
 (16) C. J. Horowitz and A. Schwenk, Phys. Lett. B638, 153 (2006).
 (17) C. J. Horowitz and A. Schwenk, Nucl. Phys. A 776, 55 (2006).
 (18) C. J. Horowitz and A. Schwenk, in preparation.
 (19) W. A. Fowler, C. A. Engelbrecht, and S. E. Woosley, ApJ 226, 984 (1978).
 (20) T. Rauscher, F.K. Thielemann, and K.L. Kratz, Phys. Rev. C 56, 1613 (1997).
 (21) A. Gilbert and A. G. W. Cameron, Canadian Journal of Physics. 43, 1446 (1965).
 (22) T. J. Mazurek, J. M. Lattimer, and G. E. Brown, Astrophysical Journal 229, 713 (1979).
 (23) C. A. Engelbrecht and J. R. Engelbrecht, Ann. Phys. 207, 1 (1991).
 (24) H. A. Bethe, Phys. Rev. 50, 332 (1936).
 (25) T. Rauscher and F.K. Thielemann, Atom. Data Nucl. Data Tab. 75, 1 (2000).
 (26) D. J. Dean, S. E. Koonin, K.H. Langanke, P. B. Radha, and Y. Alhassid, Phys. Rev. Lett. 74, 2909 (1995).
 (27) V. Paar and R. Pezer, Phys. Rev. C 55, R1637 (1997).
 (28) N. V. Brilliantov, Contrib. Plasma Phys. 38, N4 (1998).
 (29) G. Baym, H. A. Bethe, and C. J. Pethick, Nucl. Phys. A 175, 225 (1971).
 (30) D. K. Nzdyozhin and A. V. Yudin, Astro. Lett. 30, 634 (2004).
 (31) S. Goriely, M. Samyn, and J. M. Pearson, Phys. Rev. C 75, 064312 (2007).
 (32) See for a recent review, H. Duan, G. M. Fuller, and Y.Z. Qian, arXiv:1001.2799.
 (33) G. Shen, C. J. Horowitz, and S. Teige, to be published.
 (34) F. J. Fattoyev, C. J. Horowitz, J. Piekarewicz, and G. Shen, arXiv:1008.3030.
 (35) H. Koura, T. Tachibana, M. Uno, and M. Yamada, Prog. Theo. Phys. 113, 305 (2005).