Equation of State of Dense Matter from a density dependent relativistic mean field model

# Equation of State of Dense Matter from a density dependent relativistic mean field model

## Abstract

We calculate the equation of state (EOS) of dense matter, using a relativistic mean field (RMF) model with a density dependent coupling that is a slightly modified form of the original NL3 interaction. For nonuniform nuclear matter we approximate the unit lattice as a spherical Wigner-Seitz cell, wherein the meson mean fields and nucleon Dirac wave functions are solved fully self-consistently. We also calculate uniform nuclear matter for a wide range of temperatures, densities, and proton fractions, and match them to non-uniform matter as the density decreases. The calculations took over 6,000 CPU days in Indiana University’s supercomputer clusters. We tabulate the resulting EOS at over 107,000 grid points in the proton fraction range = 0 to 0.56. For the temperature range = 0.16 to 15.8 MeV we cover the density range = 10 to 1.6 fm; and for the higher temperature range = 15.8 to 80 MeV we cover the larger density range = 10 to 1.6 fm. In the future we plan to study low density, low temperature (T15.8 MeV), nuclear matter using a Virial expansion, and we will match the low density and high density results to generate a complete EOS table for use in astrophysical simulations of supernova and neutron star mergers.

###### pacs:
21.65.Mn,26.50.+x,26.60.Kp,21.60.Jz,97.60.Bw

## I Introduction

The equation of state (EOS) for hot, dense matter in massive stars relates energy and pressure to temperature, density, and composition. It has been a long-standing problem to understand the EOS at both subnuclear and supranuclear density, to which great efforts have been devoted, from laboratory heavy ion collision experiments heavyions (), computer simulations of supernova simulation1 (); simulation2 (), and theoretical many-body calculations theory (). The EOS of hot dense matter in supernovae (SN) and neutron star (NS) mergers encompass multi-scale physics. Temperature can vary from 0 to as high as 100 MeV, exciting nuclei, nucleon and possibly pion and other degrees of freedom. The density can vary from to 10 g cm, where matter can be in gas, liquid or solid phases. The proton fraction can vary from 0 to 0.6, from extremely neutron rich matter to proton rich matter. These very large parameter ranges make construction of a full EOS table difficult. It is necessary to employ different approximations for different parameter ranges. As a result, there exist only two realistic EOS tables that are in widespread use for astrophysical simulations, the Lattimer-Swesty (L-S) equation of state LS (), that uses a compressible liquid drop model with a Skyrme force, and the H. Shen, Toki, Oyamatsu and Sumiyoshi (S-S) equation of state Shen98a (); Shen98 (), that uses the Thomas Fermi and variational approximations with a relativistic mean field (RMF) model. We plan to generate a complete equation of state, employing relativistic mean field calculations for matter at intermediate and high density as described in this paper. In the future we plan to use the Virial expansion of a nonideal gas to describe matter at low density. These two parts will be matched together and we will generate a thermodynamically consistent EOS over the full range of parameters. Finally we will generate additional EOSs from RMF models with different high density symmetry energies. This will allow one to correlate features of astrophysical simulations with properties of the symmetry energy assumed for the EOS.

There are still large uncertainties in the EOS at supranuclear densities. The density dependence of the symmetry energy is poorly known and strongly influences the stiffness of the EOS. It can be constrained from measurements of NS radii and masses Lattimer04 (), precision determination of the neutron rms radius in Pb prex (), and also heavy ion collision experiments heavyions (). A stiff EOS (high pressure) at high density gives larger NS radii, while a stiff EOS at normal and low density favors a larger neutron radius in Pb HP01 (). The elliptic and transverse flow observables in heavy ion collisions are sensitive to the isospin dependence of mean fields and to the EOS at densities up to a few times nuclear saturation density. Many nuclear many-body models fall into two categories, the non-relativistic Skyrme models (See for example, Ref. Talmi () for a review) and relativistic mean field models RMFreviews (). The parameters in these models are usually fitted to nuclear properties at normal nuclear densities, afterwards they are extrapolated to study supranuclear matter. The L-S EOS uses a Skyrme model featuring a relatively soft EOS and the S-S EOS uses the RMF interaction TM1 that features a stiffer EOS. Since the symmetry energy is not well constrained, it is important to explore the effects of different symmetry energies on the EOS and SN simulations.

