Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
o2scl::eos_had_rmf Class Reference

Relativistic mean field theory EOS. More...

#include <eos_had_rmf.h>

Inheritance diagram for o2scl::eos_had_rmf:
o2scl::eos_had_temp_pres_base o2scl::eos_had_temp_base o2scl::eos_had_base o2scl::eos_base o2scl::eos_had_rmf_delta o2scl::eos_had_sym4_rmf

Public Member Functions

virtual const char * type ()
 Return string denoting type ("eos_had_rmf")
 
Compute EOS
virtual int calc_e (fermion &ne, fermion &pr, thermo &lth)
 Equation of state as a function of density. More...
 
virtual int calc_p (fermion &ne, fermion &pr, thermo &lth)
 Equation of state as a function of chemical potential. More...
 
virtual int calc_eq_p (fermion &neu, fermion &p, double sig, double ome, double rho, double &f1, double &f2, double &f3, thermo &th)
 Equation of state and meson field equations as a function of chemical potentials. More...
 
virtual int calc_eq_temp_p (fermion &ne, fermion &pr, double temper, double sig, double ome, double rho, double &f1, double &f2, double &f3, thermo &th)
 Equation of state and meson field equations as a function of chemical potentials at finite temperature. More...
 
virtual int calc_temp_p (fermion &ne, fermion &pr, double T, thermo &lth)
 Equation of state as a function of chemical potential. More...
 
int calc_temp_e (fermion &ne, fermion &pr, double T, thermo &lth)
 Equation of state as a function of densities at finite temperature.
 
Saturation properties
int fix_saturation (double guess_cs=4.0, double guess_cw=3.0, double guess_b=0.001, double guess_c=-0.001)
 Calculate cs, cw, cr, b, and c from the saturation properties. More...
 
virtual void saturation ()
 Calculate properties of nuclear matter at the saturation density. More...
 
double fesym_fields (double sig, double ome, double nb)
 Calculate symmetry energy assuming the field equations have already been solved. More...
 
double fcomp_fields (double sig, double ome, double nb)
 Calculate the compressibility assuming the field equations have already been solved. More...
 
void fkprime_fields (double sig, double ome, double nb, double &k, double &kprime)
 Calculate compressibilty and kprime assuming the field equations have already been solved. More...
 
Fields and field equations
int field_eqs (size_t nv, const ubvector &x, ubvector &y)
 A function for solving the field equations. More...
 
int field_eqsT (size_t nv, const ubvector &x, ubvector &y)
 A function for solving the field equations at finite temperature. More...
 
virtual int set_fields (double sig, double ome, double lrho)
 Set a guess for the fields for the next call to calc_e(), calc_p(), or saturation()
 
int get_fields (double &sig, double &ome, double &lrho)
 Return the most recent values of the meson fields. More...
 
Functions dealing with naturalness
void check_naturalness (eos_had_rmf &re)
 Set the coefficients of a eos_had_rmf object to their limits from naturalness. More...
 
void naturalness_limits (double value, eos_had_rmf &re)
 Provide the maximum values of the couplings assuming a limit on naturalness. More...
 
- Public Member Functions inherited from o2scl::eos_had_temp_base
virtual int calc_liqgas_dens_temp_e (fermion &n1, fermion &p1, fermion &n2, fermion &p2, double T, thermo &th1, thermo &th2)
 Compute liquid-gas phase transition densities using eos_had_temp_base::calc_temp_e() . More...
 
virtual int calc_liqgas_temp_e (fermion &n1, fermion &p1, fermion &n2, fermion &p2, double nB, double Ye, double T, thermo &th1, thermo &th2, double &chi)
 Compute the liquid-gas phase transition using eos_had_temp_base::calc_temp_e() . More...
 
virtual int calc_liqgas_beta_temp_e (fermion &n1, fermion &p1, fermion &n2, fermion &p2, double nB, double T, thermo &th1, thermo &th2, double &Ye, double &chi)
 Compute the liquid-gas phase transition in beta-equilibrium using eos_had_temp_base::calc_temp_e() . More...
 
virtual double fesym_T (double nb, double T, double delta=0.0)
 Compute the symmetry energy at finite temperature.
 
virtual double fsyment_T (double nb, double T, double delta=0.0)
 Compute the symmetry entropy at finite temperature.
 
virtual double calc_temp_mun_e (double nn, double np, double T)
 Neutron chemical potential as a function of the densities.
 
virtual double calc_temp_mup_e (double nn, double np, double T)
 Proton chemical potential as a function of the densities.
 
virtual double calc_temp_nn_p (double mun, double mup, double T)
 Neutron density as a function of the chemical potentials.
 
virtual double calc_temp_np_p (double mun, double mup, double T)
 Proton density as a function of the chemical potentials.
 
double calc_fr (double nn, double np, double T)
 Compute the free energy as a function of the temperature and the densities.
 
virtual void f_number_suscept_T (double mun, double mup, double T, double &dPdnn, double &dPdnp, double &dPdpp)
 Compute the number susceptibilities as a function of the chemical potentials, $ \partial^2 P / \partial \mu_i \mu_j $ at a fixed temperature.
 
