Class eos_had_skyrme (o2scl)¶

class eos_had_skyrme : public o2scl::eos_had_temp_eden_base ¶

Skyrme hadronic equation of state.

Equation of state of nucleonic matter based on the Skryme interaction from [Skyrme59].

Hamiltonian

The Hamiltonian is defined (using the notation of [Steiner05b] as

\[ {\cal H} = {\cal H}_{k1} + {\cal H}_{k2} + {\cal H}_{k3} + {\cal H}_{p1} + {\cal H}_{p2} + {\cal H}_{p3} + {\cal H}_{g1} + {\cal H}_{g2} \, . \]

The kinetic terms are:

\[ {\cal H}_{k1} = \frac{\tau_n}{2 m_n} + \frac{\tau_p}{2 m_p} \]

\[ {\cal H}_{k2} = n \left(\tau_n + \tau_p \right) \left[ \frac{t_1}{4} \left( 1 + \frac{x_1}{2} \right) + \frac{t_2}{4} \left( 1 + \frac{x_2}{2} \right) \right] \]

\[ {\cal H}_{k3} = \left( \tau_n n_n + \tau_p n_p \right) \left[ \frac{t_2}{4} \left( \frac{1}{2} + x_2 \right) - \frac{t_1}{4} \left( \frac{1}{2} + x_1 \right) \right] \]

where \( \tau_i \) is the Fermi gas energy density of particle \( i \) .

The potential terms are:

\[ {\cal H}_{p1} = \frac{t_0}{2} \left[ \left( 1 + \frac{x_0}{2} \right) n^2 - \left( {\textstyle \frac{1}{2}} + x_0 \right) \left( n_n^2 + n_p^2 \right) \right] \]

\[ {\cal H}_{p2} = \frac{a t_3}{6} \left[ \left( 1 + \frac{x_3}{2} \right) n^{\alpha} n_n n_p + 2^{\alpha-2} \left(1 - x_3\right) \left(n_n^{\alpha+2} + n_p^{\alpha+2}\right) \right] \]

\[ {\cal H}_{p3} = \frac{b t_3}{12} \left[ \left(1 + \frac{x_3}{2} \right) n^{\alpha+2} - \left(\frac{1}{2} + x_3 \right) n^{\alpha} \left( n_n^2+n_p^2 \right) \right] \]

Gradient and Spin-Orbit Terms

The gradient terms are displayed here for completeness even though they are not computed in the code:

\[ {\cal H}_{g1} = \frac{3}{32} \left[ t_1 \left(1 - x_1 \right) - t_2 \left(1 + x_2 \right) \right] \left[ \left( \nabla n_n\right)^2 + \left( \nabla n_p \right)^2 \right] \]

\[ {\cal H}_{g2} = \frac{1}{8} \left[ 3 t_1 \left( 1 + \frac{x_1}{2} \right) - t_2 \left(1 + \frac{x_2}{2} \right) \right] \nabla n_n \nabla n_p \]

The values \(a=0, b=1\) give the standard definition of the Skyrme Hamiltonian [Skyrme59], while \(a=1, b=0\) contains the modifications suggested by [Onsi94].

The spin-orbit term is (following [Steiner05b])

\[ {\cal H}_{J} = -\frac{W_0}{2} \left( n_n \vec{\nabla} \cdot \vec{J}_n + n_p \vec{\nabla} \cdot \vec{J}_p + n \vec{\nabla} \cdot \vec{J} \right) + \frac{t_1}{16} \left(\vec{J}_n^2 + \vec{J}_p^2 - x_1 \vec{J}^2\right) - \frac{t_2}{16} \left(\vec{J}_n^2 + \vec{J}_p^2 + x_2 \vec{J}^2\right) \]

where sometimes the \( \vec{J}^2 \) terms are not included. Alternatively, one can separate the isoscalar and isovector parts in the first term

\[ {\cal H}_{J} = - b_4 n \vec{\nabla} \cdot \vec{J} - b_4^{\prime} n_n \vec{\nabla} \cdot \vec{J}_n - b_4^{\prime} n_p \vec{\nabla} \cdot \vec{J}_p \]

then the earlier Skyrme interactions have \( b_4 = b_4^{\prime} = W_0/2 \). For example, for SLy4, \( b_4 = b_4^{\prime} = W_0/2 = 61.5~\mathrm{MeV} \).

Three quantities are defined in [Steiner05b] for use in computing the properties of matter at saturation

\[ t_3^{\prime} = \left(a + b\right) t_3 \, , \]

\[ C = \frac{3 }{10 m} \left( \frac{3 \pi^2 } {2} \right)^{2/3} \, , \]

and

\[\begin{split} \beta = \frac{M}{2} \left[ \frac{1}{4} \left( 3 t_1 + 5 t_2 \right) + t_2 x_2 \right] \, . \\ \end{split}\]

Units

Quantities which have units containing powers of energy are divided by \(\hbar c\) to ensure all quantities are in units of \( \mathrm{fm} \). The \(x_i\) and \(\alpha\) are unitless, while the original units of the \(t_i\) are:

\(t_0\) - \(\mathrm{MeV}\) \(\mathrm{fm}^3\)
\(t_1\) - \(\mathrm{MeV}\) \(\mathrm{fm}^5\)
\(t_2\) - \(\mathrm{MeV}\) \(\mathrm{fm}^5\)
\(t_3\) - \(\mathrm{MeV}\) \(\mathrm{fm}^{3(1+\alpha)}\)

These are stored internally with units of:

\(t_0\) - \(\mathrm{fm}^2\)
\(t_1\) - \(\mathrm{fm}^4\)
\(t_2\) - \(\mathrm{fm}^4\)
\(t_3\) - \(\mathrm{fm}^{2+3 \alpha}\)

Other Notes

The functions for the usual saturation properties are based partly on [Brack85].

Models are taken from the references: [Bartel79], [Beiner75], [Chabanat95], [Chabanat97], [Danielewicz09], [Dobaczewski94], [Dutta86], [Friedrich86], [Onsi94], [Reinhard95], [Reinhard99], [Tondeur84], and [VanGiai81] .

The variables \( \nu_n \) and \( \nu_p \) contain the expressions \( (-\mu_n+V_n)/T \) and \( (-\mu_p+V_p)/T \) respectively, where \( V \) is the potential part of the single particle energy for particle i (i.e. the derivative of the Hamiltonian w.r.t. density while energy density held constant). Equivalently, \( \nu_n\) is just \( -k_{F_n}^2/ 2 m_n^{*} \).

Skyrme models can be loaded using o2scl_hdf::skyrme_load() . The full list is given in the repository in o2scl/data/o2scl/skdata/model_list.

Todos

Todo

In class eos_had_skyrme:

Convert W0 to b4 and b4p everywhere
Remove use of mnuc in calparfun()?
Update reference list.

Idea for Future:

This EOS typically converges very well. One exception seems to be using calc_temp_p() at very low densities. I have had problems, for example, with mun=5.0, mup=6.5 at T=1.0/197.33.

Note

The finite temperature code does not attempt to include antiparticles and uses o2scl::fermion_nonrel_tl::calc_density(). At finite temperature, pure neutron matter implies a zero proton number density which means the proton chemical potential is \( - \infty \) and thus set to the additive inverse of the return value of the function numeric_limits<double>::infinity() The case of pure proton matter is handled similarly. Negative densities result in calling the error handler.

Subclassed by o2scl::eos_had_sym4_skyrme

Basic Skyrme model parameters

double t0¶

double t1¶

double t2¶

double t3¶

double x0¶

double x1¶

double x2¶

double x3¶

double alpha¶

double a¶

double b¶

Other parameters

double W0¶

Spin-orbit splitting (in \( \mathrm{fm}^{-1} \))

This is unused, but included for possible future use and present in the internally stored models.

double b4¶: Isoscalar spin-orbit term (in \( \mathrm{fm}^{-1} \))

double b4p¶: Isovector spin-orbit term (in \( \mathrm{fm}^{-1} \))

std::string reference¶: Bibliographic reference.

bool parent_method¶

Use eos_had_base methods for saturation properties.

This can be set to true to check the difference between the exact expressions and the numerical values from class eos_had_base.

Particle classes

fermion_nonrel nrf¶: Thermodynamics of non-relativistic fermions.

fermion_deriv_nr nrfd¶: Thermodynamics of non-relativistic fermions with derivatives.

Functions and parameters for calpar()

double fixn0¶

double fixeoa¶

double fixesym¶

double fixcomp¶

double fixmsom¶

int calparfun(size_t nv, const ubvector &x, ubvector &y)¶

int calparfun2(size_t nv, const ubvector &x, ubvector &y)¶

EOS helper functions

void hamiltonian_coeffs(double &ham1, double &ham2, double &ham3, double &ham4, double &ham5, double &ham6)¶: Compute the coefficients of the potential energy which have unique dependence on the densities.

template<class fermion_t> inline void base_thermo(fermion_t &ne, fermion_t &pr, double ltemper, thermo &locth, double term, double term2, double ham1, double ham2, double ham3, double ham4, double ham5, double ham6)¶

Compute the base thermodynamic quantities.

This function computes the energy density, pressure, entropy, and chemical potentials.

template<class fermion_t> inline void second_deriv(fermion_t &ne, fermion_t &pr, double ltemper, thermo &locth, thermo_np_deriv_helm &locthd, double term, double term2, double ham1, double ham2, double ham3, double ham4, double ham5, double ham6)¶: Compute second derivatives of the free energy.

template<class fermion_t> inline void check_input(fermion_t &ne, fermion_t &pr, double T)¶

Check the EOS input

This function checks that

the densities and temperature are finite and positive
the spin denegeracies are correct
the masses are sensible
the values of ‘non_interacting’ are false
the alpha parameter is positive
the temperature is not negative

Basic usage

eos_had_skyrme()¶: Create a blank Skyrme EOS.

inline virtual ~eos_had_skyrme()¶: Destructor.

inline eos_had_skyrme(const eos_had_skyrme &f)¶: Copy constructor.

inline eos_had_skyrme &operator=(const eos_had_skyrme &f)¶: Copy construction with operator=()

virtual int calc_temp_e(fermion &ne, fermion &pr, double temper, thermo &th)¶: Equation of state as a function of densities at finite temperature.

virtual int calc_deriv_temp_e(fermion_deriv &ne, fermion_deriv &pr, double temper, thermo &th, thermo_np_deriv_helm &thd)¶: Equation of state as a function of the densities at finite temperature and including second derivatives.

virtual int calc_e(fermion &ne, fermion &pr, thermo &lt)¶: Equation of state as a function of densities at zero temperature.

virtual int calc_deriv_e(fermion_deriv &ne, fermion_deriv &pr, thermo &th, thermo_np_deriv_helm &thd)¶: Equation of state as a function of densities at zero temperature including second derivatives.

inline virtual const char *type()¶: Return string denoting type (“eos_had_skyrme”)

Saturation properties

These calculate the various saturation properties exactly from the parameters at any density. These routines often assume that the neutron and proton masses are equal.

virtual double feoa_symm(double nb)¶

Calculate binding energy in symmetric matter.

\[ \frac{E}{A} = C n_B^{2/3} \left( 1 + \beta n_B \right) + \frac{3 t_0}{8} n_B + \frac{t_3^{\prime}}{16} n_B^{\alpha+1} \]

virtual double fmsom_symm(double nb)¶

Calculate effective mass in symmetric matter.

\[\begin{split} M^{*}/M = \left(1+ \beta n_B \right)^{-1} \\ \end{split}\]

virtual double fcomp_nuc(double nb)¶

Calculate compressibility in nuclear (isospin-symmetric matter)

\[ K = 10 C n_B^{2/3} + \frac{27}{4} t_0 n_B + 40 C \beta n_B^{5/3} + \frac{9 t_3^{\prime}}{16} \alpha \left( \alpha+1 \right) n_B^{1 + \alpha} + \frac{9 t_3^{\prime}}{8} \left( \alpha+1 \right) n_B^{1 + \alpha} \]

virtual double fesym(double nb, double alpha = 0.0)¶

Calculate symmetry energy.

If pf=0.5, then the exact expression below is used. Otherwise, the method from class eos_had_base is used.

\[ E_{sym} = \frac{5}{9} C n^{2/3} + \frac{10 C m}{3} \left[ \frac{t_2}{6} \left(1 + \frac{5}{4} x_2 \right) - \frac{1}{8} t_1 x_1 \right] n^{5/3} - \frac{t_3^{\prime}}{24} \left({\textstyle \frac{1}{2}} + x_3 \right) n^{1+\alpha} - \frac{t_0}{4} \left( {\textstyle \frac{1}{2}} + x_0 \right) n \]

virtual double fkprime_nuc(double nb)¶

Skewness in nuclear (isospin-symetric) matter.

\[ 2 C n_B^{2/3} \left(9-5/M^{*}/M\right)+ \frac{27 t_3^{\prime}}{16} n^{1+\alpha} \alpha \left(\alpha^2-1\right) \]

Compute and test Landau parameters

int check_landau(double nb, double m)¶

Check the Landau parameters for instabilities.

This returns zero if there are no instabilities.

void landau_nuclear(double n0, double m, double &f0, double &g0, double &f0p, double &g0p, double &f1, double &g1, double &f1p, double &g1p)¶

Calculate the Landau parameters for nuclear matter.

Given n0 and m, this calculates the Landau parameters in nuclear matter as given in [Margueron02].

Todo

This function, eos_had_skyrme::landau_nuclear() needs to be checked.

(Checked once on 11/05/03)

void landau_neutron(double n0, double m, double &f0, double &g0, double &f1, double &g1)¶

Calculate the Landau parameters for neutron matter.

Given n0 and m, this calculates the Landau parameters in nuclear matter as given in [Margueron02].

Todo

This function, eos_had_skyrme::landau_neutron() needs to be checked.

(Checked once on 11/05/03)

Other functions

template<class fermion_t> inline void eff_mass(fermion_t &ne, fermion_t &pr, double &term, double &term2)¶: Evaluate the effective masses for neutrons and protons.

int calpar(double gt0 = -10.0, double gt3 = 70.0, double galpha = 0.2, double gt1 = 2.0, double gt2 = -1.0)¶

Calculate \( t_0,t_1,t_2,t_3 \) and \( \alpha \) from the saturation properties.

In nuclear matter:

\( E_b=E_b(n_0,M^{*},t_0,t_3,\alpha) \) \( P=P(n_0,M^{*},t_0,t_3,\alpha) \) \( K=K(n_0,M^{*},t_3,\alpha) \) (the \( t_0 \) dependence vanishes) \( M^{*}=M^{*}(n_0,t_1,t_2,x_2) \) (the \( x_1 \) dependence cancels), \( E_{sym}=E_{sym}(x_0,x_1,x_2,x_3,t_0,t_1,t_2,t_3,\alpha) \)

To fix the couplings from the saturation properties, we take \( n_0, M^{*}, E_b, K \) as inputs, and we can fix \( t_0,t_3,\alpha \) from the first three relations, then use \( M^{*}, E_b \) to fix \( t_2 \) and \( t_1 \). The separation into two solution steps should make for better convergence. All of the x’s are free parameters and should be set before the function call.

The arguments gt0, gt3, galpha, gt1, and gt2 are used as initial guesses for skyme_eos::t0, eos_had_skyrme::t3, eos_had_skyrme::alpha, eos_had_skyrme::t1, and eos_had_skyrme::t2 respectively.

Todo

In function eos_had_skyrme::calpar():

Does this work for both ‘a’ and ‘b’ non-zero?
Compare to similar formulas in [Margueron02].

void alt_params_set(double Crr00, double Crr10, double Crr0D, double Crr1D, double Crt0, double Crt1, double CrDr0, double CrDr1, double CrnJ0, double CrnJ1, double alpha2)¶

Set using alternate parameterization.

From [Bender03] as in, e.g. [Kortelainen14]

\[\begin{split}\begin{eqnarray*} C^{\rho \rho}_{00} &=& 3 t_0/8 \nonumber \\ C^{\rho \rho}_{10} &=& -t_0/4 \left(\frac{1}{2}+x_0 \right) \nonumber \\ C^{\rho \rho}_{0D} &=& t_3/16 \nonumber \\ C^{\rho \rho}_{1D} &=& -t_3/24 \left(\frac{1}{2}+x_3\right) \nonumber \\ C^{\rho \tau}_{0} &=& 3 t_1/16+t_2/4\left(\frac{5}{4}+x_2\right) \nonumber \\ C^{\rho \tau}_{1} &=& -t_1/8 \left(\frac{1}{2}+x_1\right) + t_2/8 \left(\frac{1}{2}+x_2\right) \nonumber \\ C^{\rho \Delta \rho}_{0} &=& -9/64 t_1+t_2/16 \left(\frac{5}{4}+x_2\right) \nonumber \\ C^{\rho \Delta \rho}_{1} &=& 3/32 t_1 \left(\frac{1}{2}+x_1\right) + t_2/32 \left(\frac{1}{2}+x_2\right) \nonumber \\ C^{\rho \nabla J}_{0} &=& -b_4 -b_4^{\prime}/2 \nonumber \\ C^{\rho \nabla J}_{1} &=& -b_4^{\prime}/2 \end{eqnarray*}\end{split}\]

The parameters should have the following units

Crr00: \( \mathrm{fm}^2 \)
Crr10: \( \mathrm{fm}^2 \)
Crr0D: \( \mathrm{fm}^{3 \alpha+2} \)
Crr1D: \( \mathrm{fm}^{3 \alpha+2} \)
Crt0: \( \mathrm{fm}^4 \)
Crt1: \( \mathrm{fm}^4 \)
CrDr0: \( \mathrm{fm}^4 \)
CrDr1: \( \mathrm{fm}^4 \)
CrnJ0: \( \mathrm{fm}^{-1} \)
CrnJ1: \( \mathrm{fm}^{-1} \)
alpha2: unitless

Todo

In eos_had_skyrme::alt_params_set(): These expressions are not exactly the same as those in [Bender03], so I need to find out why and make this more clear.

void alt_params_get(double &Crr00, double &Crr10, double &Crr0D, double &Crr1D, double &Crt0, double &Crt1, double &CrDr0, double &CrDr1, double &CrnJ0, double &CrnJ1, double &alpha2)¶

Get alternate parameterization.

The parameters will have the following units

Crr00: \( \mathrm{fm}^2 \)
Crr10: \( \mathrm{fm}^2 \)
Crr0D: \( \mathrm{fm}^{3 \alpha+2} \)
Crr1D: \( \mathrm{fm}^{3 \alpha+2} \)
Crt0: \( \mathrm{fm}^4 \)
Crt1: \( \mathrm{fm}^4 \)
CrDr0: \( \mathrm{fm}^4 \)
CrDr1: \( \mathrm{fm}^4 \)
CrnJ0: \( \mathrm{fm}^{-1} \)
CrnJ1: \( \mathrm{fm}^{-1} \)
alpha2: unitless

See eos_had_skyrme::alt_params_set().

void alt_params_saturation(double n0, double EoA, double K, double Ms_star, double a, double L, double Mv_star, double CrDr0, double CrDr1, double CrnJ0, double CrnJ1)¶

Use the specified saturation properties and couplings and the function alt_params_set() to set the Skyrme coefficients.

See [Kortelainen10].

This function uses the relations in Kortelainen et al. (2010), The parameters should have the following units

n0: \( \mathrm{fm}^{-3} \)
EoA: \( \mathrm{fm}^{-1} \)
K: \( \mathrm{fm}^{-1} \)
Ms_star: unitless
a: \( \mathrm{fm}^{-1} \)
L: \( \mathrm{fm}^{-1} \)
Mv_star: unitless
CrDr0: \( \mathrm{fm}^{-3} \)
CrDr1: \( \mathrm{fm}^{-3} \)
CrnJ0: \( \mathrm{fm}^{-3} \)
CrnJ1: \( \mathrm{fm}^{-3} \)

Kortelainen et al. (2010) assumed equal neutron and proton masses, so this function uses \( \hbar^2/(2m) = \hbar^2/(m_n+m_p) \) and the neutron and proton masses in eos_had_base::def_neutron and eos_had_base::def_proton, respectively. To obtain the results in the original paper, set neutron and proton masses to ensure that \( \hbar^2/(2m) = 20.73553~\mathrm{MeV}~\mathrm{fm}^2 \) .