In this paper, we use a RMF model for non-uniform matter at intermediate density and uniform matter at high density. Low density pure neutron matter is analogous to a unitary gas HS05b (), where the neutron-neutron scattering length is much larger than both the effective range and the average inter-particle spacing. To better describe neutron-rich matter at low density, we use a density dependent scalar meson-nucleon coupling. At high density, the model reduces to the normal RMF parameter set NL3. The unit lattice of non-uniform nuclear matter is conveniently approximated by a spherical Wigner-Seitz (W-S) cell. The meson mean fields and nucleon Dirac wave functions inside the Wigner-Seitz cell are solved fully self-consistently. This is unlike the S-S EOS that used Thomas-Fermi and variational approximations and the L-S EOS that used a simple liquid drop model. The size of the W-S cell is found by minimization of the free energy per nucleon. The W-S approximation provides a framework to incorporate the best known microscopic nuclear physics Negele73 (). The nuclear shell structure effects are included automatically and it is already possible for some effects of complex nuclear pasta states to be included in spherical calculations in the form of shell states HS08 (). Full three-dimensional W-S calculations in principle could incorporate various pasta shapes Gogelein07 (); Newton09 (), which would make the transition to uniform matter more smooth. However this will demand much larger computational resources. In this work we use the spherical W-S approximation.

Our relativistic mean field calculations can accurately describe the radial shape of large neutron rich nuclei including the expected neutron rich skin. In contrast the original L-S EOS is based on a very simple liquid drop model of nuclear structure that may incorrectly describe the neutron skin. Alternatively the S-S EOS is based on a Thomas Fermi approximation that neglects shell effects. These are included in our Hartree calculations. Furthermore, the variational forms for the densities assumed by S-S may be a poor approximation for large proton numbers where the Coulomb repulsion is large. Instead our exact solutions of the radial mean field equations allow richer density distributions including shell states with central depressions HS08 (). These differences in densities may be important for neutrino interactions in Supernovae. Finally our calculations correctly reproduce the Unitary gas limit for a low density neutron gas, see below, while both the L-S and S-S EOS reduce incorrectly to the energy of free neutrons.

One can demand that any EOS be consistent with, possibly model dependent interpretations of, observations of neutron stars. For example, Klahn et. al. propose a series of tests that an EOS should satisfy to be consistent with observations klahn (). They demand that any reliable nuclear EOS be able to reproduce the recently reported high pulsar mass of 0.2 for PSR J0751+1807 nice (). However, this observation may have been retracted notnice (). Furthermore, Klahn et al. require the EOS to reproduce a large binding energy for Pulsar B in J0737-3039. However, this conclusion could be sensitive to assumptions about the system such as the amount of mass loss. Klahn et. al. go on to demand that the EOS not allow direct URCA cooling of neutron stars of mass 1 to 1.5 . We consider a more conservative approach. While many stars cool slowly, observations do suggest that at least some stars have enhanced cooling. Unfortunately observations do not directly constrain the mass that may separate enhanced from normal cooling. Indeed, there is little direct observational evidence that more massive stars cool more quickly, although this is a theoretical prejudice.

One can also use laboratory data to constrain the EOS. The neutron skin thickness of a heavy nucleus constrains the density dependence of the symmetry energy. Furthermore, there are many measurements of the skin thickness with a variety of strongly interacting probes. However, there may be important model dependence from strong interaction uncertainties. For example measurements of spin dipole strength have been used to extract neutron skin thicknesses in Sn isotopes spindipole (). For these measurements, the spin dipole strength was assumed to be proportional to the measured cross section, and the proportionality constant was arbitrarily fixed in order to reproduce the skin thickness of Sn as predicted by an old Hartree-Fock calculation oldHF (). Presumably, if a different skin thickness in Sn is fit, one would also get different answers for the skin thickness in other isotopes.