virtual void f_inv_number_suscept_T (double mun, double mup, double T, double &dednn, double &dednp, double &dedpp)
 Compute the 'inverse' number susceptibilities as a function of the densities, $ \partial^2 \varepsilon / \partial n_i n_j $ at a fixed temperature.
 
void check_en (fermion &n, fermion &p, double T, thermo &th, double &en_deriv, double &en_err)
 Check the entropy by computing the derivative numerically.
 
void check_mu_T (fermion &n, fermion &p, double T, thermo &th, double &mun_deriv, double &mup_deriv, double &mun_err, double &mup_err)
 Check the chemical potentials at finite temperature by computing the derivative numerically.
 
virtual void set_fermion_eval_thermo (fermion_eval_thermo &f)
 Computing finite-temperature integrals. More...
 
- Public Member Functions inherited from o2scl::eos_had_base
virtual double fcomp (double nb, double delta=0.0)
 Calculate the incompressibility in $ \mathrm{fm}^{-1} $ using calc_e() More...
 
virtual double fcomp_err (double nb, double delta, double &unc)
 Compute the incompressibility and its uncertainty. More...
 
virtual double feoa (double nb, double delta=0.0)
 Calculate the energy per baryon in $ \mathrm{fm}^{-1} $ using calc_e() More...
 
virtual double fesym (double nb, double delta=0.0)
 Calculate symmetry energy of matter in $ \mathrm{fm}^{-1} $ using calc_dmu_delta() . More...
 
virtual double fesym_err (double nb, double delta, double &unc)
 Calculate symmetry energy of matter and its uncertainty in $ \mathrm{fm}^{-1} $. More...
 
virtual double fesym_slope (double nb, double delta=0.0)
 The symmetry energy slope parameter in $ \mathrm{fm}^{-1} $. More...
 
virtual double fesym_curve (double nb, double delta=0.0)
 The curvature of the symmetry energy in $ \mathrm{fm}^{-1} $.
 
virtual double fesym_skew (double nb, double delta=0.0)
 The skewness of the symmetry energy in $ \mathrm{fm}^{-1} $.
 
virtual double fesym_diff (double nb)
 Calculate symmetry energy of matter as energy of neutron matter minus the energy of nuclear matter in $ \mathrm{fm}^{-1} $. More...
 
virtual double feta (double nb)
 The strength parameter for quartic terms in the symmetry energy.
 
virtual double feta_prime (double nb)
 The derivative of the strength parameter for quartic terms in the symmetry energy.
 
virtual double fkprime (double nb, double delta=0.0)
 Calculate skewness of nuclear matter in $ \mathrm{fm}^{-1} $ using calc_e() More...
 
virtual double fmsom (double nb, double delta=0.0)
 Calculate reduced neutron effective mass using calc_e() More...
 
virtual double f_effm_neut (double nb, double delta=0.0)
 Neutron (reduced) effective mass.
 
virtual double f_effm_prot (double nb, double delta=0.0)
 Proton (reduced) effective mass.
 
virtual double f_effm_scalar (double nb, double delta=0.0)
 Scalar effective mass. More...
 
virtual double f_effm_vector (double nb, double delta=1.0)
 Vector effective mass. More...
 
virtual double fn0 (double delta, double &leoa)
 Calculate saturation density using calc_e() More...
 
virtual void f_number_suscept (double mun, double mup, double &dPdnn, double &dPdnp, double &dPdpp)
 Compute the number susceptibilities as a function of the chemical potentials, $ \partial^2 P / \partial \mu_i \mu_j $.
 
virtual void f_inv_number_suscept (double mun, double mup, double &dednn, double &dednp, double &dedpp)
 Compute the 'inverse' number susceptibilities as a function of the densities, $ \partial^2 \varepsilon / \partial n_i n_j $.
 
double calc_mun_e (double nn, double np)
 Compute the neutron chemical potential at fixed density. More...
 
double calc_ed (double nn, double np)
 Compute the energy density as a function of the nucleon densities.
 
double calc_pr (double nn, double np)
 Compute the pressure as a function of the nucleon chemical potentials.
 
double calc_mup_e (double nn, double np)
 Compute the proton chemical potential at fixed density. More...
 
double calc_nn_p (double mun, double mup)
 Compute the neutron density at fixed chemical potential. More...
 
double calc_np_p (double mun, double mup)
 Compute the proton density at fixed chemical potential. More...
 
double calc_dmu_delta (double delta, double nb)
 Compute the difference between neutron and proton chemical potentials as a function of the isospin asymmetry. More...
 
double calc_musum_delta (double delta, double nb)
 Compute the sum of the neutron and proton chemical potentials as a function of the isospin asymmetry. More...
 
double calc_pressure_nb (double nb, double delta=0.0)
 Compute the pressure as a function of baryon density at fixed isospin asymmetry. More...
 
double calc_edensity_nb (double nb, double delta=0.0)
 Compute the energy density as a function of baryon density at fixed isospin asymmetry. More...
 
void const_pf_derivs (double nb, double pf, double &dednb_pf, double &dPdnb_pf)
 Compute derivatives at constant proton fraction.
 
double calc_press_over_den2 (double nb, double delta=0.0)
 Calculate pressure / baryon density squared in nuclear matter as a function of baryon density at fixed isospin asymmetry. More...
 
