Hot and Dense Homogeneous Nucleonic Matter Constrained by Observations, Experiment, and Theory
Xingfu Du Andrew W. Steiner Jeremy W. Holt






^1Department of Physics and Astronomy, University of
  Tennessee, Knoxville, TN 37996, USA
^2Physics Division, Oak Ridge National Laboratory, Oak
  Ridge, TN 37831, USA
^3Cyclotron Institute and Department of Physics and
  Astronomy, Texas A&M University, College Station, Texas 77843, USA


  We construct a new class of phenomenological equations of state for
  homogeneous matter for use in simulations of hot and dense matter in
  local thermodynamic equilibrium. We construct a functional form
  which respects experimental, observational and theoretical
  constraints on the nature of matter in various density and
  temperature regimes. Our equation of state matches (i) the virial
  coefficients expected from nucleon-nucleon scattering phase shifts,
  (ii) experimental measurements of nuclear masses and charge radii,
  (ii) observations of neutron star radii, (iii) theory results on the
  equation of state of neutron matter near the saturation density, and
  (iv) theory results on the evolution of the EOS at finite
  temperatures near the saturation density. Our analytical model
  allows one to compute the variation in the thermodynamic quantities
  based on the uncertainties in the nature of the
  nucleon-nucleon interaction. Finally, we perform a correction to
  ensure the equation of state is causal at all densities,
  temperatures, and electron fractions. 


97.60.Jd, 95.30.Cq, 26.60.-c

Hot and Dense Homogeneous Nucleonic Matter Constrained by Observations,
  Experiment, and Theory
    Jeremy W. Holt^3
    December 30, 2023
================================================================================================





§ INTRODUCTION


The equation of state (EOS) of dense nucleonic matter is a central
microphysical input required for numerical simulations of
core-collapse supernovae <cit.> and neutron star
mergers <cit.>. These simulations probe baryon densities
(n_B) up to a few nucleons per fm^3, temperatures (T) up to 100
MeV, and a wide range of electron fractions (Y_e). This
three-dimensional space is often described through an EOS table which
relates the free energy or pressure of the system to the thermodynamic
variables describing the ambient environment. These tables, in turn,
are built from modeling the interactions between nucleons, nuclei, and
potentially exotic particles that may appear at high densities.

Over this large three-dimensional (n_B,Y_e,T) space, there are
several different physical regimes each constrained by different
observables and theoretical approaches. The first regime, zero
temperature nuclear matter at nuclear saturation density, is closely
connected to nuclear masses, charge radii, giant resonances, and other
laboratory observables. Global fits to experimental data have been
performed with Skyrme (e.g., Ref. <cit.>) and covariant
mean-field models (e.g., Ref. <cit.>). The second regime,
cold neutron matter below nuclear saturation density, is difficult to
probe experimentally but is well-constrained by theoretical
calculations based on semi-phenomenological nuclear forces or
microscopic chiral effective field theory-based
interactions <cit.>. The third
regime, strongly-interacting high-temperature matter, is best
described by interactions and many-body approaches similar to those
applied to cold neutron matter near saturation
density <cit.>.
The fourth regime, low-density and high-temperature matter that is
nearly non-degenerate, is best described by the virial expansion. The
equation of state in this regime is determined from nucleon-nucleon
scattering phase shifts <cit.>. Finally,
neutron-rich matter at densities above twice saturation density is
most strongly constrained by observations of neutron star masses and
radii, particularly the observation of neutron stars with M≃
2M_⊙ <cit.>.

There are several currently available tabulated equations of state.
The first set of EOSs, developed by Lattimer and Swesty
(LS) <cit.>, was constructed in the single-nucleus
approximation and based on three different non-relativistic Skyrme
interactions. Two of the three Skyrme interactions have nuclear
incompressibilities (K) far outside of modern constraints
<cit.>. The third, with K=220 MeV, has a
combination of symmetry energy (S) and slope of the symmetry energy
(L) that are only slightly outside of current constraints
(<cit.>). The LS EOS with K=220
MeV also produces a 2 M_⊙ neutron star and is still
important for simulations of core-collapse supernovae and neutron star
mergers. The second set of EOSs came from H. Shen et al. <cit.>
(also using the single-nucleus approximation) and was based on the NL3
relativistic mean-field Lagrangian. The values of K and L for NL3
are much larger than current neutron star obervations (see the
analysis in e.g., <cit.>) and nuclear
theory <cit.> suggest.

While the single-nucleus approximation is sufficient to describe the
bulk thermodynamics, it does not in general accurately describe the
composition <cit.>
and the associated weak reaction rates. More modern EOS tables often
include a more complete nuclear distribution as a result. The third
set from G. Shen et al. <cit.> includes an EOS table based on
a more modern relativistic mean-field model,
“FSUGold” <cit.>, and goes beyond the single nucleus
approximation to include a full distribution of nuclei in nuclear
statistical equilibrium (NSE). This model has values of K, L and
S that are within recent constraints from experiment and (in later
versions) produces a neutron star maximum mass larger than 2
M_⊙. These EOSs based on FSUGold include more modern
nuclear physics input, including a proper treatment of nearly
non-degenerate matter that matches the virial expansion. The fourth
set of EOS tables was based on the work of Hempel and
others <cit.> that built upon several nucleon-nucleon
interactions which produce reasonable values of K, S and L and
generate 2 M_⊙ neutron stars, including FSUGold,
DD2 <cit.>, IUFSU <cit.>, SFHo <cit.>
and SFHx <cit.>. The latter two interactions were designed
to simultaneously match laboratory nuclei and give neutron star radii
that match astronomical observations of neutron
stars <cit.>. More recently, several EOSs have been added
to the CompOSE (CompStar Online Supernovae Equations of State)
database <cit.>, including an EOS with
hyperons <cit.>.

In this work, we construct a phenomenological free energy density that
is consistent with observational and theoretical constraints in the
five aforementioned physical regimes. This is in contrast to works
which attempt to describe matter over the entire density and
temperature range with a single detailed model of the nucleon-nucleon
interaction. Many previous works proceed this way using a Skyrme-based
or relativistic mean-field model to describe matter at all densities
and temperatures. The principal problem is that these models are
guaranteed to work well only for isospin-symmetric nuclear matter at
zero temperature. Extrapolating these models to other density and
temperature regimes may lead to inaccurate EOS results or may
introduce unphysical correlations between the nature of matter across
different regimes. For example, given a Skyrme model it is common to
observe that the nuclear incompressibility is correlated with the
maximum mass of neutron stars. Such a correlation has little physical
meaning, however, since the neutron star maximum mass is determined by
interactions in high-density matter that likely have little similarity
to nucleons in the laboratory (see a similar argument in
Ref. <cit.>). We avoid extrapolations where possible, but
some extrapolation will still be required where experimental and
theoretical guidance is lacking.

Our second advance is in the treatment of uncertainties. The most
relevant parameters which describe the uncertainties in different
density and temperature regimes are not clearly related. The virial
expansion provides a clear path forward for describing uncertainties
at low-density and high-temperature, but higher-order virial terms are
not necessarily useful for quantifying uncertainties at higher
densities. In this work, through the construction of a
phenomenological model one can vary uncertainties in different regimes
independently, without spoiling agreement elsewhere.



§ METHOD


The EOS table is constructed by combining an EOS for homogenous
nucleonic matter consisting only of neutrons and protons. The EOS is
written in the form of the Helmholtz free energy (including only the
contribution from nucleons) f_np(n_B,Y_e,T). In the
discussion below, we remove the nucleon rest mass contributions from
the free energy densities and chemical potentials and use a tilde when
these rest mass contributions are included, i.e.

    f̃_np(n_B,Y_e,T)    ≡    f_np(n_B,Y_e,T)
      
           +
      n_B [ (1-Y_e) m_n + Y_e m_p
        ],