This situation may soon change. The Lead Radius Experiment (PREX) at Jefferson Laboratory is using parity violating electron scattering to measure the neutron skin thickness in Pb prex (). Parity violation is a sensitive probe of neutrons because the weak charge of a neutron is much larger than that of a proton. Furthermore, this electro-weak reaction may have much smaller strong interaction uncertainties. Data taking for PREX should be completed by June 2010.

Instead of trying to determine, ahead of time, the best EOS to satisfy existing observational constraints, we adopt what we hope will be a more robust approach. We are calculating a number of EOSs based on different effective interactions. In this paper we present first results for the NL3 interaction with a symmetry energy that is large at high densities. In later work we will present EOSs with softer high density symmetry energies. These different EOS will allow one to correlate features of astrophysical simulations with properties of the EOS. Then one can draw conclusions based on combined information from laboratory experiments and astronomical observations.

In this paper we focus on nucleon degrees of freedom. Hyperons could play a role at high densities, see for example shen2 (). However, the contribution of hyperons could depend on uncertain hyperon interactions. In addition, there could be pion or kaon condensates or a variety of quark matter phases. See for example the review by Page and Reddy theory (). Chiral symmetry restoration and the softening of pionic or kaonic modes could be important. Finally, thermal pions and pion interactions should be very important at high temperatures. All of these effects may increase the uncertainties in the EOS.

We tabulate the equation of state for intermediate and high density nuclear matter over the range of temperatures , densities , and proton fractions given in Table 1 and described in Sec. IV. We calculate the free energy of nonuniform matter for 17021 points, and the free energy of uniform matter for 90478 points in , , and space. This took 6000 CPU days on Indiana University’s supercomputer clusters.

The paper is organized as follows: in Section II the density dependent RMF model is explained in detail. In Section III we describe the RMF parameters that we use including a density dependent coupling. We describe the computational methodology for our large computer runs in Section IV. Section V shows results for RMF calculations, including the free energy and the nucleon density distributions in the non-uniform W-S cells. Finally, Section VI presents a summary of our results and gives an outlook for future work.

## Ii Formalism

We now describe the mean field formalism that we use for non-uniform matter in Section II.1 and for uniform matter in Section II.2.

### ii.1 Non-uniform nuclear matter in Wigner-Seitz approximation

The formalism for relativistic mean field theory has been reviewed in previous papers, see eg RMFreviews (). To better describe neutron rich matter at low density we introduce a density dependent coupling between the scalar meson and the nucleon as described in Section III. We note that many previous studies of density dependent RMF models mainly focused on better descriptions of nuclear matter at supranuclear density (see for example, Typel99 (); Ban04 ()). In this section we focus on low density neutron rich matter.

The basic ansatz of the RMF theory is a Lagrangian density where nucleons interact via the exchange of sigma- (), omega- (), and rho- () mesons, and also photons ().

 L = ¯¯¯¯ψ[iγμ∂μ−m−Γσσ−gωγμωμ (1) − gργμ→τ⋅→\boldmathρμ−eγμ1+τ32Aμ]ψ + 12∂μσ∂μσ−12m2σσ2−13g2σ3 − 14g3σ4 − 14ωμνωμν+12m2ωωμωμ+14c3(ωμωμ)2 − 14→\boldmathρμν⋅→\boldmathρμν+12m2ρ→\boldmathρμ⋅→\boldmathρμ− 14AμνAμν

We note that ( and is nucleonic current) is the density dependent coupling between the sigma meson and the nucleon. Here the field tensors of the vector mesons and the electromagnetic field take the following forms:

 ωμν = ∂μων−∂νωμ, Aμν = ∂μAν−∂νAμ, →\boldmathρμν = ∂μ→\boldmathρν−∂ν→\boldmathρμ−gρ→\boldmathρμ×→\boldmathρν. (2)

In charge neutral nuclear matter composed of neutrons, , protons, , and electrons, , there are equal numbers of electrons and protons. Electrons can be treated as a uniform Fermi gas at high densities 4. They contribute to the Coulomb energy of the matter and serve as one source of the Coulomb potential.

The variational principle leads to the following equations of motion

 [\slα⋅\slp + V(r) + \slβ(m + S(r))]ψi = εiψi (3)

