Class eos_had_base (o2scl)

O2scl : Class List

class eos_had_base : public o2scl::eos_base

Hadronic equation of state [abstract base].

Denote the number density of neutrons as \( n_n \), the number density of protons as \( n_p \), the total baryon density \( n_B = n_n + n_p \), the asymmetry \( \delta \equiv (n_n-n_p)/n_B \), the nuclear saturation density as \( n_0 \approx 0.16~\mathrm{fm}^{-3} \), and the quantity \( \epsilon \equiv (n_B-n_0)/3n_0 \). (Note that some authors define \( \delta \) as \( n_n - n_p \), which is not the same as the definition above.) Then the energy per baryon of nucleonic matter can be written as an expansion around \( \epsilon =\delta = 0 \)

\[ E(n_B,\delta) = -B + \frac{\tilde{K}}{2!} {\epsilon}^2 + \frac{Q}{3!} {\epsilon}^3 + \delta^2 \left(S + L \epsilon + \frac{K_{\mathrm{sym}}}{2!} {\epsilon}^2 + \frac{Q_{\mathrm{sym}}}{3!} \epsilon^3 \right) + E_4(n_B,\delta) + {\cal O}(\delta^6) \qquad \left(\mathrm{Eq.}~1\right) \]

where \( E_4 \) represents the quartic terms

\[ E_4(n_B,\delta) = \delta^4 \left(S_4 + L_4 \epsilon + \frac{K_4}{2!} {\epsilon}^2 + \frac{Q_4}{3!} \epsilon^3 \right) \qquad \left(\mathrm{Eq.}~2\right) \]

From this, one can compute the energy density of nuclear matter \( \varepsilon(n_B,\delta) = n_B E(n_B,\delta) \), the chemical potentials \( \mu_i \equiv (\partial \varepsilon) / (\partial n_i ) \) and the pressure \( P = -\varepsilon + \mu_n n_n + \mu_p n_p \). This expansion motivates the definition of several separate terms. The binding energy \( B \) of symmetric nuclear matter ( \( \delta = 0 \)) is around \( -16~\mathrm{MeV} \) .

See [Piekarewicz09] for an excellent discussion on which some of this documentation is based.

The compression modulus is usually defined by \( \chi = -1/V (dV/dP) = 1/n (dP/dn)^{-1} \) . In nuclear physics it has become common to use the incompressibility (or bulk) modulus with an extra factor of 9 (motivated by the 3 in the denominator in the definition of \( \epsilon \)), \( K=9/(n \chi) \) and refer to \( K \) simply as the incompressibility. Here, we define the function

\[ K(n_B,\delta) \equiv 9 \left( \frac{\partial P}{\partial n_B} \right) = 9 n_B \left(\frac{\partial^2 \varepsilon} {\partial n_B^2}\right) \]

This quantity is computed by the function fcomp() by computing the first derivative of the pressure, which is more numerically stable than the second derivative of the energy density (since most EOSs compute the pressure exactly). This function is typically evaluated at the point \( (n_B=n_0,\delta=0) \) and is stored in comp by the function saturation(). This quantity is not always the same as \( \tilde{K} \), defined here as

\[ \tilde{K}(n_B,\delta) = 9 n_B^2 \left(\frac{\partial^2 E}{\partial n_B^2}\right) = K(n_B,\delta) - \frac{1}{n_B} 18 P(n_B,\delta) \]

We denote \( K \equiv K(n_B=n_0,\delta=0) \) and similarly for \( \tilde{K} \), the quantity in Eq. 1 above. In nuclear matter at saturation, the pressure is zero and \( K = \tilde{K} \).

See [Chabanat97] for further discussion of the distinction between \(K\) and \(\tilde{K}\).

The symmetry energy, \( S(n_B,\delta), \) can be defined as

\[ S(n_B,\delta) \equiv \frac{1}{2 n_B}\frac{\partial^2 \varepsilon} {\partial \delta^2} \]

and the parameter \( S \) in Eq. 1 is just \( S(n_0,0) \). Using

\[ \left(\frac{\partial \varepsilon}{\partial \delta}\right)_{n_B} = \frac{\partial \varepsilon}{\partial n_n} \left(\frac{\partial n_n}{\partial \delta}\right)_{n_B} + \frac{\partial \varepsilon}{\partial n_p} \left(\frac{\partial n_p}{\partial \delta}\right)_{n_B} = \frac{n_B}{2} \left(\mu_n - \mu_p \right) \]