where m_n and m_p are the neutron and proton masses. When
electrons are included, their rest mass contribution to the free
energy is also included. We ignore muons and exotic charged particles
at higher densities. Thus the proton fraction, x_p, and electron
fraction, Y_e are always equal.



 §.§ Virial expansion and homogeneous nucleonic matter


The virial expansion is a model-independent way of computing the
pressure of matter at low densities and high
temperatures <cit.>. It is an expansion in
powers of the fugacity, z_i, of particle i defined by

    z_i = exp( μ_i,vir / T ) ,

where μ_i,vir denotes the nucleon chemical
potential. In matter consisting only of
neutrons and protons, the first-order terms in the pressure,
proportional to z_n and z_p, consist of the classical
non-interacting contribution to the pressure. The coefficients of the
second-order terms (second-order virial coefficients) in the virial
expansion can be obtained directly from nucleon-nucleon scattering
phase shifts. Third-order virial coefficients are not well known.
Nuclear statistical equilibrium implies that the fugacity of nuclei
can be written in terms of the neutron and proton fugacities

    z_(Z,N)∝ z_n^N z_p^Z  .

Thus when the neutron and proton fugacities are nearly equal the
contribution of deuterons comes at second order in the virial
expansion and the contribution from alpha particles comes at fourth
order. Two-body scattering between nucleons and alpha particles
contributes at fifth order in the virial expansion. In this work,
because the third-order virial coefficients are not well-known,
third-order and higher terms are ignored. 

Second-order terms in the virial expansion affect the description
of homogeneous nucleonic matter. In order to ensure that
the free energy matches the virial result at low densities
and high temperatures, the free energy density is written as

    f_np(n_B,x_p,T)    =    f_virial(n_n,x_p,T) g 
           + f_deg(n_B,x_p,T) (1-g)   ,

where f_virial is the virial free energy density,
f_deg is the free energy density when either
the neutrons or protons are sufficiently degenerate so
that the virial expansion is a poor approximation. The
function g is defined by

    g ≡ 1/(1+3z_n^2 + 3z_p^2)   .

This definition ensures that (1-g)f_deg appears as a
third- or higher-order correction to the free energy density in the
virial expansion as long as f_deg is at least linear in
the fugacity at low densities (we verify this below). The value of g
is 1 only when z_n and z_p are both sufficiently small. The
numerical coefficient 3 was chosen to ensure a positive entropy in
the entire region in (n_B,Y_e,T) space for which the pressure of our
EOS is positive. The relationship between the fugacities and the
densities is

    n_n    =    2 λ^-3[ z_n + 2 z_n^2 b_n(T)+ 2 z_n z_p
        b_pn(T)] 
    
      n_p    =    2 λ^-3[ z_p + 2 z_p^2 b_n(T) + 2 z_n z_p
        b_pn(T)]   .

These equations are solved for the fugacities in order to compute the
free energy density from the virial expansion. The quantity λ≡[4 π/(m_n T+m_p T)]^1/2 is the average nucleon
thermal wavelength, the quantity b_n(T) is the second neutron virial
coefficient, and the quantity b_pn(T) is the virial coefficient
describing the interaction between neutrons and protons.

The virial coefficients b_n(T) and b_pn(T) are determined by
scattering phase shifts. Analytical fits can be employed, similar to
those in Ref. <cit.>, but previous fits for b_n employ
functional forms which diverge for T→ 0. We perform an
alternate fit, constraining the zero temperature behavior to match
that expected from a finite-range expansion, which will be correct
when the density is sufficiently small. The values of the virial
coefficients at high temperature are not well-known, so we arbitrarily
constrain the fits so that the virial coefficients give the value
expected for noninteracting fermions at T=150 MeV. While very hot
and nearly nondegenerate matter is present in simulations, it is
unlikely to strongly affect the dynamics.

For the neutron matter virial coefficient, we use the data given in
Refs. <cit.> and add three points at T=0.1,
0.5 and 150 MeV, with virial coefficients of 0.207, 0.272, and
2^-5/2 respectively. The first two are determined from an
effective range expansion to the phase shift with scattering length
-18.9 fm and effective range 2.75 fm as determined from
Ref. <cit.>. The last value at T=150 MeV is the
non-interacting result. We fit this data to a 10-parameter functional
form

    b_n(T)    =     b_0 + b_1 T + b_2 T^2 + b_3 T^3 +
      b_4 e^-b_5 (T- b_6)^2
          + b_7 e^-b_8 (T- b_9)  .

We find that the parameter set b_0=0.28745,
b_1=2.2006 × 10^-3 MeV^-1, b_2= -2.6210  × 10^-5
MeV^-2, b_3= 6.0617  × 10^-8 MeV^-3, b_4= 1.0595
 × 10^-2, b_5= 5.6734  × 10^-2 MeV^-2, b_6=
3.4925 MeV, b_7=-2.7106  × 10^-3, b_8= 3.1405 MeV^-1,
b_9=1.2010 MeV matches the data. The data and the
fit are shown in the top panel of Fig. <ref>. 

The contribution from the deuteron binding energy is typically
included in b_pn(T), but in this work the deuteron binding energy
is removed (it will be added by the nuclear statistical equilibrium
part of the free energy in later work). For the low-temperature result
(T<1 MeV), both ^1 S_0 and ^3 S_1 phase shifts contribute
at low energy, while due to the factor e^-E/2T
from Eq. (22) in Ref. <cit.>, it is reasonable to
ignore the higher-order phase shift contributions. We use the
scattering length -23.74 fm and effective range 2.77
fm <cit.> for ^1 S_0 channel and the scattering length
5.418 fm and effective range 1.833 fm for the ^3 S_1 channel from
Ref. <cit.>. An alternate fit for b_pn is

    b_pn(T) = c_0 e^-c_1(T+c_2)^2 + c_3 e^-c_4 (T+ c_5),

where c_0=1.5273, c_1=1.7488  × 10^-4 MeV^-2,
c_2=1.7550  × 10^1 MeV, c_3= 0.45104, c_4= 0.27513
MeV^-1, c_5= -1.1250 MeV. The data and the fit are shown in the
bottom panel of Fig. <ref>.





 §.§ First derivatives of the free energy


From Eq. (<ref>), we can obtain the chemical
potentials and entropy directly

    μ_i     =    μ_i,vir g +
      f_virial∂g/∂ n_i
           + ∂ f_deg/∂ n_i(1-g)-f_deg∂ g/∂ n_i

for i=n,p 
where

    μ_i,vir≡∂ f_virial/∂ n_i

and for the entropy

    s    =    -∂ f_virial/∂ Tg-
      f_virial∂ g/∂ T
           -∂f_deg/∂ T(1-g )+f_deg∂ g/∂ T  .

In order to compute the derivatives of g with respect to the
densities, one can differentiate Eqs. (<ref>) with respect to
n_n and n_p and then solve the resulting four equations for the
quantities ∂μ_i,vir/∂ n_j (for
i,j=n,p).



 §.§ Matter near nuclear saturation density


Experimentally-measured nuclear masses are well-described by Skyrme
energy density functionals and thus it is expected that the energy
density of nuclear matter at zero temperature is also well-described
by the Skyrme model. In Ref. <cit.>, the parameters of
the Skyrme model were fit to several nuclear masses, charge radii, and
pairing energies using Bayesian inference. We use a set of 1000 Skyrme
parameterizations selected from the posterior distribution computed in
Ref. <cit.> to describe isospin symmetric matter.