double calc_edensity_delta (double delta, double nb)
 Calculate energy density as a function of the isospin asymmetry at fixed baryon density. More...
 
int nuc_matter_p (size_t nv, const ubvector &x, ubvector &y, double nn0, double np0)
 Solve for the chemical potentials given the densities. More...
 
int nuc_matter_e (size_t nv, const ubvector &x, ubvector &y, double mun0, double mup0)
 Solve for the densities given the chemical potentials. More...
 
virtual void set_mroot (mroot<> &mr)
 Set class mroot object for use in calculating chemical potentials from densities. More...
 
virtual void set_sat_root (root<> &mr)
 Set class mroot object for use calculating saturation density. More...
 
virtual void set_sat_deriv (deriv_base<> &de)
 Set deriv_base object to use to find saturation properties.
 
virtual void set_sat_deriv2 (deriv_base<> &de)
 Set the second deriv_base object to use to find saturation properties. More...
 
virtual void set_n_and_p (fermion &n, fermion &p)
 Set neutron and proton.
 
void gradient_qij (fermion &n, fermion &p, thermo &th, double &qnn, double &qnp, double &qpp, double &dqnndnn, double &dqnndnp, double &dqnpdnn, double &dqnpdnp, double &dqppdnn, double &dqppdnp)
 Calculate coefficients for gradient part of Hamiltonian. More...
 
void check_mu (fermion &n, fermion &p, thermo &th, double &mun_deriv, double &mup_deriv, double &mun_err, double &mup_err)
 Check the chemical potentials by computing the derivatives numerically.
 
void check_den (fermion &n, fermion &p, thermo &th, double &nn_deriv, double &np_deriv, double &nn_err, double &np_err)
 Check the densities by computing the derivatives numerically.
 
- Public Member Functions inherited from o2scl::eos_base
virtual void set_thermo (thermo &th)
 Set class thermo object.
 
virtual const thermoget_thermo ()
 Get class thermo object.
 

Public Attributes

Other data members
size_t calc_e_steps
 The number of separate calls to the solver that the calc_e functions take (default 20) More...
 
bool calc_e_relative
 If true, solve for relative densities rather than absolute densities (default false) More...
 
bool zm_mode
 Modifies method of calculating effective masses (default false)
 
int verbose
 Verbosity parameter. More...
 
bool err_nonconv
 If true, throw exceptions when the function calc_e() does not converge (default true)
 
Masses
double mnuc
 The scale $ M $. More...
 
double ms
 $ \sigma $ mass (in $ \mathrm{fm}^{-1} $ )
 
double mw
 $ \omega $ mass (in $ \mathrm{fm}^{-1} $ )
 
double mr
 $ \rho $ mass (in $ \mathrm{fm}^{-1} $ )
 
Standard couplings (including nonlinear sigma terms)
double cs
 
double cw
 
double cr
 
double b
 
double c
 
Quartic terms for omega and rho.
double zeta
 
double xi
 
Additional isovector couplings
double a1
 
double a2
 
double a3
 
double a4
 
double a5
 
double a6
 
double b1
 
double b2
 
double b3
 
- Public Attributes inherited from o2scl::eos_had_temp_base
fermion_eff def_fet
 Default fermion thermodynamics object.
 
- Public Attributes inherited from o2scl::eos_had_base
double eoa
 Binding energy (without the rest mass) in $ \mathrm{fm}^{-1} $.
 
double comp
 Compression modulus in $ \mathrm{fm}^{-1} $.
 
double esym
 Symmetry energy in $ \mathrm{fm}^{-1} $.
 
double n0
 Saturation density in $ \mathrm{fm}^{-3} $.
 
double msom
 Effective mass (neutron)
 
double kprime
 Skewness in $ \mathrm{fm}^{-1} $.
 
fermion def_neutron
 The defaut neutron. More...
 
fermion def_proton
 The defaut proton. More...
 
deriv_gsl def_deriv
 The default object for derivatives. More...
 
deriv_gsl def_deriv2
 The second default object for derivatives. More...
 
mroot_hybrids def_mroot
 The default solver. More...
 
root_cern def_sat_root
 The default solver for calculating the saturation density. More...
 
- Public Attributes inherited from o2scl::eos_base
thermo def_thermo
 The default thermo object.
 

Protected Member Functions

int fix_saturation_fun (size_t nv, const ubvector &x, ubvector &y)
 The function for fix_saturation()
 
virtual int zero_pressure (size_t nv, const ubvector &ex, ubvector &ey)
 Compute matter at zero pressure (for saturation())
 
virtual int calc_e_solve_fun (size_t nv, const ubvector &ex, ubvector &ey)
 The function for calc_e()
 
virtual int calc_temp_e_solve_fun (size_t nv, const ubvector &ex, ubvector &ey)
 The function for calc_temp_e()
 
int calc_cr (double sig, double ome, double nb)
 Calculate the cr coupling given sig and ome at the density 'nb'. More...
 
- Protected Member Functions inherited from o2scl::eos_had_temp_base
int nuc_matter_temp_e (size_t nv, const ubvector &x, ubvector &y, double nn0, double np0, double T)
 Solve for nuclear matter at finite temperature given density.
 