this can be rewritten

\[ S(n_B,\delta) = \frac{1}{4} \frac{\partial}{\partial \delta} \left(\mu_n - \mu_p\right) \]

where the dependence of the chemical potentials on \( n_B \) and \( \delta \) is not written explicitly. This quantity is computed by function fesym(). Note that many of the functions in this class are written in terms of the proton fraction \( x_p = (1-\delta)/2 \) denoted as 'pf' instead of as functions of \( \delta \). Frequently, \( S(n_B,\delta) \) is evaluated at \( \delta=0 \) to give a univariate function of the baryon density. It is sometimes also evaluated at the point \( (n_B=n_0, \delta=0) \), and this value is denoted by \( S \) above and is typically stored in esym. Alternatively, one can define the symmetry energy by

\[ \tilde{S}(n_B) = E(n_B,\delta=1)-E(n_B,\delta=0) \]

which is computed by function fesym_diff() . The functions \( S(n_B,\delta=0) \) and \( \tilde{S}(n_B) \) are equal when \( {\cal O}(\delta^4) \) terms are zero. In this case, \( \mu_n - \mu_p \) is proportional to \( \delta \) and so

\[ S(n_B) = \tilde{S}(n_B) = \frac{1}{4} \frac{(\mu_n-\mu_p)}{\delta} \, . \]

These functions can also be generalized to finite temperature

\[ S(n_B,\delta,T) = \frac{1}{4} \frac{\partial}{\partial \delta} \left[\mu_n(n_B,\delta,T) - \mu_p(n_B,\delta,T)\right] \, , \]

and

\[ \tilde{S}(n_B,T) = F(n_B,\delta=1,T)-F(n_B,\delta=0,T) \, . \]

The symmetry energy slope parameter \( L \), can be defined by

\[ L(n_B,\delta) \equiv 3 n_B \frac{\partial S(n_B,\delta)} {\partial n_B} = 3 n_B \frac{\partial}{\partial n_B} \left[ \frac{1}{2 n_B} \frac{\partial^2 \varepsilon}{\partial \delta^2} \right] \]

This can be rewritten as

\[ L(n_B,\delta) = \frac{3 n_B}{4} \frac{\partial}{\partial n_B} \frac{\partial}{\partial \delta} \left(\mu_n - \mu_p\right) \]

(where the derivatives can be evaluated in either order) and this is the method used to compute this function in fesym_slope(). Alternatively, using

\[ \left(\frac{\partial \varepsilon}{\partial n_B}\right)_{\delta} = \frac{\partial \varepsilon}{\partial n_n} \left(\frac{\partial n_n}{\partial n_B}\right)_{\delta} + \frac{\partial \varepsilon}{\partial n_p} \left(\frac{\partial n_p}{\partial n_B}\right)_{\delta} = \frac{1}{2} \left(\mu_n + \mu_p \right) \]

\( L \) can be rewritten

\[\begin{split}\begin{eqnarray*} L(n_B,\delta) &=& 3 n_B \left[\frac{-1}{2 n_B^2} \frac{\partial^2 \varepsilon}{\partial \delta^2} + \frac{1}{4 n_B} \frac{\partial^2}{\partial \delta^2} \left(\mu_n + \mu_p\right)\right] \\ &=& \frac{3}{4}\frac{\partial^2}{\partial \delta^2} \left(\mu_n + \mu_p\right) - 3 S(n_B,\delta) \, . \end{eqnarray*}\end{split}\]

The third derivative with respect to the baryon density is sometimes called the skewness. Here, we define

\[ Q(n_B,\delta) = 27 n_B^3 \frac{\partial^3}{\partial n_B^3} \left(\frac{\varepsilon}{n_B}\right) = 27 n_B^3 \frac{\partial^2}{\partial n_B^2} \left(\frac{P}{n_B^2}\right) \]

and this function is computed in fkprime() by computing the second derivative of the pressure.

The second derivative of the symmetry energy with respect to the baryon density is