Because there is relatively little information from theory or
experiment in some regimes, the Skyrme model will also be extrapolated
to higher densities and temperatures from below. This means that,
however, we cannot use all of the parameterizations because some of
them have a nucleon effective mass which becomes negative for
densities below n_B<2 fm^-3. We remove such Skyrme models
from consideration.

Nuclear mass measurements are restricted to relatively
isospin-symmetric nuclei, thus neutron matter is not necessarily
accurately described by Skyrme models (see e.g. the discussion
regarding large fluctuations in the isovector channel in
Ref. <cit.>). Zero-temperature neutron matter up to
nuclear saturation density is tractable in quantum Monte
Carlo <cit.> and many-body perturbation
theory <cit.>. It has thus become common to fit to
neutron matter calculations as well as nuclear mass
data <cit.>. However, this practice presumes that
the Skyrme functional is well-suited to describing pure neutron
matter, an assumption that is not necessarily valid. Thus, for pure
neutron matter we use the four-parameter expression based on quantum
Monte Carlo results from Ref. <cit.>,

    ε_QMC(n_B) =
      f_QMC(n_B)    
        = n_B[ a
        ( n_B/n_0)^α+ b (n_B/n_0)^β]      .

The range for the parameters 0.47 < a < 0.53 and 12 MeV <
α < 13 MeV is chosen as in Ref. <cit.> to
enclose the limits in Ref. <cit.>.

The symmetry energy implied by many of the Skyrme fits, when combined
with the quantum Monte Carlo results for neutron matter, naturally
implies bound neutron matter at subsaturation densities. Similarly,
much of the range for S and L implied by the Skyrme
parameterizations is outside the allowed range from
Ref. <cit.>. Thus we ignore the values for S and L from the
Skyrme models and limit L between 44 and 65 MeV, and S between
29.5 and 36.1 MeV as in Ref. <cit.>. These bounds are
consistent with recent microscopic constraints <cit.> on the
density dependence of the symmetry energy from chiral effective field
theory. The prescription (9.17 S - 266 MeV) < L < (14.3 S
- 379 MeV) is used to ensure that S and L are
correlated. The coefficients b and β are determined by

    b= S - a + (E/A)_sky


    β=1/b(L/3 - α a)

where (E/A)_sky is the binding energy per
particle of nuclear matter from the Skyrme interaction. Finally, we
combine the nuclear matter and neutron matter results by assuming the
symmetry energy is quadratic in x_p. This choice ensures that
nuclear matter is representative of experimental results on nuclear
masses while neutron matter agrees with modern theory results.

We note that the free energy density of matter from the QMC results
above is always at least linear in the density, and by
Eq. (<ref>) at least linear in the fugacity at low densities.
This also holds for the Skyrme model, since the kinetic part of the
energy density is proportional to k_F^5 and the potential energy
part is proportional to at least one power of the density. Thus our
function g in Eq. (<ref>) above is defined so that
f_deg will leave the second-order virial coefficients
unchanged from the values determined by experiment in the low-density
limit.



 §.§ Matter at high densities


Above nuclear saturation density, there are two principal sources for
constraints on the EOS of matter: heavy ion collisions and neutron
star observations. Constraints on the EOS from heavy ion collisions
near the saturation density do not yet contradict results from Skyrme
fits. On the other hand, constraints from heavy ion collisions on the
EOS at higher densities do not yet provide a clear picture. Until the
results from heavy ion collisions are more definitive, the Skyrme
model from above is extrapolated to higher densities to describe
isospin-symmetric nuclear matter.

Neutron star mass and radius observations constrain the equation of
state of neutron-rich matter at high densities, in particular, the
pressure as a function of the energy density <cit.>.
Unfortunately, neutron star observations do not yet currently
constrain the proton fraction of neutron star matter.
We find that the form

    ε_NS(n_B)
      = f_NS(n_B)    =    p_0 n_B √(n_B) +
      p_1 n_B^2 + p_2 n_B^2 √(n_B)
           + p_3 n_B^3 + p_4 n_B^4

provides a good fit to the results from Ref. <cit.>. We
randomly select EOSs from a Markov chain constructed in
Ref. <cit.>, and fit them to Eq. (<ref>).
Ref. <cit.> constructed several Markov chains, and we use
the chain which was constructed using GCR (as above; from
Ref. <cit.>), Model A (which models high-density matter
using polytropes), and includes all of the mass and radius data from
photospheric radius expansion X-ray bursts and quiescent low-mass
X-ray binaries.

Astrophysical simulations can probe densities larger than those
constrained by the neutron star data in Ref. <cit.>.
Between a transition density, n_Bf and the largest
baryon density we consider, n_B = 2 fm^-3, we
implement a simple EOS adapted from Constantinou and
Prakash <cit.>.
The procedure for matching these EOSs begins by 
setting n_Bf equal to the highest density specified
by the Monte Carlo data in Ref. <cit.>. We decrease
this transition density as necessary to ensure that the EOS from
Eq. (<ref>) is causal for densities lower than this
transition density. We add an additional parameter, ϕ,
which is equal to the speed of sound at the largest density we
consider, n_B = 2 fm^-3.

The EOS between n_B=n_Bf and n_B = 2 fm^-3 is chosen
depending on the relative magnitude of the speed of sound at these two
endpoints. If the speed of sound is increasing with increasing baryon
density, then, we choose

    c_s^2 = 1-a_1 +a_1 a_2 n_B^a_1/1+a_2 n_B^a_1

and determine a_1 and a_2 by matching the
boundaries ≡ c_s^2 (n_B f) and
ϕ=c_s^2 ( n_B=2 fm^-3) thus
ensuring β→ 1 as n_B →∞.
The energy density above n_B = n_Bf is 

    ε_NS = -m_n n_B + c_1 ( 1/2 a_2 n_B^2+ n_B^2-a_1/2-a_1) + c_2,

where

    c_1    =   ε_f+m_n n_Bf+ P_f/n_Bf^2
        (a_2+n_Bf^-a_1)
    
      c_2    =   1/2[(ε_f +m_n n_Bf
        -P_f )+a_1 ε_f +m_n n_Bf+P_f/(a_1-2)(1+a_2 n_Bf^a_1)]

Alternatively, if the speed of sound is decreasing with increasing
density, then we set

    c_s^2 = a_1 - a_1 a_2 n_B^a_1/1+a_2 n_B^a_1

and match the boundaries as before
ensuring β→ 0 as n_B →∞. The corresponding energy density is 

    ε_NS = c_1 n_B  _2 F_1 (1, -1/a_1,
        1-1/a_1, -n_B^-a_1/a_2)/a_2+c_2 - m_n n_B

where _2 F_1 is a hyper-geometric function with Pfaff's
transformation and the constants c_1 and c_2 are

    c_1    =    n_B_f^-a_1-1(a_2 n_B_f^a_1+1)
      (ε_f + m_n n_B_f
      + P_f )
      
    
      c_2    =     n_B_f^-a_1[a_2 n_B_f^a_1(ε_f + m_n n_B_f) .
        
           - (a_2 n_B_f^a_1 +1) _2 F_1 (1, -1/a_1,
        1-1/a_1, -n_B^-a_1/a_2)
        
          . (ε_f + m_n n_B_f + 
        P_f ) ]
      1/a_2  .