int nuc_matter_temp_p (size_t nv, const ubvector &x, ubvector &y, double mun0, double mup0, double T)
 Solve for nuclear matter at finite temperature given mu.
 
int liqgas_dens_solve (size_t nv, const ubvector &x, ubvector &y, fermion &n1, fermion &p1, fermion &n2, fermion &p2, double T, thermo &th1, thermo &th2)
 Solve for the liquid gas phase transition as a function of the densities.
 
int liqgas_solve (size_t nv, const ubvector &x, ubvector &y, fermion &n1, fermion &p1, fermion &n2, fermion &p2, double nB0, double Ye0, double T, thermo &th1, thermo &th2)
 Solve for the liquid-gas phase transition at fixed baryon density and electron fraction.
 
int liqgas_beta_solve (size_t nv, const ubvector &x, ubvector &y, fermion &n1, fermion &p1, fermion &n2, fermion &p2, double nB0, double T, thermo &th1, thermo &th2, fermion &e)
 Solve for the liquid-gas phase transition in beta-equilibrium.
 
double calc_entropy_delta (double delta, double nb, double T)
 Compute the entropy.
 
double calc_dmu_delta_T (double delta, double nb, double T)
 Compute the difference between the neutron and proton chemical potentials.
 
- Protected Member Functions inherited from o2scl::eos_had_base
double t1_fun (double barn)
 Compute t1 for gradient_qij().
 
double t2_fun (double barn)
 Compute t2 for gradient_qij().
 

Protected Attributes

double n_baryon
 Temporary baryon density.
 
double n_charge
 Temporary charge density. More...
 
double fe_temp
 Temperature for solving field equations at finite temperature.
 
bool ce_neut_matter
 For calc_e(), if true, then solve for neutron matter.
 
bool ce_prot_matter
 For calc_e(), if true, then solve for proton matter.
 
bool guess_set
 True if a guess for the fields has been given.
 
mroot< mm_funct11, ubvector, jac_funct11 > * sat_mroot
 The solver to compute saturation properties.
 
double ce_temp
 Temperature storage for calc_temp_e()
 
The meson fields
double sigma
 
double omega
 
double rho
 
- Protected Attributes inherited from o2scl::eos_had_temp_base
fermion_eval_thermofet
 Fermion thermodynamics (default is def_fet)
 
- Protected Attributes inherited from o2scl::eos_had_base
mrooteos_mroot
 The EOS solver.
 
rootsat_root
 The solver to compute saturation properties.
 
deriv_basesat_deriv
 The derivative object for saturation properties.
 
deriv_basesat_deriv2
 The second derivative object for saturation properties.
 
fermionneutron
 The neutron object.
 
fermionproton
 The proton object.
 
- Protected Attributes inherited from o2scl::eos_base
thermoeos_thermo
 A pointer to the thermo object.
 

Solver

mroot_hybrids def_sat_mroot
 The default solver for calculating the saturation density. More...
 
virtual int set_sat_mroot (mroot< mm_funct11, ubvector, jac_funct11 > &mrx)
 Set class mroot object for use calculating saturation density.
 

Additional Inherited Members

- Public Types inherited from o2scl::eos_had_base
typedef boost::numeric::ublas::vector< double > ubvector
 

Detailed Description

This class computes the properties of nucleonic matter using a mean-field approximation to a field-theoretical model.

Before sending neutrons and protons to these member functions, the masses should be set to and the degeneracy factor should be set to 2. Some models which can be loaded using o2scl_hdf::rmf_load() expect that the neutron and proton masses are set to the value stored in mnuc.

Note
Since this EOS uses the effective masses and chemical potentials in the o2scl::part class, the values of o2scl::part::non_interacting for neutrons and protons are set to false in many of the functions.
Matter at two different densities can have the same chemical potentials, so the behavior of the function o2scl::eos_had_rmf::calc_temp_p() is ambiguous. This arises because the field equations have more than one solution for a specified chemical potential. Internally, o2scl::eos_had_rmf::calc_temp_p() either uses the initial guess specified by a call to o2scl::eos_had_rmf::set_fields(), or uses hard-coded initial guess values typical for saturation densities. In order to ensure that the user gets the desired solution to the field equations, it may be necessary to specify a sufficiently accurate initial guess. There is no ambiguity in the behavior of o2scl::eos_had_rmf::calc_eq_temp_p(), however.
This class can fail to solve the meson field equations or fail to solve for the nucleon densities. By default the error handler is called when this happens. If err_nonconv is false, then functions which don't converge (which also return int) will return a non-zero value. Note that the solvers (in def_sat_mroot and o2scl::eos_had_base::def_mroot) also has its own data member indicating how to handle nonconvergence o2scl::mroot::err_nonconv which is separate.


Background

The full Lagragian can be written as a sum of several terms

\[ {\cal L} = {\cal L}_{\mathrm{Dirac}} + {\cal L}_{\sigma} + {\cal L}_{\omega} + {\cal L}_{\rho} + {\cal L}_{\mathrm{int}} \, . \]

The part for the nucleon fields is