\[ K_{\mathrm{sym}}(n_B,\delta) = 9 n_B^2 \frac{\partial^2}{\partial n_B^2} S(n_B,\delta) \]

and this function is computed in fesym_curve().

The third derivative of the symmetry energy with respect to the baryon density is

\[ Q_{\mathrm{sym}}(n_B,\delta) = 27 n_B^3 \frac{\partial^3}{\partial n_B^3} S(n_B,\delta) \]

and this function is computed in fesym_skew(). Note that the numerical evaluation of higher derivatives can make eos_had_base::fesym_curve() and eos_had_base::fesym_skew() inaccurate.

Note that assuming terms of order \( \epsilon^3 \) and higher are zero and solving for the baryon density for which \( P=0 \) gives, to order \( \delta^2 \),

\[ n_{B,\mathrm{sat}} = n_0 \left[ 1 - \frac{3 L \delta^2}{K} \right] \]

this implies a new ‘incompressibility’ around the saturation point, i.e.

\[ K(n_B=n_{B,\mathrm{sat}},\delta)= K+\delta^2 \left( K_{\mathrm{sym}}-6 L- \frac{L Q}{K} \right) + {\cal O}\left(\delta^4\right) \]

The quantity in parenthesis is referred to by some authors as \( K_{\tau} \). Note that, because one is evaluating this at \( n_B=n_{B,\mathrm{sat}} \), this is distinct from

\[ \tilde{K}_{\tau} \equiv \frac{1}{2} \frac{\partial^2 K(n_B,\delta)} {\partial \delta^2} \]

which is equal to \( K_{\mathrm{sym}} + 6 L \) to lowest order in \( \delta \) at \( n_B = n_0 \).

The quartic symmetry energy \( S_4(n_B,\delta) \) can be defined as

\[ S_4(n_B,\delta) \equiv \frac{1}{24 n_B}\frac{\partial^4 \varepsilon} {\partial \delta^4} \]

However, fourth derivatives are difficult numerically, and so an alternative quantity is preferable. Instead, one can evaluate the extent to which \( {\cal O}(\delta^4) \) terms are important from

\[ \eta(n_B) \equiv \frac{E(n_B,1)-E(n_B,1/2)} {3 \left[E(n_B,1/2)-E(n_B,0)\right]} \]

(see [Steiner06] and [Piekarewicz09]).

This function can be expressed either in terms of \( \tilde{S} \) or \( S_4 \)

\[ \eta(n_B) = \frac{5 \tilde{S}(n_B) - S(n_B,0)} {\tilde{S}(n_B) + 3 S(n_B,0)} = \frac{5 S_4(n_B,0) + 4 S(n_B,0)} {S_4(n_B,0) + 4 S(n_B,0)} \]

Alternatively, \( S_4 \) can be written

\[ 4 S(n_B) \left[ \frac{1-\eta(n_B)}{\eta(n_B)-5} \right] \, . \]

Evaluating this function at the saturation density gives

\[ \eta(n_0) = \frac{4 S + 5 S_4}{4 S + S_4} \]

(Note that \(S_4\) is referred to as \(Q\) in [Steiner06].)

Sometimes it is useful to separate out the kinetic and potential parts of the energy density when computing \( \eta \), and the class eos_had_sym4_base is useful for this purpose.

The function \( L_4 \) can also be rewritten in \( \eta^{\prime} \) (now suppressing the dependence on \( n_B \)),

\[ \eta^{\prime} = \frac{16 \left( L_4 S - L S_4 \right)} {3 n_B \left(4 S +S_4 \right)^2} \]

then using the expression for \( S_4 \),

\[ \eta^{\prime} = \frac{\left(\eta -5\right)}{48 n_B S } \left[ L_4 \left(\eta -5\right) + 4 L \left(\eta -1\right)\right] \]

Idea for Future:

Replace fmsom() with f_effm_scalar(). This has to wait until f_effm_scalar() has a sensible definition when mn is not equal to mp

Could write a function to compute the “symmetry free energy” or the “symmetry entropy”

Compute the speed of sound or the number susceptibilities?

A lot of the numerical derivatives here might possibly request negative number densities for the nucleons, which may cause exceptions, espescially at very low densities. Since the default EOS objects are GSL derivatives, one can get around these issues by setting the GSL derivative object step size, but this is a temporary solution.