Although in practice ϕ is chosen randomly so this is
rare, if c_s^2 (n_B f)=
ϕ=c_s^2 ( n_B=2 fm^-3), then we ensure
c_s^2 is constant at high densities.
The corresponding energy density is 

    ε_NS   =    -m_n n_B + (ε_f+
        m_n n_B_f + P_f)/(1+C_s_f^2)( n_B/n_B_f)^1+C_s_f^2
           + C_s_f^2 ( ε_f +
        m_n n_B_f)-P_f/1+C_s_f^2.

This speed of sound correction ensures that neutron star matter is
causal, but an additional correction (described below) will be
required to ensure that the speed of sound is not larger than the
speed of light at all temperatures and electron fractions.

In order to combine information from QMC near the saturation
density and information from neutron star observations at
higher densities, we define a function h,

    h = 1/1+expγ(n_B - 3/2 n_0 )

where γ is 20.0 fm^3. This function is
used to interpolate between the two density regimes.



 §.§ Hot matter near the saturation density


Nuclear two- and three-body forces based on chiral effective theory have
shown great progress in computing the EOS of matter near nuclear
saturation. The Kohn-Luttinger-Ward perturbation series can be used to
compute the EOS of matter at finite temperature as described in
Refs. <cit.>. The resulting EOS can then be 
fitted with a Skyrme interaction, as done for example in Ref. <cit.>.
However, it is difficult to use these results to 
quantify the uncertainties in these EOS calculations for matter at
T=0 where large cancellations between attractive and repulsive 
interactions lead to large theoretical errors.

To attempt to address this, we refit only the finite-temperature
correction from the chiral EOS,

    Δ f_hot(n_B,x_p,T) ≡
      f_hot(n_B,x_p,T) - f_hot(n_B,x_p,T=0)

and add these finite temperature corrections on top of our EOS. The
EOS for neutron matter (x_p=0) and nuclear matter (x_p=1/2) is
obtained from the perturbation series (including the non-interacting
contribution) and fitted with a single Skyrme model. The resulting
parameter set is given in Table <ref>. We assume that these
finite temperature corrections are quadratic in the isospin asymmetry,
δ. The EOS is not fully
quadratic <cit.>, but the quadratic
approximation is good enough in comparison to the uncertainties in the
nature of the strong interaction above the saturation density. We do
not expect this Skyrme interaction to give a reasonable decription of
nuclei or saturated nuclear matter, because we only employ it to
describe the finite temperature part of the EOS. There are some
regions, especially at large densities, for which the EOS is unstable,
i.e. ds/dT < 0, but these regions most often result in an acausal
EOS and are thus fixed by the speed of sound correction described
below.





 §.§ The full combined EOS


First, we define the symmetry energy to
include a zero-temperature contribution which combines
the QMC EOS near saturation density, the neutron star fit at higher
densities, and the Skyrme interaction for isospin-symmetric matter

    ε_sym(n_B)    =   
      h (n_B) ε_QMC(n_B) +
      [ 1-h(n_B) ] ε_NS(n_B)
      
           - f_Skyrme(n_B,x_p=1/2,T=0)   .

Defining the isospin asymmetry δ=1-2 x_p, we can combine this
with the model described above to obtain the free energy density of
degenerate matter

  
    f_deg(n_B,x_p,T)    =   
        f_Skyrme(n_B,x_p=1/2,T=0)
        + δ^2 ε_sym(n_B)
        + 
           + δ^2 Δ f_hot(n_B,x_p=0,T)
        + (1-δ^2) Δ f_hot(n_B,x_p=1/2,T).


Finally, we ensure that the total nucleonic
free energy gives the result from the virial expansion at
high temperatures using Eq. (<ref>).
When we need to include the electrons, positrons,
and photons, we define the free energy density

    f_npeγ≡ f_np +f_e^-+f_e^++f_γ  .
 

Using this formalism, the chemical potentials and entropy
can be computed directly:

    ∂ f_deg/∂ n_n   =   1/2μ_n,Skyrme(n_B,x_p=1/2,T=0)
     
          + 1/2μ_p,Skyrme(n_B,x_p=1/2,T=0)
      
           + δ^2
      ∂ε_sym/∂ n_B +
    2 δ (1-δ)/n_Bε_sym
          
      + 2 δ (1-δ)/n_BΔ f_hot(n_B,x_p=0,T)
      
            + δ^2 Δμ_n,hot(n_B,x_p=0,T) 
      
          
      - 2 δ (1-δ)/n_BΔ f_hot(n_B,x_p=1/2,T) 
      
            +  ( 1-δ^2 ) Δμ_n,hot
        (n_B,x_p=1/2,T),


    ∂ f_deg/∂ n_p   =   1/2μ_p,Skyrme(n_B,x_p=1/2,T=0)
      
            +1/2μ_n,Skyrme(n_B,x_p=1/2,T=0)
      
           + δ^2
      ∂ε_sym/∂ n_B -
      2 δ( 1 + δ)/n_Bε_sym
          
      - 2 δ( 1 + δ)/n_BΔ f_hot(n_B,x_p=0,T) 
      
            + δ^2 Δμ_p,hot(n_B,x_p=0,T) 
      
          
      + 2 δ( 1 + δ)/n_BΔ f_hot(n_B,x_p=1/2,T) 
      
            + ( 1-δ^2 ) Δμ_p,hot
        (n_B,x_p=1/2,T),

and

    ∂ f_deg/∂ T   =   
      -  δ^2 s_hot(n_B, x_p=0,T ) - 
      
          ( 1-δ^2 ) 
      s_hot(n_B, x_p=1/2,T ),

where

    ∂ε_sym/∂ n_B   =   
      h^'(n_B) ε_QMC(n_B)
      + h(n_B) ε^'_QMC(n_B)
      
          
      - h^'(n_B) ε_NS(n_B)
      + [1-h(n_B)] ε^'_NS(n_B) +
      
          
      - 1/2[μ_n,Skyrme(n_B,x_p=1/2,T=0) .
        
          
        + . μ_p,Skyrme(n_B,x_p=1/2,T=0) ].


In summary, we have 5 parameters: (i,ii) the values of a and
α which determine sub-saturation neutron matter, (iii,iv) the
values of S and L which determine the symmetry energy and its
density dependence, and (v) the value, ϕ, of the speed of sound
in neutron star matter at n_B=2 fm^-3. In addition, we
have two indexes which enumerate random samples from posterior
distributions including (i) the index of the Skyrme parameterization
from Ref. <cit.> and (vi) the index of the neutron star
EOS from the Markov chain generated in Ref. <cit.>.



 §.§ Enforcing causality at high density



Since our phenomenological EOS does not have manifest Lorentz
covariance, it has the potential to become acausal at high-densities.
At every electron fraction and temperature, there may be a baryon
density, n_B^*, above which the EOS becomes acausal.
Because our phenomenological EOS (as all other EOS tables) operate
as functions of the densities and temperatures, it is useful to
rewrite the speed of sound in terms of derivatives of the
Helmholtz free energy. This is done in Appendix I below for a
general system with any number of conserved charges (though here we
only have two, baryon number and electric charge). 

When our phenomenological EOS becomes acausal above some baryon
density, n_B^*, we replace the EOS with a causal EOS,
f̃_C&P, following the prescription in
Ref. <cit.>. We construct a modified free energy
density with the following

    f̃_all=f̃_npeγΘ(n_B^*-n_B)+
      f̃_C&PΘ(n_B-n_B^*),