for the nucleon spinors, with vector and scalar potentials

 V(r)= β{gωto0.0pt/ωμ+gρ→τ⋅→to0.0pt/ρμ+e(1+τ3)2 to0.0pt/Aμ+ΣR},S(r)= Γσσ, (4)

where

 ΣR=γμjμn∂Γσ∂nρsσ, (5)

is the rearrangement term due to the density dependent coupling between the sigma meson and the nucleon, and is the scalar density of nucleons to be defined below.

The equations of motion for the mesons and photons are,

 (m2σ − ∇2) σ= −Γσρs− g2σ2 − g3σ3,(m2ω − ∇2)ωμ= gωjμ−c3ωμ(ωνων),(m2ρ − ∇2)→ρμ= gρ→jμ,− ∇2Aμ= e(jμp−jμe), (6)

where the electrons are included as a source of Coulomb potential. The nucleon spinors provide the relevant source terms:

 ρs= ∑i¯¯¯¯ψiψini,jμ= ∑i¯¯¯¯ψiγμψini,→jμ= ∑i¯¯¯¯ψiγμ→τψini,jμp= ∑i¯¯¯¯ψiγμ1+τ32ψini.

At finite temperature, Fermi-Dirac statistics imply the occupations, , of protons and neutrons are:

 ni = 1eβ(εi−μ)+1, (8)

where is the chemical potential for neutron (proton). In the calculation, we include all levels with , where is the degeneracy of the level.

Since the systems under consideration have temperatures of, at most, tens of MeV, we neglect the contribution of negative energy states, i.e., the so-called no sea approximation. In a spherical nucleus, there are no currents in the nucleus and the spatial vector components of , and vanish. One is left with the timelike components, , and . Charge conservation guarantees that only the 3-component of the isovector field survives. The above non-linear equations are solved by iteration within the context of the mean field approximation whereby the meson field operators are replaced by their expectation values.

The spherical Wigner-Seitz approximation is used to describe non-uniform matter. One W-S cell has one nucleus. In this approximation it is important to include lattice Coulomb corrections between neighboring W-S cells. The detailed treatments we use have been discussed in a previous paper HS08 () and we will not repeat them here.

We solve for the meson mean fields and the nucleon Dirac wave functions self-consistently inside a W-S cell of radius , for a given baryon density , proton fraction and temperature . In our RMF model nucleons (proton and neutron) are the only baryons. The nucleon number inside the W-S cell is and the proton number is . The internal energy of a W-S cell, including the approximate lattice Coulomb energy correction, is,

 Eb = Enucleon+Eσ+Eρ+Eω+ECoul−mA, (9) = ∑iϵini− ∫d3rj0(r)∂Γσ∂j0ρs(r)σ(r)− 12∫d3r{Γσσρs(r) + 13g2σ3 + 12g3σ4}, (11) − 12∫d3rgρρ0,3j0,3(r) − 12∫d3r{gωω0 j0(r)− 12c3ω40}− e2∫(ρp+ρe)A0(r)d3r + dw−mA,

where is the approximate Coulomb correction for a bcc latticeOyamatsu (), and is the volume of W-S cell.

The nucleon contribution to the entropy is given by the usual formula,

 Sb=−kB∑igi[niln(ni) +(1−ni)ln(1−ni)], (12)

where is given in Eq. (8). With Eqs. (9) and (12), it is easy to obtain the nucleon contribution to the free energy per nucleon ,

 F=Fb/A=(Eb − TSb)/A. (13)

### ii.2 Uniform nuclear matter

To make the paper self-contained, we give the formulas for uniform matter in RMF model. As we show below, at high temperatures or high densities the matter is uniform. We include anti-nucleon terms which make a small contribution at very high temperatures.

The energy density of uniform nuclear matter is,

 ε = ∑i=N,Pεikin + 12[m2σσ2 + m2ωω20 + m2ρρ20,3] (15) + 13g2σ3 + 14g3σ4 + 34c3ω40,

where

 εikin = 2(2π)3∫d3kE∗(k)[nk(T)+¯nk(T)], (16)