Note

This class and all of its children expect the neutron and proton properties to be specified in powers of \( \mathrm{fm} \) so that masses and chemical potentials are in powers of \( \mathrm{fm}^{-1} \) and energy densities and pressures are in powers of \( \mathrm{fm}^{-4} \) , except where otherwise specified.

Subclassed by o2scl::eos_had_eden_base, o2scl::eos_had_pres_base, o2scl::eos_had_temp_base

Default solvers, derivative classes, and fermions

deriv_gsl def_deriv

The default object for derivatives.

The value of deriv_gsl::h is set to \( 10^{-3} \) in the eos_had_base constructor.

deriv_gsl def_deriv2

The second default object for derivatives.

The value of deriv_gsl::h is set to \( 10^{-3} \) in the eos_had_base constructor.

mroot_hybrids def_mroot

The default solver.

Used by calc_e() to solve nuc_matter_p() (2 variables) and by calc_p() to solve nuc_matter_e() (2 variables).

mroot_hybrids def_sat_mroot

The default solver for calculating the saturation density.

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

fermion def_neutron

The defaut neutron.

By default this has a spin degeneracy of 2 and a mass of o2scl_mks::mass_neutron converted to \( \mathrm{fm}^{-1} \). Also the value of part::non_interacting is set to false.

fermion def_proton

The defaut proton.

By default this has a spin degeneracy of 2 and a mass of o2scl_mks::mass_proton converted to \( \mathrm{fm}^{-1} \). <Also the value of part::non_interacting is set to false.

Consistency checks

mroot_hybrids beta_mroot

Beta-equilibrium solver.

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.

Miscellaneous objects [protected]

fermion *neutron

The neutron object.

fermion *proton

The proton object.

double t1_fun(double barn)

Compute t1 for gradient_qij().

double t2_fun(double barn)

Compute t2 for gradient_qij().

virtual int solve_beta_eq_T0(size_t nv, const ubvector &x, ubvector &y, const double &nB, eos_leptons &elep)

Equation for solving for beta-equilibrium at T=0.

This function is very similar to o2scl::nstar_cold::solve_fun().

Numerical methods [protected]

mroot *eos_mroot

The EOS solver.

mroot *sat_mroot

The solver to compute saturation properties.

deriv_base *sat_deriv

The derivative object for saturation properties.

deriv_base *sat_deriv2

The second derivative object for saturation properties.

Equation of state

virtual int calc_p(fermion &n, fermion &p, thermo &th) = 0

Equation of state as a function of the chemical potentials.

virtual int calc_e(fermion &n, fermion &p, thermo &th) = 0

Equation of state as a function of density.

EOS properties

virtual double fcomp(double nb, double delta = 0.0)

Calculate the incompressibility in \( \mathrm{fm}^{-1} \) using calc_e()

This function computes \( K (n_B,\delta) = 9 n_B \partial^2 \varepsilon /(\partial n_B^2) = 9 \partial P / (\partial n_B) \) . The value of \( K(n_0,0) \), often referred to as the “compressibility”, is stored in comp by saturation() and is about 240 MeV at saturation density.

virtual double fcomp_err(double nb, double delta, double &unc)

Compute the incompressibility and its uncertainty.

This function works like fcomp(), except that it also returns the uncertainty in unc.

virtual double feoa(double nb, double delta = 0.0)

Calculate the energy per baryon in \( \mathrm{fm}^{-1} \) using calc_e()

This function computes the energy per baryon of matter without the nucleon rest masses at the specified baryon density, nb, and isospin asymmetry delta.

virtual double fesym(double nb, double delta = 0.0)

Calculate symmetry energy of matter in \( \mathrm{fm}^{-1} \) using calc_dmu_delta() .

This function computes the symmetry energy,

\[ \left(\frac{1}{2 n_B}\frac{d^2 \varepsilon}{d \delta^2} \right) = \frac{1}{4} \frac{\partial}{\partial \delta} \left(\mu_n - \mu_p \right) \]

at the value of \( n_B \) given in nb and \( \delta \) given in delta. The symmetry energy at \( \delta=0 \) at the saturation density and is stored in esym by saturation().

virtual double fesym_err(double nb, double delta, double &unc)