where contributions from electrons, positrons and photons are included
in f̃_npeγ. To be more concise, we suppress the
subscripts npeγ in the following. Using ε for
energy density (including rest mass energy density), S for entropy,
s for entropy density, and s̃ for entropy per baryon, the
C&P speed of sound is

    c_s^2 = ( dP/dε)_s̃,N_B,N_e
    = ( dP/dε)_s̃,N_B,Y_e.

Note that s̃=S/N_B = s/n_B, where N_B is the number of
baryons. The C&P derivation begins by noting that

    P = - ε + 
    n_B (∂ε/∂
    	n_B)_s̃,N_B,Y_e  .

To see this we can write

    (∂ε/∂
    	n_B)_s̃,N_B,Y_e      =
    [∂ (E/V)/∂
    	V]_S,N_B,Y_e[∂ (N_B/V)/∂
    	V]_S,N_B,Y_e^-1
          =(-P/V - E/V^2)
    (-N_B/V^2)^-1 = (P+ε)/n_B.

Taking the derivative of Eq. (<ref>), we can also rewrite the
pressure as a second derivative

    (∂ P/∂ n_B)_s̃,N_B,Y_e   =   
    -(∂ε/∂
    	n_B)_s̃,N_B,Y_e +
    (∂ε/∂
    	n_B)_s̃,N_B,Y_e
          + n_B
    (∂^2 ε/∂
    	n_B^2)_s̃,N_B,Y_e
       =   n_B ( ∂^2 ε/∂ n_B^2)_s̃,N_B,Y_e.

Thus we can proceed as C&P do,

    ( ∂^2 ε/∂ n_B^2)_s̃,N_B,Y_e - c_s^2/n_B( ∂ε/∂ n_B)_s̃,N_B,Y_e = 0.


Following the analytical continuation, at every value of s̃,
N_B, and Y_e in the acausal region, we can use the C&P solution

    ε_C&P(s̃,n_B,N_B,Y_e)       =
    [ε^*(s̃,N_B,Y_e)+P^*(s̃,N_B,Y_e)/β+1]
    
          ×[ n_B/n_B^*
        (s̃,N_B,N_e)]^β+1
          +[
    βε^*(s̃,N_B,Y_e) - P^*(s̃,N_B,Y_e) /β+1]   .

If we assume that ε^*, P^*, and n_B^* are
volume independent, then since they are determined at a fixed
value of n_B,
they cannot separately depend on N_B. Thus the full energy density
also does not depend on N_B.
To show this explicitly, we start from

    E   =   μ_B N_B+μ_L N_B Y_e +T S - P V,
    
    
    dE   =   ( μ_B+μ_L Y_e +T s̃) dN_B +μ_L N_B dY_e
    
          +T N_B ds̃ -P dV.

where μ_L ≡μ_p + μ_e - μ_n. On the right hand side,

    dV=1/n_BdN_B-N_B/n_B^2 d n_B

and on the left hand side,

    dE= d(εN_B/n_B)
    =N_B/n_B dε+
    ε/n_B dN_B -ε N_B/n_B^2 dn_B.

Substituting the above two equations back into Eq. (<ref>),
we find

    dε= μ_L n_B dY_e
    +T n_B ds̃ + (μ_B + μ_L Y_e + T s̃) d n_B

and the N_B dependence disappears. Therefore 
ε_C&P(s̃,n_B,N_B,Y_e)=
ε_C&P(s̃,n_B, Y_e):

    ε_C&P(s̃,n_B,Y_e)    =   [ε^*(s̃,Y_e)+P^*(s̃,Y_e)/β+1]
    [ n_B/n_B^*(s̃,Y_e)]^β+1
           +[
    βε^*(s̃,Y_e) - P^*(s̃,Y_e) /β+1]   .

We also find 

    ( ∂ε_C&P/∂ Y_e)_s̃,n_B=μ_L n_B,


    ( ∂ε_C&P/∂s̃)_n_B,Y_e=T n_B,


    ( ∂ε_C&P/∂ n_B)_s̃,Y_e=μ_B+μ_L Y_e + T s̃.

We choose to use this solution above the value of n_B, denoted
n_B^*(s̃,Y_e), for which c_s^2 = β.
From Eq. (<ref>) the temperature is 

    T = 1/n_B(∂ε_C&P/∂s̃)_n_B,Y_e.

Thus computing the temperature using the method in
Ref. <cit.> requires evaluating derivatives of the form
(∂ n_B^*/∂s̃)_Y_e. These derivatives are
computed along the surface for which the speed of sound is equal to
its largest possible value (which we set to 0.9). Because the speed of
sound requires two derivatives of the free energy, these derivatives
(derivatives along a line of constant c_s^2) require third
derivatives of the free energy. The chemical potentials require a
similar calculation based on Eqs. (<ref>) and (<ref>) 
which we do not explicitly show here.

The numerical evaluation of third derivatives introduces quite a bit of
noise. Our phenomenological formalism ensures that these derivatives
can be computed analytically. We leave this calculation to future
work. In the mean time, we can more easily compute the temperature
by an indirect approach. Presuming we would like to compute the
EOS at some fixed values of the baryon density n_B,
electron fraction Y_e,
and temperature T, one needs to solve the five
equations

    c_s,DS^2(n_B1^*,Y_e,T_1^*)    =   β  , 
    
    
    c_s,DS^2(n_B2^*,Y_e,T_2^*)    =   β  ,
    
    s̃_DS( n_B1^*,Y_e,
          T_1^*)    =   s̃  ,
    
    s̃_DS( n_B2^*,Y_e,
          T_2^*)    =   s̃ + δs̃  ,
     
    1/n_Bϵ_C&P (s̃ + δs̃,
      Y_e,n_B)
      - ϵ_C&P (s̃, Y_e,n_B)/δs̃   =   T

for the five values n_B1^*, n_B2^*, T_1^*, T_2^*, and
s̃ given some small fixed stepsize in the entropy per baryon,
δs̃. Here “DS” denotes our EOS while “C&P” denotes
Eq. (<ref>) from Ref. <cit.>. The required
numerical derivatives to compute the temperature in the last equation
of Eq. (<ref>) can then be computed directly from finite
differences. By solving the five equations above, we get
s̃_C&P( n_B,Y_e, T), and then

    f̃_C&P(n_B,Y_e, T) =
    ε_C&P - Tn_Bs̃_C&P( n_B,Y_e, T).




§ RESULTS




Fig. <ref> shows how our full EOS behaves in the (i)
non-degenerate limit, (ii) the limit of zero-temperature neutron
matter, and (iii) the limit of high-temperature and high-density
nuclear matter. Only one parameter set is chosen and the same
parameter set is chosen for each of the three panels. In the
non-degenerate limit, our full EOS smoothly matches on to the virial
EOS as determined by Eq. (<ref>). The middle panel shows that, in
the limit of zero-temperature neutron matter, our full EOS matches the
QMC result at low density and remains close to the neutron star EOS at
moderate densities. At higher densities, the free energy does not
increase too quickly with density because our correction for causality
begins to start becoming important. The bottom panel compares
our full EOS (for this parameterization) with the (T=0) Skyrme EOS,
showing a small modification in the EOS due to the finite temperature
correction from the chiral EOS.

We can construct a figure similar to Fig. <ref> for any
physical combination of our model parameters: α, a, S, L,
ϕ and Skyrme model from Ref. <cit.> (which
represents a 12-dimensional space of possible Skyrme models) and any
neutron star model from Ref. <cit.> (which represents a
6-dimensional space of high-density EOSs).