\[ {\cal L}_{\mathrm{Dirac}} = \bar{\Psi}_i \left[ i {{\partial}\!\!\!{\slash}} - g_{\omega} {{\omega}\!\!\!{\slash}} - \frac{g_{\rho}}{2} {{\vec{\rho}}\!\!\!{\slash}} \vec{\tau} - M_i + g_{\sigma} \sigma - e q_i A\!\!\!{\slash} \right] \Psi_i \]

where $ \Psi $ is the nucleon field and $ \sigma, \omega $ and $ \rho $ are the meson fields. The meson masses are $ m_{\sigma}, m_{\omega} $ and $ m_{\rho} $ and meson-nucleon couplings are $ g_{\sigma}, g_{\omega} $ and $ g_{\rho} $ . The couplings cs, cw, and cr are related to $ g_{\sigma}, g_{\omega} $ and $ g_{\rho} $ by

\[ c_{\sigma} = g_{\sigma}/m_{\sigma} \quad c_{\omega} = g_{\omega}/m_{\omega} \quad \mathrm{and} \quad c_{\rho} = g_{\rho}/m_{\rho} \]

The nucleon masses are in $ M_i $ and stored in part::m and $ q_i $ just represents the charge (1 for protons and 0 for neutrons). The Coulomb field, $ A_{\mu} $, is ignored in this class, but used in o2scl::nucleus_rmf.

The part for the $ \sigma $ field is

\[ {\cal L}_{\sigma} = {\textstyle \frac{1}{2}} \left( \partial_{\mu} \sigma \right)^2 - {\textstyle \frac{1}{2}} m^2_{\sigma} \sigma^2 - \frac{b M}{3} \left( g_{\sigma} \sigma\right)^3 - \frac{c}{4} \left( g_{\sigma} \sigma\right)^4 \, . \]

where $ m_{\sigma} $ is the meson mass, $ b $ and $ c $ are unitless couplings and $ M $ is a dimensionful scale, ususally taken to be 939 MeV (which need not be equal to $ M_i $ above). The coefficients $ b $ and $ c $ are related to the somewhat standard $ \kappa $ and $ \lambda $ by:

\[ \kappa=2 M b \quad \lambda=6 c; \]

The part for the $ \omega $ field is

\[ {\cal L}_{\omega} = - {\textstyle \frac{1}{4}} f_{\mu \nu} f^{\mu \nu} + {\textstyle \frac{1}{2}} m^2_{\omega}\omega^{\mu}\omega_{\mu} + \frac{\zeta}{24} g_{\omega}^4 \left(\omega^\mu \omega_\mu\right)^2 \]

where $ m_{\omega} $ is the meson mass.

The part for the $ \rho $ field is

\[ {\cal L}_{\rho} = - {\textstyle \frac{1}{4}} \vec{B}_{\mu \nu} \cdot \vec{B}^{\mu \nu} + {\textstyle \frac{1}{2}} m^2_{\rho} \vec{\rho}^{~\mu} \cdot \vec{\rho}_{~\mu} + \frac{\xi}{24} g_{\rho}^4 \left(\vec{\rho}^{~\mu}\right) \cdot \vec{\rho}_{~\mu} \]

Finally, additional meson interactions are

\[ {\cal L}_{\mathrm{int}} = g_{\rho}^2 f (\sigma, \omega) \vec{\rho}^{~\mu} \cdot \vec{\rho}_{~\mu} \nonumber \\ \]

The function $ f $ is the coefficient of $ g_r^2 \rho^2 $ $ f(\sigma,\omega) = b_1 \omega^2 + b_2 \omega^4 + b_3 \omega^6 + a_1 \sigma + a_2 \sigma^2 + a_3 \sigma^3 + a_4 \sigma^4 + a_5 \sigma^5 + a_6 \sigma^6 $ where the notation from Horowitz01 is: $ f(\sigma,\omega) = \lambda_4 g_s^2 \sigma^2 + \lambda_v g_w^2 \omega^2 $ This implies $ b_1=\lambda_v g_w^2 $ and $ a_2=\lambda_4 g_s^2 $

The couplings, cs, cw, and cr all have units of $ \mathrm{fm} $, and the couplings b, c, zeta and xi are unitless. The additional couplings from Steiner05b, $ a_i $ have units of $ \mathrm{fm}^{(i-2)} $ and the couplings $ b_j $ have units of $ \mathrm{fm}^{(2j-2)} $ .

When the variable zm_mode is true, the effective mass is fixed using the approach of Zimanyi90 .

The expressions for the energy densities are often simplified in the literature using the field equations. These expressions are not used in this code since they are only applicable in infinite matter where the field equations hold, and are not suitable for use in applications (such as to finite nuclei in o2scl::nucleus_rmf) where the spatial derivatives of the fields are non-zero. Notice that in the proper expressions for the energy density the similarity between terms in the pressure up to a sign. This procedure allows one to verify the thermodynamic identity even if the field equations are not solved and allows the user to add gradient terms to the energy density and pressure.

See also Muller96 .


Field equations

The field equations are:

\[ 0 = m_{\sigma}^2 \sigma - g_{\sigma} \left( n_{s n} + n_{s p} \right) + b M g_{\sigma}^3 \sigma^2 + c g_{\sigma}^4 \sigma^3 - g_{\rho}^2 \rho^2 \frac{\partial f}{\partial \sigma} \]