with effective mass , and . The occupation probabilities for particles and antiparticles are,

 nk(T) = 1exp(E∗(k)+gωω0+gρτ3ρ0,3+∂Γσ∂nρsσ−μ)/T+1, (17) ¯nk(T) = 1exp(E∗(k)−gωω0−gρτ3ρ0,3−∂Γσ∂nρsσ+μ)/T+1. (18)

The pressure of uniform nuclear matter is,

 P = ∑i=N,PPikin − 12m2σσ2−13g2σ3−14g3σ4 (20) +∂Γσ∂nρsσn+ 12m2ωω20 + 14c3ω40 + 12m2ρρ20,3,

where

 Pikin = 23(2π)3∫d3kk2√k2+m∗2[nk(T)+¯nk(T)]. (22)

The entropy density of uniform nuclear matter is,

 s=−2kB(2π)3∫d3k[nk(T)lnnk(T) + (1−nk(T))ln(1−nk(T)) + ¯nk(T)ln¯nk(T) + (1−¯nk(T))ln(1−¯nk(T))]. (23)

Using Eqs. (15) and (23) one can obtain the free energy density per nucleon for uniform matter,

 F = (ε−Ts)/nB. (24)

## Iii Parameter set with density dependent coupling

In this work we use the NL3 effective interaction NL3 () that has been successful in reproducing ground state properties of stable nuclei and the saturation properties of symmetric nuclear matter. The values of parameters in the NL3 effective interaction are listed in Table 2.

As is well known, the mean field approach for pure neutron matter is problematic at low densities because long range correlations are important. Neutron matter at low density is very close to a unitary gas HS05b (), since the scattering length is much larger than the inter-particle spacing, which is also larger than the effective range of nuclear interaction. To describe neutron matter phenomenologically in the RMF framework, without losing its success for the properties of nuclear matter, we introduce a density-dependent scalar meson-nucleon coupling,

 Γσ = ⎧⎪⎨⎪⎩Γ0σ,n>n0Γ0σ1+α[(n+n02n)16+α],n≤n0. (25)

The two free parameters and are determined by matching the energy of neutron matter to that of a unitary gas at zero temperature unitarygas (),

 EU = ξ⋅35k2F2m≃0.44⋅35k2F2m, (26)

where is the neutron Fermi momentum. The best fitted values are fm and .

In Fig. 1, the energy of pure neutron matter at = 0 is shown for the original NL3 set, the modified NL3 set with a density dependent -N coupling as in Eq. (25), and the unitary gas calculated by Eq. (26). The unitary gas gives lower energy than the original NL3 result by about 0.2 MeV per particle, due to the strong wave attractive interactions. This energy difference is very relevant for matching a Virial expansion to the mean field calculations, since the Virial expansion includes the long range two-body neutron-neutron attractive interaction SHT10b (); HS05 () while the normal mean field calculation does not. In the density range shown in the figure, NL3 also gives a lower energy than the TM1 or FSUGold Fsugold () RMF parameter sets. However, the density dependent NL3 can fit the unitary gas result by tuning the coupling strength in the attractive scalar meson channel. Therefore, the density dependent NL3 set describes successfully the properties of both neutron rich matter and nuclear matter. In the following when we refer to NL3 set, we mean the density dependent NL3 set unless otherwise specified.

## Iv Computational Methodology

In this section we describe our strategy for evaluating the equation of state. We calculate the equation of state for the following partitioning of , , and parameter space, see Table 1:

• We use a step of 0.2 in for from -0.8 to 0, a step of 0.1 for from 0 to 1.1 and a step of 0.05 for above 1.1. We have a total of 32 points for from 0.16 to 80 MeV.

• For temperatures below 15.8 MeV (where matter can be nonuniform) we use a step of 0.1 in for from -4.0 to 0.2. We have a total of 43 points for from to 1.6 fm.

• For temperatures above 15.8 MeV (where matter is uniform) we use a step of 0.1 in for from -8.0 to 0.2. We have a total of 83 points for from to 1.6 fm.

• We use a step of 0.01 in proton fraction for from 0.05 to 0.56. We aslo include . This gives a total of 53 points for from 0.0 to 0.56.