Calculate symmetry energy of matter and its uncertainty in \( \mathrm{fm}^{-1} \).

This estimates the uncertainty due to the numerical differentiation, assuming that difference betwen the neutron and proton chemical potentials is computed exactly by calc_dmu_delta() .

virtual double fesym_slope(double nb, double delta = 0.0)

The symmetry energy slope parameter in \( \mathrm{fm}^{-1} \).

This returns the value of the “slope parameter” of the symmetry energy as a function of baryon density nb and isospin asymmetry delta. This ranges between about zero and 200 MeV for most equations of state.

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} \).

This function returns the energy per baryon of neutron matter minus the energy per baryon of nuclear matter. This will deviate significantly from the results from fesym() only if the dependence of the symmetry energy on \( \delta \) is not quadratic.

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

The skewness is defined to be \( 27 n^3 d^3 (\varepsilon/n)/(d n^3) = 27 n^3 d^2 (P/n^2) / (d n^2) \) and is denoted ‘kprime’. This definition seems to be ambiguous for densities other than the saturation density and is not quite analogous to the compression modulus.

virtual double fmsom(double nb, double delta = 0.0)

Calculate reduced neutron effective mass using calc_e()

Neutron effective mass (as stored in part::ms) divided by vacuum mass (as stored in part::m) in nuclear matter at saturation density. Note that this simply uses the value of n.ms from calc_e(), so that this effective mass could be either the Landau or Dirac mass depending on the context. Note that this may not be equal to the reduced proton effective mass.

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.

See e.g. [Farine01].

Given the reduced nucleon effective masses, \( m_n^{*} \) and \( m_p^{*} \), the scalar and vector effective masses are defined by

\[ \frac{1}{m^{*}_n} = (1+\delta) \frac{1}{m^{*}_s} - \delta \frac{1}{m^{*}_v} \]

\[ \frac{1}{m^{*}_p} = (1-\delta) \frac{1}{m^{*}_s} + \delta \frac{1}{m^{*}_v} \]

this implies

\[ m_{\mathrm{scalar}}^{*} = \frac{2 m^{*}_n m^{*}_p}{m^{*}_n+m^{*}_p} \]

and

\[ m_{\mathrm{vector}}^{*} = \frac{2 m^{*}_n m^{*}_p \delta}{m^{*}_n - m^{*}_p + \delta(m^{*}_n + m^{*}_p)} \]

virtual double f_effm_vector(double nb, double delta = 1.0)

Vector effective mass.

See documentation for eos_had_base::f_effm_scalar().

Note that the vector effective mass diverges when \( m^{*}_n = m^{*}_p \) and \( \delta = 0 \), but many models have vector effective masses which are independent of \( \delta \). For now, we set \( \delta =1 \) to be the default value, corresponding to neutron matter.

virtual int fn0(double delta, double &nb, double &leoa)

Calculate saturation density using calc_e()

This function finds the baryon density for which the pressure vanishes.

virtual int saturation()

Calculates some of the EOS properties at the saturation density.

This computes the saturation density, and the incompressibility, the symmetry energy, the binding energy, the reduced neutron effective mass at the saturation density, and the skewness in isospin-symmetric matter. The results are stored in n0, comp, esym, eoa, msom, and kprime, respectively.

Idea for Future:

It would be great to provide numerical uncertainties in the saturation properties.

Note

If the saturation density is less than \( 0.08~\mathrm{fm}^{-3} \) or larger than \( 0.24~\mathrm{fm}^{-3} \) then either the error handler is called, or if err_nonconv is false, then a non-zero value is returned.

Nucleonic matter functions

double calc_mun_e(double nn, double np)

Compute the neutron chemical potential at fixed density.

This function uses neutron, proton, eos_base::eos_thermo, and calc_e() .

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.

This function uses neutron, proton, eos_base::eos_thermo, and calc_e() .

double calc_nn_p(double mun, double mup)

Compute the neutron density at fixed chemical potential.

This function uses neutron, proton, eos_base::eos_thermo, and calc_e() .

double calc_np_p(double mun, double mup)

Compute the proton density at fixed chemical potential.

This function uses neutron, proton, eos_base::eos_thermo, and calc_e() .

double calc_dmu_delta(double delta, double nb)