\[ 0 = m_{\omega}^2 \omega - g_{\omega} \left(n_n+n_p\right) + \frac{\zeta}{6} g_{\omega}^4 \omega^3 + g_{\rho}^2 \rho^2 \frac{\partial f}{\partial \omega} \]

\[ 0 = m_{\rho}^2 \rho + \frac{1}{2} g_{\rho} \left(n_n-n_p\right) + 2 g_{\rho}^2 \rho f + \frac{\xi}{6} g_{\rho}^4 \rho^3 \]


Saturation properties

Defining

\[ U(\sigma)=\frac{1}{2} m_\sigma^2\sigma^2+\frac{b M}{3}(g_\sigma\sigma)^3 +\frac{c}{4}(g_\sigma\sigma)^4\;, \]

the binding energy per particle in symmetric matter at equilibrium is given by

\[ \frac{E}{A} = \frac{1}{n_0} \left[U(\sigma_0)+ \frac{1}{2} m_\omega\omega_0^2+ \frac{\zeta}{8}(g_\omega\omega_0)^4+\frac{2}{\pi^2} \int\limits_0^{k_F} dk k^2\sqrt{k^2+M^{*2}} \right] \]

where the Dirac effective mass is $ M^{*}_i = M_i - g_{\sigma}\sigma_0 $ . The compressibility is given by

\[ K=9\frac{g_\omega^2}{m_\omega^2}n_0+3\frac{k_F^2}{E_F^*} -9n_0\frac{M^{*2}}{E_F^{*2}}\left[\left(\frac{1}{g_\sigma^2} \frac{\partial^2}{\partial\sigma_0^2}+\frac{3}{g_\sigma M^*} \frac{\partial}{\partial\sigma_0}\right) U(\sigma_0)-3\frac{n_0}{E_F^*}\right]^{-1}\;. \]

The symmetry energy of bulk matter is given by

\[ E_{sym} = \frac{k_F^2}{6 E_F^{*}} + \frac{ n } {8 \left(g_{\rho}^2/m_{\rho}^2 + 2 f (\sigma_0, \omega_0) \right)} \, . \]

In the above equations, the subscipt $ 0 $ denotes the mean field values of $ \sigma $ and $ \omega $ . For the case $ f=0 $ , the symmetry energy varies linearly with the density at large densities. The function $ f $ permits variations in the density dependence of the symmetry energy above nuclear matter density.


Todo:
  • The functions fcomp_fields(), fkprime_fields(), and fesym_fields() are not quite correct if the neutron and proton masses are different. For this reason, they are currently unused by saturation().
  • The fix_saturation() and calc_cr() functions use mnuc, and should be modified to allow different neutron and proton masses.
  • Check the formulas in the "Background" section
  • Make sure that this class properly handles particles for which inc_rest_mass is true/false
  • The calc_e() function fails to converge at lower densities. See the testing code which has trouble with NL3 and RAPR
Idea for Future:
  • Finish putting the err_nonconv system into calc_p(), calc_temp_e() and fix_saturation(), etc.
  • It might be nice to remove explicit reference to the meson masses in functions which only compute nuclear matter since they are unnecessary. This might, however, demand redefining some of the couplings.
  • Fix calc_p() to be better at guessing
  • The number of couplings is getting large, maybe new organization is required.
  • Overload eos_had_base::fcomp() with an exact version
  • It would be nice to analytically compute the Jacobian of the field equations for the solver

Definition at line 295 of file eos_had_rmf.h.

Member Function Documentation

◆ calc_cr()

int o2scl::eos_had_rmf::calc_cr ( double  sig,
double  ome,
double  nb 
)
protected

Used by fix_saturation().

◆ calc_e()

virtual int o2scl::eos_had_rmf::calc_e ( fermion ne,
fermion pr,
thermo lth 
)
virtual

Initial guesses for the chemical potentials are taken from the user-given values. Initial guesses for the fields can be set by set_fields(), or default values will be used. After the call to calc_e(), the final values of the fields can be accessed through get_fields().

This is a little more robust than the standard version in the parent eos_had_base.

Idea for Future:
Improve the operation of this function when the proton density is zero.

Reimplemented from o2scl::eos_had_temp_pres_base.

Reimplemented in o2scl::eos_had_rmf_delta.

◆ calc_eq_p()

virtual int o2scl::eos_had_rmf::calc_eq_p ( fermion neu,
fermion p,
double  sig,
double  ome,
double  rho,
double &  f1,
double &  f2,
double &  f3,
thermo th 
)
virtual

This calculates the pressure and energy density as a function of $ \mu_n,\mu_p,\sigma,\omega,rho $ . When the field equations have been solved, f1, f2, and f3 are all zero.

The thermodynamic identity is satisfied even when the field equations are not solved.

Idea for Future:
Probably best to have f1, f2, and f3 scaled in some sensible way, i.e. scaled to the fields?

◆ calc_eq_temp_p()

virtual int o2scl::eos_had_rmf::calc_eq_temp_p ( fermion ne,
fermion pr,
double  temper,
double  sig,
double  ome,
double  rho,
double &  f1,
double &  f2,
double &  f3,
thermo th 
)
virtual