For any baryon density, electron fraction, and temperature, we can
compute a probability distribution for the free energy per baryon over
our entire parameter space. The magnitude and shape of the uncertainty
in the free energy changes depending on the relevant physics in that
region which our model describes. Fig. <ref> (??The
titles on Figs. 3-7 should have N_B → n_B, and I would
remove the “/197.327” in the temperature??) shows the variation in
the free energy in nondegenerate matter where the EOS is dominated by
the virial expansion. The small remaining uncertainty here originates
in the Skyrme model that is selected, and the non-smooth nature in the
probability distribution is a relic of the limited sampling size of
Skyrme interactions.



Fig. <ref> shows the variation in the free energy per
baryon at nuclear saturation density in nearly pure neutron matter and
in the limit of zero temperature. The results for LS220, SFHo, SFHx,
and IUFSU are also shown. The distribution is centered around 16 MeV,
corresponding to a symmetry energy of 32 MeV, and values lower than 13
MeV are excluded in our model because they correspond to symmetry
energies lower than 29.5 MeV. Symmetry energies this small seem to be
excluded from Quantum Monte Carlo and chiral effective field theory
calculations of pure neutron
matter <cit.>. See also
Ref. <cit.> for a more general result that gives a similar
lower limit for S of 28 MeV.



Fig. <ref> shows the distribution in the free energy
per baryon at a slightly larger temperature, and the entropy
contribution drops the free energy per baryon in each case. In
comparison to the results from Fig. <ref>, the LS220
free energy per particle drops more than SFHo or SFHx because its
effective mass (equal to the nucleon mass) is larger than that in
SFHo/x (about 0.7 times the nucleon mass). The effective mass can be
computed from the chiral interaction directly
<cit.>, and close to the Fermi surface it is
found to be nearly equal to the free-space nucleon mass. Second-order
many-body perturbation theory contributions, however, produce a strong
momentum dependence in this region, and averaging around k=k_F the
effective mass is about 0.85 times the free-space nucleon mass. The
distribution in Fig. <ref> is thus larger than LS220
because of the smaller symmetry energy (S=28.6 MeV) and larger
effective mass in the LS220 EOS.



Fig. <ref> shows the probability distribution for the
low-temperature neutron-rich matter free energy per baryon
at higher densities. There is clearly a much larger spread in the
free energy per baryon, corresponding to our larger ignorance
regarding the nature of matter at higher density. While all models
LS220, SFHo, SFHx and IUFSU are inside the region suggested by
our parameterization, our distribution leans towards
smaller values of the free energy because of the constraint
from relatively small neutron star radii.



Fig. <ref> shows low-temperature nuclear matter at
higher densities. The non-smooth nature of the distribution is due to
the small statistics afforded by the limited number of Skyrme models
we have employed. IUFSU suggests a larger free energy here because it
originates in a relativistic mean field model which tends to give
larger pressures than the non-relativistic models like Skyrme. This
region of parameter space is almost entirely unconstrained by
experiment, since it is not possible to make cold isospin symmetric
matter at this density. However, dense isospin-symmetric matter is not
as relevant for this work since simulations are typically neutron-rich
at high density.



Fig. <ref> shows the free energy per
baryon for one of our parameterizations as a function of baryon
density and temperature for two electron fractions. Using the
formalism presented in this work, thousands of similar density
plots can be generated with alternate parameterizations. The
largest variation between parameterizations is in the
free energy per particle at high density which changes with
the “stiffness” of the EOS.



Since accurate neutron star radii have not yet been measured for large
mass neutron stars, the speed of sound of matter at the highest
densities probed in supernova and merger simulations is not
constrained by experiment. We have parameterized this variation with
ϕ. However, the speed of sound must increase at moderate
densities in order to reproduce the observation of a two solar mass
neutron star. Fig. <ref> shows the behavior of the speed of
sound in neutron-rich matter between n_B=0.1 fm^-3 and
n_B=2 fm^-3 and demonstrates these two regimes. The
speed of sound must increase quickly below 1 fm^-3 to
ensure that the neutron star maximum mass is sufficiently large, and
the speed of sound at higher densities varies considerably depending
on the value of ϕ. We restrict the speed of sound to be less than
c√(0.9) to ensure finite-precision errors in simulations do
not create unphysical sound speeds.



In high-density isospin-symmetric matter, the speed of sound is
dominated by the Skyrme model used for isospin-symmetric matter near
the saturation density. This region is principally constrained to have
a speed of sound less than c√(0.9) by our implementation of the
prescription from Ref. <cit.> as described in
section <ref>. This is shown in Fig. <ref>. There is
also a slight residual impact from the modification in the speed of
sound from ϕ at values of Y_e which are nearly but not exactly
equal to 1/2, so there are some slight kinks in the curves in
Fig. <ref> near n_B=1 fm^-3.

At sufficiently high density, the entropy from the chiral EOS begins
decreasing with increasing temperature. Because of the presence of
ds/dT in the speed of sound (see f_TT in the denominator of
Eq. (<ref>)), this unstable region implies a large speed of
sound. Our use of the C&P prescription thus cures this instability in
the extrapolated form of the finite-temperature corrections from the
chiral EOS. This is demonstrated in Fig. <ref>, which shows
contours of fixed ds/dT. This derivative becomes negative in the
upper right region, but this is always at a density larger than
n_B^* where the C&P EOS takes over.



As a final demonstration that our implementation of the C&P
prescription generates a continuous EOS, we show the entropy as a
function of density for several fixed temperatures. Our EOS above is
used for densities below n_B^* (indicated by the red dot) and the
C&P EOS is used for densities above n_B^*. It appears in this
plot that n_B^* is independent of temperature, but this is not
exactly true as there is a weak temperature dependence as shown in
Fig. <ref>. We have found that solving Eqs. (<ref>) is
numerically challenging because of the numerical derivatives involved
in computing the speed of sound (we use exact expressions for the
entropy and chemical potentials but numerical differentiation for the
number susceptibilities and other second derivatives of the free
energy). The combination of the numerical derivatives plus the
numerical noise in the Newton-Raphson method used to solve
Eqs. (<ref>) leads to a bit of noise in the entropy at large
densities. Future work will use exact expressions for second
derivatives of the free energy and thus facilitate the correction to
the speed of sound.





§ DISCUSSION


While we have made an attempt to explore the allowed EOS space as much
as possible, there are several regions in which our parameterization
is limited. Variations in the functions g and h could be explored,
but modifications of these functions cannot be too large or they are
likely to create small regions where the EOS has an acausal speed of
sound. A better quantification of the uncertainties in the finite
temperature part of the chiral EOS will be performed in future work.
We have also naively extrapolated the Skyrme model used in
isospin-symmetric matter at zero temperature near the saturation
density to higher densities. There is some experimental constraint on
matter in this region from heavy ion collisions <cit.>, but
the associated systematic uncertainties are not well-understood. In
any case, dense and isospin-symmetric matter is not often explored in
the neutron-rich matter encountered in neutron stars.

Uncertainties in the EOS may be relevant for the timescale for
core-collapse supernovae to explode after bounce <cit.> and
also, for example, in determining the amount of r-process
nucleosynthesis which occurs in the neutrino-driven
wind <cit.>. EOS uncertainties are also relevant for
neutron star mergers, as they dictate the lifetime of hypermassive
neutron star remnants and also the amount of r-process material
ejected <cit.>.

Because we use the Markov chain from “Model A” in
Ref. <cit.>, our EOS specifically prefers more moderate
phase transitions, which is appropriate to our assumption that matter
consists only of nucleons and no exotic matter. Thus our
uncertainties at high-density may be underestimated because strong
phase transitions are disfavored. One of the advantages of our
analytical form for the EOS is that our work can be easily generalized
to an EOS which includes exotic matter at high densities if desired.