Compute the difference between neutron and proton chemical potentials as a function of the isospin asymmetry.

This function uses neutron, proton, eos_base::eos_thermo, and calc_e() .

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.

This uses neutron, proton, eos_base::eos_thermo, and calc_e() .

double calc_pressure_nb(double nb, double delta = 0.0)

Compute the pressure as a function of baryon density at fixed isospin asymmetry.

Used by fcomp().

int calc_pressure_nb_mroot(size_t nv, const ubvector &x, ubvector &y, double delta)

double calc_edensity_nb(double nb, double delta = 0.0)

Compute the energy density as a function of baryon density at fixed isospin asymmetry.

This uses neutron, proton, eos_base::eos_thermo, and calc_e() .

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.

Used by fkprime().

This uses neutron, proton, eos_base::eos_thermo, and calc_e() .

double calc_edensity_delta(double delta, double nb)

Calculate energy density as a function of the isospin asymmetry at fixed baryon density.

Used by fesym().

This function calls eos_had_base::calc_e() with the internally stored neutron and proton objects.

int nuc_matter_p(size_t nv, const ubvector &x, ubvector &y, double nn0, double np0)

Solve for the chemical potentials given the densities.

The neutron and proton chemical potentials should be stored in x[0] and x[1] and the neutron and proton densities should be stored in pa[0] and pa[1].

Because this function is designed to be used in a solver, it returns exc_efailed without calling the error handler if the densities are not finite.

This function is used by eos_had_pres_base::calc_e().

int nuc_matter_e(size_t nv, const ubvector &x, ubvector &y, double mun0, double mup0)

Solve for the densities given the chemical potentials.

The neutron and proton densities should be stored in x[0] and x[1] and the neutron and proton chemical potentials should be stored in pa[0] and pa[1].

Because this function is designed to be used in a solver, it returns exc_efailed without calling the error handler if the chemical potentials are not finite.

This function is used by eos_had_eden_base::calc_p().

virtual int beta_eq_T0(ubvector &nB_grid, ubvector &guess, eos_leptons &elep, std::shared_ptr<table_units<>> results)

Compute the EOS in beta-equilibrium at zero temperature.

This solves the function solve_beta_eq_T0(). This function is different from nstar_cold because it is a generic interface which works for non-hadronic EOSs.

virtual void f_number_suscept(double mun, double mup, double &dPdnn, double &dPdnp, double &dPdpp)

Compute (numerically) the number susceptibilities as a function of the chemical potentials, \( \partial^2 P / \partial \mu_i \mu_j \).

This function works by taking derivatives of calc_nn_p and calc_np_p using the object specified in set_sat_deriv().

Todo

• This function, eos_had_base::f_number_suscept(), should be overloaded for Skyrme with derivatives

virtual void f_inv_number_suscept(double mun, double mup, double &dednn, double &dednp, double &dedpp)

Compute (numerically) the ‘inverse’ number susceptibilities as a function of the densities, \( \partial^2 \varepsilon / \partial n_i n_j \).

Todo

• This function, eos_had_base::f_inv_number_suscept(), should be overloaded for Skyrme with derivatives

Set auxilliary objects

virtual void set_mroot(mroot<> &mr)

Set class mroot object for use in calculating chemical potentials from densities.

Note

While in principle this allows one to use any mroot object, in practice some of the current EOSs require mroot_hybrids because it automatically avoids regions where the equations are undefined.

virtual void set_sat_mroot(mroot<> &mr)

Set class mroot object for use calculating saturation density.

Note

While in principle this allows one to use any mroot object, in practice some of the current EOSs require mroot_hybrids because it automatically avoids regions where the equations are undefined.

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.

Computing the slope of the symmetry energy at the saturation density requires two derivative objects, because it has to take an isospin derivative and a density derivative. Thus this second deriv_base object is used in the functions fesym_slope(), fesym_curve(), and fesym_skew().

virtual void set_n_and_p(fermion &n, fermion &p)

Set neutron and proton.

Other functions

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 \part of Hamiltonian.

We want the \part of the Hamiltonian in the form

\[ {\cal H}_{\mathrm{grad}} = \frac{1}{2} \sum_{i=n,p} \sum_{j=n,p} Q_{ij} \vec{\nabla} n_i \cdot \vec{\nabla} n_j \]