Analogous to calc_eq_p() except at finite temperature.

◆ calc_p()

virtual int o2scl::eos_had_rmf::calc_p ( fermion ne,
fermion pr,
thermo lth 
)
virtual

Solves for the field equations automatically.

Idea for Future:
It may be possible to make the solver for the field equations more robust

Implements o2scl::eos_had_temp_pres_base.

◆ calc_temp_p()

virtual int o2scl::eos_had_rmf::calc_temp_p ( fermion ne,
fermion pr,
double  T,
thermo lth 
)
virtual

Solves for the field equations automatically.

Implements o2scl::eos_had_temp_pres_base.

◆ check_naturalness()

void o2scl::eos_had_rmf::check_naturalness ( eos_had_rmf re)
inline

As given in Muller96 .

The definition of the vector-isovector field and coupling matches what is done here. Compare the Lagrangian above with Eq. 10 from the reference.

The following couplings should all be of the same size:

\[ \frac{1}{2 c_s^2 M^2}, \frac{1}{2 c_v^2 M^2} \frac{1}{8 c_{\rho}^2 M^2},~\mathrm{and}~\frac{ \bar{a}_{ijk} M^{i+2 j+2 k-4}}{2^{2 k}} \]

which are equivalent to

\[ \frac{m_s^2}{2 g_s^2 M^2}, \frac{m_v^2}{2 g_v^2 M^2} \frac{m_{\rho}^2}{8 g_{\rho}^2 M^2},~\mathrm{and}~\frac{ a_{ijk} M^{i+2 j+2 k-4}}{g_s^i g_v^{2 j} g_{\rho}^{2 k} 2^{2 k}} \]

The connection the $ a_{ijk} $ 's and the coefficients that are used here is

\begin{eqnarray*} \frac{b M}{3} g_{\sigma}^3 \sigma^3 &=& a_{300}~\sigma^3 \nonumber \\ \frac{c}{4} g_{\sigma}^4 \sigma^4 &=& a_{400}~\sigma^4 \nonumber \\ \frac{\zeta}{24} g_{\omega}^4 \omega^4 &=& a_{020}~\omega^4 \nonumber \\ \frac{\xi}{24} g_{\rho}^4 \rho^4 &=& a_{002}~\rho^4 \nonumber \\ b_1 g_{\rho}^2 \omega^2 \rho^2 &=& a_{011}~\omega^2 \rho^2 \nonumber \\ b_2 g_{\rho}^2 \omega^4 \rho^2 &=& a_{021}~\omega^4 \rho^2 \nonumber \\ b_3 g_{\rho}^2 \omega^6 \rho^2 &=& a_{031}~\omega^6 \rho^2 \nonumber \\ a_1 g_{\rho}^2 \sigma^1 \rho^2 &=& a_{101}~\sigma^1 \rho^2 \nonumber \\ a_2 g_{\rho}^2 \sigma^2 \rho^2 &=& a_{201}~\sigma^2 \rho^2 \nonumber \\ a_3 g_{\rho}^2 \sigma^3 \rho^2 &=& a_{301}~\sigma^3 \rho^2 \nonumber \\ a_4 g_{\rho}^2 \sigma^4 \rho^2 &=& a_{401}~\sigma^4 \rho^2 \nonumber \\ a_5 g_{\rho}^2 \sigma^5 \rho^2 &=& a_{501}~\sigma^5 \rho^2 \nonumber \\ a_6 g_{\rho}^2 \sigma^6 \rho^2 &=& a_{601}~\sigma^6 \rho^2 \nonumber \end{eqnarray*}

Note that Muller and Serot use the notation

\[ \frac{\bar{\kappa} g_s^3 }{2} = \frac{\kappa}{2} = b M g_s^3 \qquad \mathrm{and} \qquad \frac{\bar{\lambda} g_s^4}{6} = \frac{\lambda}{6} = c g_s^4 \]

which differs slightly from the "standard" notation above.

We need to compare the values of

\begin{eqnarray*} &\frac{m_s^2}{2 g_s^2 M^2}, \frac{m_v^2}{2 g_v^2 M^2} \frac{m_{\rho}^2}{8 g_{\rho}^2 M^2},\frac{b}{3}, \frac{c}{4} & \nonumber \\ &\frac{\zeta}{24}, \frac{\xi}{384}, \frac{b_1}{4 g_{\omega}^2}, \frac{b_2 M^2}{4 g_{\omega}^4}, \frac{b_3 M^4}{4 g_{\omega}^6}, \frac{a_1}{4 g_{\sigma} M},& \nonumber \\ &\frac{a_2}{4 g_{\sigma}^2}, \frac{a_3 M}{4 g_{\sigma}^3}, \frac{a_4 M^2}{4 g_{\sigma}^4}, \frac{a_5 M^3}{4 g_{\sigma}^5},~\mathrm{and}~\frac{a_6 M^4} {4 g_{\sigma}^6}\, .& \end{eqnarray*}

These values are stored in the variables cs, cw, cr, b, c, zeta, xi, b1, etc. in the specified eos_had_rmf object. All of the numbers should be around 0.001 or 0.002.

For the scale $ M $, mnuc is used.