This article is the first step towards a full quantification of how
microphysical uncertainties may affect core-collapse supernovae and
neutron star mergers. The next step is a full description of nuclei in
the dense matter environment with uncertainties that properly reflect
the relationship between nuclear structure and the underlying
nucleon-nucleon interaction. One way to include nuclei on top of our
EOS for homogeneous matter is to use the framework recently developed
in Ref. <cit.>. Finally, the EOS uncertainties must be
propagated through to the neutrino opacities. As this uncertainty
quantification matures, the comparison of simulations with data
points, such as GW170817 <cit.> and future nearby
core-collapse supernovae will provide more insight into what models
might be ruled out.



§ ACKNOWLEDGEMENTS


The authors would like to thank D. Higdon for suppling the samples
from the posterior generated in Ref. <cit.>. XD and AWS
were supported by DOE SciDAC grant DE-SC0018232 and NSF grant PHY
1554876, and JWH was supported by NSF grant PHY 1652199. Portions of
this research were conducted with the advanced computing resources
provided by Texas A&M High Performance Research Computing.



§ APPENDIX I - SPEED OF SOUND FOR A MULTICOMPONENT SYSTEM


Using ε for energy density, S for entropy,
s for entropy density, and s̃
for entropy per baryon, and assuming neutrinos are not trapped,
the speed of sound is (all chemical potentials and energy densities
below include the rest mass contribution even though not explicitly
indicated)

    c_s^2 = ( ∂ P/∂ε)_s̃,{ N_i }  .

In infinite matter, it is useful to rewrite this derivative in
terms of fixed volume rather than fixed number. 

    c_s^2 = ( ∂ P/∂ε)_S,{ N_i } =
    ( ∂ P/∂ V)_S,{ N_i }( ∂ε/∂ V)_S,{ N_i }^-1,

where the second derivative on the right-hand-side of this expression is

    ( ∂ε/∂ V)_S,{ N_i }   =   [ ∂  (E/V)/∂ V]_S,{ N_i } =
    -1/V P - E/V^2
       =    - P+ε/V
    = - T s + ∑_i μ_i n_i/V

and the first derivative on the right-hand side is

    ( ∂ P/∂ V)_S,{ N_j }   =    -
    ( ∂ε/∂ V)_S,{ N_j} +
    S [ ∂ (T/V)/∂ V]_S,{ N_j } +
    
          ∑_i 
    N_i [ ∂  (μ_i/V)/∂ V]_S,{ N_j }
       =    -
    ( ∂ε/∂ V)_S,{ N_j }
           +S [ -T/V^2 + ( ∂ T/∂ V)_S,{ N_j }] 
    
           + ∑_i 
     N_i [ -μ_i/V^2 +
    ( ∂μ_i/∂ V)_S,{ N_j }]
    
       =   P + ε/V +
    S [ -T/V^2 - ( ∂ P/∂ S)_{N_j},V] 
    
           +∑_i N_i [ -μ_i/V^2 -
    ( ∂ P/∂ N_i)_S,{N_j≠ i},V] 
       =    - S ( ∂ P/∂ S)_{n_j},V
    - ∑_i N_i ( ∂ P/∂ N_i)_S,{n_j≠ i},V

Putting these two results together gives

    c_s^2    =   [s ( ∂ P/∂ s)_{n_j},V +
    ∑_i n_i ( ∂ P/∂ n_i)_S,{n_j≠ i},V]
    
          ( T s + ∑_i μ_i n_i )^-1.

To re-express this in terms of derivatives of the free energy,

    c_s^2    =   {s [ ∂ (∑_i μ_i n_i - f)/∂ s]_{n_j},V.
    
          . +∑_i 
    n_i[ ∂  ( ∑_k μ_k n_k - f)/∂ n_i]_s,{n_j≠ i},V}(
    T s + ∑_i μ_i n_i )^-1.

For the sum over k,
all densities are constant except for n_i, thus 

    ∑_i 
    n_i ∂/∂ n_i( ∑_k μ_k n_k - f )_s,{n_j≠ i},V
       =   ∑_i n_i ∂/∂ n_i( ∑_k≠ iμ_k n_k + μ_i n_i -f
    )_s,{n_j≠ i},V
       =   ∑_i [ ∑_k n_k (∂μ_k /∂ n_i)_s,{n_j≠ i},V + μ_i .
    
          . -(∂ f/∂ n_i)_s,{n_j≠ i},V].


To compute this we need

    (∂ f/∂ n_i)_s,{n_j≠ i},V   =   (∂ f/∂ n_i)_{n_j≠ i},T,V
           +(∂ f/∂ T)_n_B,{n_j≠ i},V(∂ T/∂ n_i)_{n_j≠ i},s,V
       =   μ_i - s (∂ T/∂ n_i)_{n_j≠ i},s,V
    (∂μ_k/∂ n_i)_s,{n_j≠ i},V   =   (∂μ_k/∂ n_i)_{n_j≠ i},T,V
           +(∂μ_k/∂ T)_n_i,{n_j≠ i},V(∂ T/∂ n_i)_{n_j≠ i},s,V
           = f_n_i n_k + f_n_k T(∂ T/∂ n_i)_{n_j≠ i},s,V

which requires

    (∂ T/∂ n_i)_{n_j≠ i},s,V   =    -(∂ s/∂ n_i)_{n_j≠ i},T,V(∂ s/∂ T)_{n},V^-1
       =    -f_n_i T/f_TT

Finally, we get

    c_s^2    =   {
    - (s/f_TT) ( ∑_i n_i f_n_i T+s ) .
    
          . + ∑_i n_i [ ∑_k n_k (f_n_i n_k- f_n_k T
    f_n_i T f_TT^-1) 
    - s f_n_i T f_TT^-1]
    }
          ( T s + ∑_i μ_i n_i )^-1
       =   [
    ∑_i ∑_k n_i n_k (f_n_i n_k- f_n_k T
    f_n_i T f_TT^-1) .
    
           .- 2∑_i s n_i f_n_i T f_TT^-1
    - s^2 f_TT^-1] (
    T s + ∑_i μ_i n_i )^-1.




 §.§ Appendix I - Speed of Sound


It is important to ensure that the speed of sound at high densities is
smaller than the speed of light. High-density matter is nearly in
beta-equilibrium and has a nearly constant entropy per baryon
(s̃≡ S/N_B=s/n_B). We also choose to parameterize
high-density matter in terms of the chemical potential for electron
lepton number, μ_L, where cold beta-equilibrium corresponds to
μ_L=0. (When neutrinos are trapped, μ_L varies with density
and can be as large as 200 MeV <cit.>.) The relevant
thermodynamic derivative is

    c_s^2 = (dP/d ε)_S,μ_L,N_B

(since the baryon number, N_B, is constant, constant S is
equivalent to constant s̃). For fixed s̃ and fixed
μ_L, the phenomenological EOS becomes acausal above some baryon
density, n^*. Above this density, we replace the EOS with a causal
EOS following the prescription in <cit.>

    f_all=f_npeγΘ(n_B^*-n)+ f_causalΘ(n-n_B^*)


In the results below, we need to include the electrons, positions,
and photons, so we define

    f_npeγ≡ f_total +f_e^-+f_e^++f_γ

so our equation becomes 

    f_npeγ   =   μ_n n_n+ μ_p n_p+ μ_e^- n_e^-+ μ_e^+ n_e^+
          -P_npeγ

