Class eos_had_rmf (o2scl)¶

class eos_had_rmf : public o2scl::eos_had_temp_pres_base ¶

Relativistic mean field theory EOS.

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.

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}\!\!\!{\backslash}} - g_{\omega} {{\omega}\!\!\!{\backslash}} - \frac{g_{\rho}}{2} {{\vec{\rho}}\!\!\!{\backslash}} \vec{\tau} - M_i + g_{\sigma} \sigma - e q_i A\!\!\!{\backslash} \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

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

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.

Additional couplings are from [Steiner05b],

They are: \( 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

In class eos_had_rmf:

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

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.

Note

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.

Note

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 eos_had_base::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.

Subclassed by o2scl::eos_had_rmf_delta, o2scl::eos_had_rmf_hyp, o2scl::eos_had_sym4_rmf

Other data members

size_t calc_e_steps¶

The number of separate calls to the solver that the calc_e functions take (default 20)

Values larger than about \( 10^4 \) are probably not useful.

bool calc_e_relative¶

If true, solve for relative densities rather than absolute densities (default false)

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

bool zm_mode¶: Modifies method of calculating effective masses (default false)

int verbose¶

Verbosity parameter.

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.

Masses

double mnuc¶

The scale \( M \).

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

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¶

eos_had_rmf()¶

inline virtual ~eos_had_rmf()¶: Destructor.

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

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

Functions dealing with naturalness

double n_baryon¶: Temporary baryon density.

double n_charge¶

Temporary charge density.

Idea for Future:: Should use eos_had_base::proton_frac instead?

inline void check_naturalness(eos_had_rmf &re)¶

Set the coefficients of a eos_had_rmf object to their limits from naturalness.

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{split}\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*}\end{split}\]

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{split}\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*}\end{split}\]

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.

inline void naturalness_limits(double value, eos_had_rmf &re)¶

Provide the maximum values of the couplings assuming a limit on naturalness.

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.

The meson fields

double sigma¶

double omega¶

double rho¶

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.

double ce_temp¶: Temperature storage for calc_temp_e()

int fix_saturation_fun(size_t nv, const ubvector &x, ubvector &y)¶: The function for fix_saturation()

int fix_saturation2_fun(size_t nv, const ubvector &x, ubvector &y, double fix_n0, double fix_eoa, double fix_comp, double fix_esym, double fix_msom)¶: The function for fix_saturation2()

virtual int zero_pressure(size_t nv, const ubvector &ex, ubvector &ey, o2scl::thermo &th)¶: Compute matter at zero pressure (for saturation())

virtual int calc_e_solve_fun(size_t nv, const ubvector &ex, ubvector &ey, o2scl::thermo &th)¶: The function for calc_e()

virtual int calc_temp_e_solve_fun(size_t nv, const ubvector &ex, ubvector &ey, o2scl::thermo &th)¶: 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’.

Used by fix_saturation().

Compute EOS

virtual int calc_e(fermion &ne, fermion &pr, thermo &lth)¶

Equation of state as a function of density.

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.

virtual int calc_p(fermion &ne, fermion &pr, thermo &lth)¶

Equation of state as a function of chemical potential.

Solves for the field equations automatically.

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

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.

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?

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.

Analogous to calc_eq_p() except at finite temperature.

virtual int calc_temp_p(fermion &ne, fermion &pr, double T, thermo &lth)¶

Equation of state as a function of chemical potential.

Solves for the field equations automatically.

virtual 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.

Note that the meson masses and mnuc must be specified before calling this function. The neutron and proton masses must both be equal to mnuc, and the neutron and proton objects must have inc_rest_mass=false.

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.

This function uses the solver sat_mroot .

Todo:

Fix this for zm_mode=true
Ensure solver is more robust

int fix_saturation2(double guess_cs = 4.0, double guess_cw = 3.0, double guess_cr = 3.0, double guess_b = 0.001, double guess_c = -0.001)¶: Calculate cs, cw, cr, b, and c from the saturation properties.

virtual int saturation()¶

Calculate properties of nuclear matter at the saturation density.

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.

double fesym_fields(double sig, double ome, double nb)¶

Calculate symmetry energy assuming the field equations have already been solved.

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

double fcomp_fields(double sig, double ome, double nb)¶

Calculate the compressibility assuming the field equations have already been solved.

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

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.

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.

Fields and field equations

int field_eqs(size_t nv, const ubvector &x, ubvector &y, o2scl::thermo &th)¶

A function for solving the field equations.

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.

int field_eqsT(size_t nv, const ubvector &x, ubvector &y, o2scl::thermo &th)¶

A function for solving the field equations at finite temperature.

inline 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()

inline int get_fields(double &sig, double &ome, double &lrho)¶

Return the most recent values of the meson fields.

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

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

Public Functions

int check_derivs(double &dPds, double &dPdw, double &dPdr, fermion &ne, fermion &pr, double sig, double ome, double lrho)¶