Todo:
I may have ignored some signs in the above, which are unimportant for this application, but it would be good to fix them for posterity.

Definition at line 700 of file eos_had_rmf.h.

◆ fcomp_fields()

double o2scl::eos_had_rmf::fcomp_fields ( double  sig,
double  ome,
double  nb 
)

This may only work at saturation density and may assume equal neutron and proton masses.

◆ fesym_fields()

double o2scl::eos_had_rmf::fesym_fields ( double  sig,
double  ome,
double  nb 
)

This may only work at saturation density and may assume equal neutron and proton masses.

◆ field_eqs()

int o2scl::eos_had_rmf::field_eqs ( size_t  nv,
const ubvector x,
ubvector y 
)

The values x[0], x[1], and x[2] should be set to $ \sigma, \omega $ , and $ \rho $ on input (in $ \mathrm{fm}^{-1} $ ) and on exit, y[0], y[1] and y[2] contain the field equations and are zero when the field equations have been solved.

◆ field_eqsT()

int o2scl::eos_had_rmf::field_eqsT ( size_t  nv,
const ubvector x,
ubvector y 
)

The values x[0], x[1], and x[2] should be set to $ \sigma, \omega $ , and $ \rho $ on input (in $ \mathrm{fm}^{-1} $ ) and on exit, y[0], y[1] and y[2] contain the field equations and are zero when the field equations have been solved.

◆ fix_saturation()

int o2scl::eos_had_rmf::fix_saturation ( double  guess_cs = 4.0,
double  guess_cw = 3.0,
double  guess_b = 0.001,
double  guess_c = -0.001 
)

Note that the meson masses and mnuc must be specified before calling this function.

This function does not give correct results when bool zm_mode is true.

guess_cs, guess_cw, guess_b, and guess_c are initial guesses for cs, cw, b, and c respectively.

Todo:
  • Fix this for zm_mode=true
  • Ensure solver is more robust

◆ fkprime_fields()

void o2scl::eos_had_rmf::fkprime_fields ( double  sig,
double  ome,
double  nb,
double &  k,
double &  kprime 
)

This may only work at saturation density and may assume equal neutron and proton masses.

Todo:
This function, o2scl::eos_had_rmf::fkprime_fields() is currently untested.

◆ get_fields()

int o2scl::eos_had_rmf::get_fields ( double &  sig,
double &  ome,
double &  lrho 
)
inline

This returns the most recent values of the meson fields set by a call to saturation(), calc_e(), or calc_p(fermion &, fermion &, thermo &).

Definition at line 574 of file eos_had_rmf.h.

◆ naturalness_limits()

void o2scl::eos_had_rmf::naturalness_limits ( double  value,
eos_had_rmf re 
)
inline

The limits for the couplings are function of the nucleon and meson masses, except for the limits on b, c, zeta, and xi which are independent of the masses because of the way that these four couplings are defined.

Definition at line 737 of file eos_had_rmf.h.

◆ saturation()

virtual void o2scl::eos_had_rmf::saturation ( )
virtual

This function first constructs an initial guess, increasing the chemical potentials if required to ensure the neutron and proton densities are finite, and then uses eos_had_rmf::sat_mroot to solve the field equations and ensure that the neutron and proton densities are equal and the pressure is zero. The quantities eos_had_base::n0, eos_had_base::eoa, and eos_had_base::msom can be computed directly, and the compressibility, the skewness, and the symmetry energy are computed using the functions fkprime_fields() and fesym_fields(). This function overrides the generic version in eos_had_base.

If verbose is greater than zero, then then this function reports details on the initial iterations to get the initial guess for the solver.

Reimplemented from o2scl::eos_had_base.

Reimplemented in o2scl::eos_had_rmf_delta.

Member Data Documentation

◆ calc_e_relative

bool o2scl::eos_had_rmf::calc_e_relative

Setting this to true makes calc_temp_e() and calc_e() more accurate at low densities.

Definition at line 315 of file eos_had_rmf.h.

◆ calc_e_steps

size_t o2scl::eos_had_rmf::calc_e_steps

Values larger than about $ 10^4 $ are probably not useful.

Definition at line 307 of file eos_had_rmf.h.

◆ def_sat_mroot

mroot_hybrids o2scl::eos_had_rmf::def_sat_mroot

Used by fn0() (which is called by saturation()) to solve saturation_matter_e() (1 variable).

Definition at line 600 of file eos_had_rmf.h.

◆ mnuc

double o2scl::eos_had_rmf::mnuc

This need not be exactly equal to the neutron or proton mass, but provides the scale for the coupling b.

Definition at line 344 of file eos_had_rmf.h.

◆ n_charge

double o2scl::eos_had_rmf::n_charge
protected
Idea for Future:
Should use eos_had_base::proton_frac instead?

Definition at line 779 of file eos_had_rmf.h.

◆ verbose

int o2scl::eos_had_rmf::verbose

If this is greater than zero, then some functions report on their progress.

  • The function saturation() reports progress towards computing the properties of nuclear matter near saturation.
  • The functions calc_e() and calc_temp_e() report progress on solving for matter at a fixed density.

Definition at line 329 of file eos_had_rmf.h.


The documentation for this class was generated from the following file:

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).