Define μ_B=μ_n, n_e=n_e^--n_e^+,
μ_Q=μ_n-μ_p, μ_L=-μ_Q+μ_e.
We have μ_e^-=-μ_e^+. Through lepton conservation
and electrical neutrality n_p=n_e^+-n_e^-,
we obtain

    f_npeγ=μ_B n_B+ μ_L n_e-P_npeγ

with μ_L being the out-of-whackness vector.

Now the speed of sound can expressed as a series of elementary
arithmetic of thermodynamics derivatives.

    C_s^2    =   ( ∂ P/∂ϵ)_s̃, Y_e N_B
          =∂(P,s̃,Y_e,N_B )/∂(V,T,Y_e,N_B) /
      ∂(ϵ,s̃,Y_e,N_B)/∂(V,T,Y_e,N_B )

The numerator can be computed through a series of Jacobian
transformations

    ∂(P,s̃,Y_e,N_B )/∂(V,T,Y_e,N_B)   =   ( ∂ P/∂ V)_T, Y_e, N_B( ∂s̃/∂ T)_V, Y_e. N_B
          -(∂ P/∂ T)_V, Y_e, N_B( ∂s̃/∂ V)_T,Y_e, N_B

and similarly the denominator is given by 

    ∂(ϵ,s̃,Y_e,N_B)/∂(V,T,Y_e,N_B )   =   ∂( ϵ, s̃, Y_e, N_B)/∂( V, s̃, Y_e, N_B)∂( V, s̃, Y_e, N_B)/∂( V, T, Y_e, N_B)
       =    -1/V(P+ϵ)
      ( ∂s̃/∂ T)_V,Y_e, N_B

We have correspondingly,


    ( ∂s̃/∂ T)_V, Y_e, N_B =
      1/n_B( ∂ s/∂ T)_n_e,n_B



    (∂s̃/∂ V)_T, Y_e, N_B
    =-1/V( ( ∂ s/∂ n_B)_n_e, T
    + ( ∂ s/∂ n_e)_n_B, T Y_e -
    s/n_B)



    (∂ P/∂ T)_V,Y_e, N_B= 
    ( ∂ P/∂ T)_n_e, n_B



    (∂ P/∂ V)_T,Y_e, N_B= 
    -1/V( (∂ P/∂ n_B)_n_e,T n_B + (∂ P/∂ n_e)_n_B,T n_e
    )

with 

    ( ∂ P/∂ n_B)_n_e,T
    =f_n_B,n_B n_B +f_n_B,n_e n_e


    ( ∂ P/∂ n_e)_n_B,T
    =f_n_B,n_e n_B +f_n_e,n_e n_e


    ( ∂ P/∂ T)_n_B,n_e
    =f_n_B,T n_B +f_n_e,T n_e + s

Finally, we can express C_s^2 as


    C_s^2=-( f_n_B,n_Bn_B^2+2 f_n_B,n_e n_B n_e +
        f_n_e,n_en_e^2 ) f_T,T 
    	+( f_n_B,T n_B +f_n_e,T n_e + s)^2/-(p+ϵ) f_T,T


Since rest mass should be included here, then 

    ϵ   =    ε+m_n n_n+m_p n_p +m_e^- n_e^-+m_e^+ n_e^+

with ε being the internal energy density.

The f_causal can be obtained by analytical
continuation at the boundary following the prescription of
Ref. <cit.>.
Based on thermodynamics, 

    E_npeγ   =   μ_B N_B+μ_L N_e +T S - P_npeγ V
      
    
      dE_npeγ   =   ( μ_B+μ_L Y_e +
      T s̃) dN_B +μ_L N_B dY_e
      
          +T N_B ds̃_npeγ -P dV


Since 

    P   =   -V (∂ε_npeγ/∂ V)_s̃,Y_e,N_B-ε_npeγ
       =    n_B ( ∂ε_npeγ/∂ n_B)
      _s̃,Y_e,N_B -ε_npeγ

By making use of 

    ( ∂ P/∂ n_B)_s̃,Y_e,N_B= n_B ( ∂^2
        ε_npeγ/∂ n_B^2)
      _s̃,Y_e,N_B

Then we obtain a second-order partial differential equation
by substitution of 

    ( ∂ P/∂ n_B)_s̃,Y_e,N_B   =   β( ∂ϵ_epnγ/∂ n_B)_s̃,Y_e,N_B
       =   β( ∂ε_epnγ/∂ n_B)_s̃,Y_e,N_B
      +β m_n

Thus

    β( ∂ε_npeγ/∂ n_B)_s̃,Y_e,N_B
      +β m_n
      =n_B ( ∂^2  ε_npeγ/∂ n_B^2)
      _s̃,Y_e,N_B

Denote ξ=( ∂ε_npeγ/∂ n_B)_s̃,Y_e,N_B=μ_B+μ_L Y_e +T s̃,
The equation above reduces to a first-order differential equation,

    ( ∂ξ/∂ n_B)_s̃,Y_e,N_B
      - β/n_Bξ 
      =β  m_n/n_B

Take β as 0.9 and denote n_B^*, ε^*,
P^* as the boundary condition when speed of sound become superluminal.
Then we have 

    ε_causal      
       =   - m_n n_B + ε^*+m_n n_B^* +P^*/β+1(n_B/n_B^*)^β+1
          +β(ε^*+m_n n_B^* )-P^*/β+1


    P_causal   =   β/β+1(ε^*+m_n n_B^*+P^*)
      (n_B/n_B^*)^β+1
           - β(ε^*+m_n n_B^* )-P^*/β+1


    f_causal=ϵ_causal-T s^*



    μ_n_causal   =   μ_B
     
       =   
      -m_n+ε^*+m_n n_B^* + P^* /n_B^*(n_B/n_B^*)^β



    μ_p_causal=μ_n_causal+ μ_L^* -μ_e^*





 §.§ Old stuff


Derivative wrt n_n, the equations become

    1     =    2/λ^3(
      z_n ∂μ_n/∂ n_n1/T +
      4 z_n^2 b_n ∂μ_n/∂ n_n1/T +
      2 z_p z_n b_pn∂μ_p/∂ n_n1/T. 
          . + 2 z_p z_n b_pn∂μ_n/∂ n_n1/T)
      
      
      0     =    2/λ^3(
      z_p ∂μ_p/∂ n_n1/T+
      4 z_p^2 b_n ∂μ_p/∂ n_n1/T +
      2 z_p z_n b_pn∂μ_p/∂ n_n1/T +
      . 
          . +2 z_p z_n b_pn∂μ_n/∂ n_n1/T) 
      
       =   2/λ^3(
      z_p ∂μ_p/∂ n_n1/T+
      4 z_p^2 b_n ∂μ_p/∂ n_n1/T +
      2 z_p z_n b_pn∂μ_p/∂ n_n1/T +
      . 
          . +2 z_p z_n b_pn∂μ_n/∂ n_n1/T)
 
Derivative wrt n_p, the equations become


    0     =    2/λ^3(
      z_n ∂μ_n/∂ n_p1/T +
      4 z_n^2 b_n ∂μ_n/∂ n_p1/T +
      2 z_p z_n b_pn∂μ_p/∂ n_p1/T. 
          . + 2 z_p z_n b_pn∂μ_n/∂ n_p1/T)
      
      
      1     =    2/λ^3(
      z_p ∂μ_p/∂ n_p1/T+
      4 z_p^2 b_n ∂μ_p/∂ n_p1/T +
      2 z_p z_n b_pn∂μ_p/∂ n_p1/T +
      . 
          . +2 z_p z_n b_pn∂μ_n/∂ n_p1/T)



apsrev
 