See e.g. [Pethick95ti].

The expression for the terms from Pethick et al. (1995) is

\[\begin{split}\begin{eqnarray*} {\cal H}_{\mathrm{grad}}&=&-\frac{1}{4} \left(2 P_1 + P_{1;f}-P_{2;f}\right) \nonumber \\ && +\frac{1}{2} \left( Q_1+Q_2 \right) \left(n_n \nabla^2 n_n+n_p \nabla^2 n_p\right) \nonumber \\ && + \frac{1}{4}\left( Q_1-Q_2 \right) \left[\left(\nabla n_n\right)^2+\left(\nabla n_p\right)^2 \right] \nonumber \\ && + \frac{1}{2} \frac{d Q_2}{d n} \left( n_n \nabla n_n + n_p \nabla n_p \right) \cdot \nabla n \end{eqnarray*}\end{split}\]

This can be rewritten

\[\begin{split}\begin{eqnarray*} {\cal H}_{\mathrm{grad}}&=&\frac{1}{2}\left(\nabla n\right)^2 \left[ \frac{3}{2} P_1+n \frac{d P_1}{d n}\right] \nonumber \\ && - \frac{3}{4} \left[ \left( \nabla n_n\right)^2 + \left( \nabla n_p \right)^2 \right] \nonumber \\ && -\frac{1}{2} \left[ \right] \cdot \nabla n \frac{d Q_1}{d n} \nonumber \\ && - \frac{1}{4} \left( \nabla n\right)^2 P_2 - \frac{1}{4} \left[ \left( \nabla n_n\right)^2 + \left( \nabla n_p \right)^2 \right] Q_2 \end{eqnarray*}\end{split}\]

or

\[\begin{split}\begin{eqnarray*} {\cal H}_{\mathrm{grad}}&=&\frac{1}{4} \left( \nabla n\right)^2 \left[3 P_1 + 2 n \frac{d P_1}{d n}-P_2\right] \nonumber \\ && - \frac{1}{4} \left( 3 Q_1+Q_2 \right) \left[ \left( \nabla n_n\right)^2 + \left( \nabla n_p \right)^2 \right] \nonumber \\ && - \frac{1}{2} \frac{d Q_1}{d n} \left[ n_n \nabla n_n + n_p \nabla n_p \right] \cdot \nabla n \end{eqnarray*}\end{split}\]

or

\[\begin{split}\begin{eqnarray*} {\cal H}_{\mathrm{grad}}&=&\frac{1}{4} \left( \nabla n\right)^2 \left[3 P_1 + 2 n \frac{d P_1}{d n}-P_2\right] \nonumber \\ && - \frac{1}{4} \left( 3 Q_1+Q_2 \right) \left[ \left( \nabla n_n\right)^2 + \left( \nabla n_p \right)^2 \right] \nonumber \\ && - \frac{1}{2} \frac{d Q_1}{d n} \left[ n_n \left( \nabla n_n \right)^2 + n_p \left( \nabla n_p \right)^2 + n \nabla n_n \cdot \nabla n_p \right] \end{eqnarray*}\end{split}\]

Generally, for Skyrme-like interactions

\[\begin{split}\begin{eqnarray*} P_i &=& \frac{1}{4} t_i \left(1+\frac{1}{2} x_i \right) \nonumber \\ Q_i &=& \frac{1}{4} t_i \left(\frac{1}{2} + x_i \right) \, . \end{eqnarray*}\end{split}\]

for \( i=1,2 \) .

This function uses the assumption \( x_1=x_2=0 \) to calculate \( t_1 \) and \( t_2 \) from the neutron and proton effective masses assuming the Skyrme form. The values of \( Q_{ij} \) and their derivatives are then computed.

The functions set_n_and_p() and set_thermo() will be called by gradient_qij(), to facilitate the use of the n, p, and th parameters.

Note

This is still somewhat experimental.

inline virtual const char *type()

Return string denoting type (“eos_had_base”)

Public Types

typedef boost::numeric::ublas::vector<double> ubvector

The uBlas vector type.

Public Functions

eos_had_base()

inline virtual ~eos_had_base()

Public Members

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} \).

bool err_nonconv

If true, call the error handler if msolve() or msolve_de() does not converge (default true)