This partitioning gives a total of 40,248 points in the nonuniform Hartree region. However for matter at higher temperatures, but still MeV, and/or lower proton fractions, Hartree results give higher free energies than corresponding results for uniform matter. By roughly estimating the phase boundary, and keeping enough points to cross the transition density, see Section V , we calculate the free energy for a reduced number of points in the nonuniform Hartree region, see Sec. II.1, which includes a total of 17,021 points. We also calculate free energies for uniform nuclear matter at a total of 90,478 points, see Sec. II.2.

The most time is spent evaluating (temperature , proton fraction , density ) points in the non-uniform Hartree mean field region. For each point we need to minimize the free energy of the W-S cell with respect to the cell radius which typically requires evaluation at 40 to 100 cell radii. This minimization can be complicated by the existance of local minima. For each cell size, we need to solve the mean fields self consistently. We have already developed the code for this minimization (see Ref.  HS08 ()) which is slightly modified in this work to accomodate the density dependent coupling in RMF.

The mean fields provide potentials for the individual nucleons in the W-S cell which obey the Dirac equation. The Dirac equation is solved by a fourth order Runge-Kutta method with shooting techniques. For nuclear matter at finite temperature, there could be hundreds of nucleons that populate thousands of levels according to Fermi-Dirac statistics. For each level, we need to solve the Dirac equation. The potentials for the nucleons in the Dirac Eq. (4) are various meson mean fields which obey the extended Klein-Gordon (K-G) equation. The source terms for the K-G equations are provided by various nucleon density terms in Eq. (II.1). Given the nucleon density terms, the K-G equations are solved by a Green’s function method, which updates the meson mean fields. The updated mean fields can now be used to solve the nuclear levels and nucleon densities again. This process is repeated until full self-consistency is achieved in both the mean fields and the nuclear levels.

Computationally, the problem to be solved is embarrassingly parallel because each point of density, temperature and proton fraction is independent of the others. A total number of 17,000 independent tasks must be run, where each task calculates the required quantities at a single point in the phase space. Unfortunately, the run time on an individual task varies from a few minutes to more than 24 hours, depending on the number of iterated cell radii, and the number of nucleon energy levels included.

Each point in the phase 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 phase space, and a single character (A=available, R=running, r=Re-running, C=complete, T=time-limited and F=failed) that gives the status of calculations for that job index. An 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 phase space and it evaluates the required quantities.

Initially, rank zero of the MPI job

• locks the job listing file runlist

• 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 time-limits 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 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 phase space. A wall clock limit (48 hours) larger than the average run time is used. Each rank of the MPI job can run a series of points via 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. Using this scheme we have achieved 85% efficiency in CPU usage. Specifically, 85% of all jobs ended in the ”C” state rather than the ”T” or ”R” state. After the runlist has been searched once, the remaining jobs have ”R” or ”T” state. Then these remaining jobs are resubmitted via the MPI wrapper code, requesting longer time limit (typically 7 days) but fewer CPU cores. This procedure allows us to calculate 99 of the points in the runlist file.

## V Results

In this section we discuss our results for various regions of parameter space. First, the uniform matter EOS at zero temperature is presented. Second, we discuss the free energy per nucleon for mean field calculations of nonuniform matter. Finally, we show the density distributions of neutrons and protons inside W-S cells.

### v.1 Uniform matter at zero temperature

Figure 2 shows the equation of state of uniform matter at zero temperature with different proton fractions = 0, 0.1, 0.2, 0.3, 0.4, and 0.5. The solid curves are for the NL3 parameter set, which is used in our RMF calculations. The dashed curves show results for the TM1 interaction, which is used in the equation of state obtained by H. Shen et al Shen98 (). The two sets agree to a great extent for densities below 0.20.25 fm, depending the value of . Above these densities, NL3 gives a much stiffer equation of state for uniform matter. This serves as one motivation for our choice of the NL3 parameters: to explore the equation of state with a stiffer symmetry energy.

### v.2 Free energy and phase boundaries

In Fig. 3, the free energy per nucleon is shown as a function of density at = 1, 3.16, 6.31 and 10 MeV. At intermediate densities, is calculated from Eq. (13) for W-S cells using Hartree mean field calculations. At high densities, is calculated from Eq. (24) for uniform matter. The transition (as the density grows) is found at the density where uniform matter gives a lower free energy. In each panel, the (red) solid curves give the transition densities to uniform matter. The transition densities increase as proton fraction grows. Non-uniform matter can exist until higher densities in more symmetric nuclear matter. At density around 0.16 fm, there is always a minimum in the free energy per nucleon, as long as the proton fraction is not too small and temperature not too high. This is the manifestation of saturation density in nuclear matter.

### v.3 Density distributions inside Wigner-Seitz cells

The Hartree mean field calculation provides detailed wavefunctions for nucleons in the non-uniform phase. In our W-S approximation at intermediate densities, we find a “spherical pasta” phase where the proton density distribution forms a shell state with a reduced density in the center. This reduces the large coulomb repulsion between protons and was discussed in HS08 (). In this section we discuss density distributions inside the W-S cells, both for normal nuclei and for these shell states.

In Fig. 4, neutron and proton distributions, inside the W-S cell, are shown for four different baryon densities, with = 1 MeV and = 0.1. At very low density = 0.002 fm, the Hartree calculation has a minimum for protons and nucleons. Most of the neutrons are located within 10 fm of the cell center, although the W-S cell radius is around 31 fm. A small fraction of the neutrons extend to the edge of cell since this is an extremely neutron rich system. As the density rises to 0.02 fm, the W-S cell has , , and the neutron density at large becomes much greater. The W-S cell radius drops to 17.5 fm because the lattice becomes more closely packed as the density increases. At a density of 0.05 fm, the W-S cell has and and forms a shell state with both inside and outside surfaces. This has been discussed in our early paper HS08 (). As a result, the W-S cell radius become larger. At the higher density of 0.063 fm, the system becomes uniform.

In Fig. 5, the distribution of neutrons and protons are shown for = 0.020 fm and proton fraction =0.3. At low temperature = 1 MeV, where , , the density distributions are similar to those for normal isolated nuclei. As the temperature rises to 3.16 and 6.31 MeV, the size of W-S cell remains nearly fixed, but the neutron density increases with temperature at large radius. This is due to excitation of states with high angular momentum and/or large main quantum number as the temperature rises. When the temperature rises to 10 MeV, the proton density also rises at large , accompanied by an increase of the W-S cell size, with , . At even higher temperature, the nucleus melts and uniform matter appears.

Similar to Fig. 5 but at higher = 0.050 fm and proton fraction = 0.45, the density distribution of neutrons and protons are shown for four different temperatures in Fig. 6. Here a shell state exists up to high temperatures. At low temperatures = 1, 3.16 MeV, , and the shell state has inside and outside voids. As the temperature rises, nucleons populate both the inside and outside voids, due to thermal excitations. Finally, the size of the shell state shrinks at high temperature so that , at =10 MeV.

## Vi Summary and Outlook

In this paper we present large scale relativistic mean field calculations for nuclear matter at intermediate and high densities. We use a density dependent modification of the NL3 interaction in a spherical Wigner-Seitz approximation. Nuclear shell effects are included. We calculate the free energy, and tabulate the resulting equation of state at over 107,000 grid points in the proton fraction range = 0 to 0.56. For low temperatures = 0.16 to 15.8 MeV we calculate for the density range = 10 to 1.6 fm. For high temperatures = 15.8 to 80 MeV, where the matter is uniform, we calculate for the larger density range = 10 to 1.6 fm. These calculations took over 6000 CPU days.

We solve for the nucleon Dirac wave functions and meson mean fields self-consistently. This allows us to study how the distribution of neutrons and protons inside a Wigner Seitz cell evolve with density and temperature. We find a large variety of possible sizes and shapes for these distributions.

This paper provides part of our results for an equation of state which will cover a broad range of temperatures, densities, and proton fractions. In the future, we plan to study low density nuclear matter using a Virial expansion for a non-ideal gas consisting of nucleons and thousands of species of nuclei. Then, we will generate a complete thermodynamically consistent equation of state by matching the low density and higher density results. This equation of state avoids the Thomas Fermi and variational approximations of the H. Shen et al. equation of state and is exact in the low density limit. It can be used in supernova and neutron star merger simulations. Finally, in future work we will generate equations of state using other modern relativistic mean field interactions such as FSUGold Fsugold ().

## Vii 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 DE-FG02-87ER40365. This material is based upon work supported by the National Science Foundation under Grants No. ACI-0338618l, OCI-0451237, OCI-0535258, OCI-0504075 and CNS-0521433. 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.

### Footnotes

1. e-mail: gshen@indiana.edu
2. e-mail: horowit@indiana.edu
3. e-mail: steige@indiana.edu
4. It needs electron density g/cm, which is easily surmounted in the regime of mean field results.

### References

1. See for example C. Hartnack et al., Phys. Rev. Lett. 96, 012302 (2006). P. Danielewicz, R. Lacey, and W. G. Lynch, Science 22, 1592 (2002).
2. M. Liebendörfer, M. Rampp, H.-Th. Janka, A. Mezzacappa, Astrophys. J. 620, 840 (2005).
3. R. Walder, A. Burrows, C.D. Ott, E. Livne, M. Jarrah, Astrophys. J. 626, 317 (2005).
4. D. Page and S. Reddy, Ann. Rev. Nucl. Part. Sci. 56, 327 (2006).
5. J. M. Lattimer and F. D. Swesty, Nucl. Phys. A 535, 331 (1991).
6. H. Shen, H. Toki, K. Oyamatsu and K. Sumiyoshi, Nucl. Phys. A 637, 435 (1998).
7. H. Shen, H. Toki, K. Oyamatsu and K. Sumiyoshi, Prog. Theo. Phys. 100, 1013 (1998).
8. J. M. Lattimer and M. Prakash, Science 304, 536 (2004).
9. C. J. Horowitz, S. J. Pollock, P. A. Souder, and R. Michaels, Phys. Rev. C 63, 025501 (2001).
10. C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001).
11. I. Talmi, Adv. Nucl. Phys. 27, 1 (2003).
12. B. D. Serot and J. D. Walecka, Adv. Nucl. Phys. 16, 1 (1986); P.-G. Reinhard, Rep. Prog. Phys. 52, 439 (1989); P. Ring, Prog. Part. Nucl. Phys. 37, 193 (1996).
13. C. J. Horowitz and A. Schwenk, Phys. Lett. B638, 153 (2006).
14. J. W. Negele and D. Vautherin, Nucl. Phys. A 207, 298 (1973).
15. C. J. Horowitz and G. Shen, Phys. Rev. C 78, 015801 (2008).
16. P. Gögelein, H. Müther, Phys. Rev. C 76, 024312 (2007).
17. W. G. Newton and J. R. Stone, Phys. Rev. C 79, 055801 (2009).
18. T. Klahn et al., Phys. Rev. C 74, 035802 (2006).
19. D. J. Nice et al., Astrophys. J. 634, 1242 (2005).
20. D. J. Nice, Am. Astron. Soc., AAS Meeting #211, #102.06; Bulletin of the AAS 39, 918 (2007).
21. A. Krasznahorkay et al., Phys. Rev. Lett. 82, 3216 (1999).
22. I. Angeli et al., J. Phys. G 6, 303 (1980).
23. H. Shen, Phys. Rev. C 65, 035802 (2002).
24. S. Typel and H. H. Wolter, Nucl. Phys A 656, 331 (1999).
25. S. F. Ban, J. Li, S. Q. Zhang, H. Y. Jia, J. P. Sang, and J. Meng, Phys. Rev. C 69, 045805 (2004).
26. K. Oyamatsu, Nucl. Phys. A 561, 431 (1993).
27. G. A. Lalazissis, J. König, and P. Ring, Phys. Rev. C55, 540 (1997).
28. G. Shen, C. J. Horowitz and S. Teige, to be submitted.
29. See for example, J. Carlson, S.-Y. Chang, V. R. Pandharipande, and K. E. Schmidt, Phys. Rev. Lett. 91, 050401 (2003).
30. C. J. Horowitz and A. Schwenk, Nucl. Phys. A 776, 55 (2006).
31. B. G. Todd-Rutel and J. Piekarewicz, Phys. Rev. Lett. 95, 122501 (2005).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters