EOS Model Documentation

model

class model

Base class for an EOS parameterization.

Subclassed by bamr::qmc_fixp, bamr::qmc_neut, bamr::qmc_threep, bamr::qmc_twolines, bamr::quark_star, bamr::two_polytropes

Return codes for each point

const int ix_success = 0

const int ix_mr_outside = 2

const int ix_r_outside = 3

const int ix_press_dec = 4

const int ix_nb_problem = 5

const int ix_nb_problem2 = 6

const int ix_crust_unstable = 7

const int ix_mvsr_failed = 8

const int ix_tov_failure = 9

const int ix_small_max = 10

const int ix_tov_conv = 11

const int ix_mvsr_table = 12

const int ix_acausal = 13

const int ix_acausal_mr = 14

const int ix_param_mismatch = 15

const int ix_neg_pressure = 16

const int ix_no_eos_table = 17

const int ix_eos_solve_failed = 18

const int ix_trans_invalid = 19

const int ix_SL_invalid = 20

Grids

o2scl::uniform_grid<double> nb_grid

o2scl::uniform_grid<double> e_grid

o2scl::uniform_grid<double> m_grid

Functions for model parameters fixed during the MCMC run

virtual void setup_params(o2scl::cli &cl)

Setup model parameters.

virtual void remove_params(o2scl::cli &cl)

Remove model parameters.

virtual void copy_params(model &m)

Copy model parameters.

Public Functions

virtual void compute_star(const ubvector &pars, std::ofstream &scr_out, int &success, model_data &dat)

Tabulate EOS and then use in cold_nstar.

model(std::shared_ptr<const settings> s, std::shared_ptr<const ns_data> n)

virtual ~model()

virtual void compute_mr(const ubvector &pars, int &success, std::ofstream &scr_out, model_data &dat)

Compute the mass-radius relation.

virtual int compute_point(const ubvector &pars, std::ofstream &scr_out, double &weight, model_data &dat)

Compute the EOS corresponding to parameters in e and put output in tab_eos.

virtual void initial_point(ubvector &pars)

Specify the initial point.

virtual void get_param_info(std::vector<std::string> &names, std::vector<std::string> &units, ubvector &low, ubvector &high)

Set parameter information.

Public Members

double nb_n1

The fiducial baryon density.

double nb_e1

The fiducial energy density.

size_t n_eos_params

Number of EOS parameters.

bool has_eos

True if the model has an EOS.

bool has_esym

True if the model provides S and L.

double nt_low

Lower limit for baryon density of core-crust transition.

double nt_high

Upper limit for baryon density of core-crust transition.

Protected Functions

virtual void compute_eos(const ubvector &pars, int &success, std::ofstream &scr_out, model_data &dat)

Compute the EOS corresponding to parameters in e and put output in tab_eos.

Protected Attributes

nstar_cold2 cns

TOV solver.

The value of o2scl::nstar_cold::nb_start is set to 0.01 by the constructor

o2scl::tov_solve ts

TOV solver.

double schwarz_km

Schwarzchild radius (set in constructor)

o2scl::prob_dens_gaussian nt_a

Gaussians for core-crust transition.

o2scl::prob_dens_gaussian nt_b

o2scl::prob_dens_gaussian nt_c

o2scl::eos_tov_interp teos

EOS interpolation object for TOV solver.

std::shared_ptr<const settings> set

Settings object.

std::shared_ptr<const ns_data> nsd

Mass-radius data.

two_polytropes

class two_polytropes : public bamr::model

Two polytropes (8 parameters) from Steiner10te

Based on the model from [Steiner10te]. The original limits on the parameters are maintained here. This model is referred to as Model A in [Steiner13tn] and was also used in [Lattimer14ns] (where it was the “Base” model) and in [Lattimer14co] .

The EOS from o2scl::eos_had_schematic is used for the EOS near the saturation density. The first parameter is o2scl::eos_had_base::comp (comp), the second is o2scl::eos_had_base::kprime (kprime), the third is used to fix the sum (esym) of o2scl::eos_had_schematic::a and o2scl::eos_had_schematic::b, and the fourth parameter is o2scl::eos_had_schematic::gamma (gamma). The value of o2scl::eos_had_schematic::a defaults to \( 17.0~\mathrm{MeV}/(\hbar c) \), and can be changed by setting the parameter named kin_sym at run time. This EOS is used up to the transition energy density specified by the fifth parameter (trans1). The first polytrope is used with an index specified by the sixth parameter (index1), up to an energy density specified by the seventh parameter (trans2). Finally, the second polytrope is used with an index specified by the eighth parameter (index2).

For a polytrope \( P = K \varepsilon^{1+1/n} \) beginning at a pressure of \( P_1 \), an energy density of \( \varepsilon_1 \) and a baryon density of \( n_{B,1} \), the baryon density along the polytrope is

\[ n_B = n_{B,1} \left(\frac{\varepsilon}{\varepsilon_1}\right)^{1+n} \left(\frac{\varepsilon_1+P_1}{\varepsilon+P}\right)^{n} \, . \]

Similarly, the chemical potential is

\[ \mu_B = \mu_{B,1} \left(1 + \frac{P_1}{\varepsilon_1}\right)^{1+n} \left(1 + \frac{P}{\varepsilon}\right)^{-(1+n)} \, . \]

The expression for the baryon density can be inverted to determine \( \varepsilon(n_B) \)

\[ \varepsilon(n_B) = \left[ \left(\frac{n_{B,1}} {n_B \varepsilon_1} \right)^{1/n} \left(1+\frac{P_1}{\varepsilon_1}\right)-K\right]^{-n} \, . \]

Sometimes the baryon susceptibility is also useful

\[ \frac{d \mu_B}{d n_B} = \left(1+1/n\right) \left( \frac{P}{\varepsilon}\right) \left( \frac{\mu_B}{n_B}\right) \, . \]

Subclassed by bamr::alt_polytropes, bamr::fixed_pressure, bamr::generic_quarks

Functions for model parameters fixed during the MCMC run

virtual void setup_params(o2scl::cli &cl)

Setup model parameters.

virtual void remove_params(o2scl::cli &cl)

Remove model parameters.

virtual void copy_params(model &m)

Copy model parameters.

Public Functions

two_polytropes(std::shared_ptr<const settings> s, std::shared_ptr<const ns_data> n)

Create a model object.

virtual ~two_polytropes()

virtual void get_param_info(std::vector<std::string> &names, std::vector<std::string> &units, ubvector &low, ubvector &high)

Set parameter information [pure virtual].

virtual void compute_eos(const ubvector &e, int &success, std::ofstream &scr_out, model_data &dat)

Compute the EOS corresponding to parameters in e and put output in tab_eos.

virtual void initial_point(ubvector &e)

Function to compute the initial guess.

Public Members

o2scl::eos_had_schematic se

Low-density EOS.

Protected Attributes

o2scl::cli::parameter_double p_kin_sym

Parameter for kinetic part of symmetry energy.

o2scl::fermion neut

Neutron for se.

o2scl::fermion prot

Proton for se.

alt_polytropes

class alt_polytropes : public bamr::two_polytropes

Alternate polytropes from Steiner13tn (8 parameters)

This class is referred to as Model B in [Steiner13tn]. This model is just as in bamr::two_polytropes, but in terms of the exponents instead of the polytropic indices. The lower limit on exp1 is 1.5, as in [Steiner13tn], but softer EOSs could be allowed by setting this to zero. This would not change the the final results in [Steiner13tn], because the lowest pressure EOSs came from bamr::fixed_pressure anyway.

For a polytrope \( P = K \varepsilon^{\Gamma} \) beginning at a pressure of \( P_1 \), an energy density of \( \varepsilon_1 \) and a baryon density of \( n_{B,1} \), the baryon density along the polytrope is

\[ n_B = n_{B,1} \left(\frac{\varepsilon}{\varepsilon_1} \right)^{\Gamma/(\Gamma-1)} \left(\frac{\varepsilon+P} {\varepsilon_1+P_1}\right)^{1/(1-\Gamma)} \]

Public Functions

alt_polytropes(std::shared_ptr<const settings> s, std::shared_ptr<const ns_data> n)

virtual ~alt_polytropes()

virtual void get_param_info(std::vector<std::string> &names, std::vector<std::string> &units, ubvector &low, ubvector &high)

Set parameter information [pure virtual].

virtual void compute_eos(const ubvector &e, int &success, std::ofstream &scr_out, model_data &dat)

Compute the EOS corresponding to parameters in e and put output in tab_eos.

virtual void initial_point(ubvector &e)

Function to compute the initial guess.

fixed_pressure

class fixed_pressure : public bamr::two_polytropes

Fix pressure on a grid of energy densities from Steiner13tn (8 parameters)

This class is referred to as Model C in [Steiner13tn] and was also used in [Lattimer14ns] (where it was the model labeled “Exo”).

This model is computed as in two_polytropes, but instead of using polytropes at high densities, pressures are linearly interpolated on a fixed grid of energy densities. The schematic EOS (o2scl::eos_had_schematic) is used up to an energy density of \( 1~\mathrm{fm^{-4}} \). The last four parameters are pressures named pres1 through pres4. Then the line segments are defined by the points

\[ P(2~\mathrm{fm}^{-4}) - P(1~\mathrm{fm}^{-4}) = \mathrm{pres1}; \quad P(3~\mathrm{fm}^{-4}) - P(2~\mathrm{fm}^{-4}) = \mathrm{pres2}; \quad P(5~\mathrm{fm}^{-4}) - P(3~\mathrm{fm}^{-4}) = \mathrm{pres3}; \quad P(7~\mathrm{fm}^{-4}) - P(5~\mathrm{fm}^{-4}) = \mathrm{pres4} \]

The final line segment is extrapolated up to \( \varepsilon = 10~\mathrm{fm^{-4}} \)

For a linear EOS, \( P = P_1 + c_s^2 (\varepsilon-\varepsilon_1) \) , beginning at a pressure of \( P_1 \) , an energy density of \( \varepsilon_1 \) and a baryon density of \( n_{B,1} \), the baryon density is

\[ n_B = n_{B,1} \left\{ \frac{\left[\varepsilon+ P_1+c_s^2(\varepsilon-\varepsilon_1)\right]} {\varepsilon_1+P_1} \right\}^{1/(1+c_s^2)} \]

Public Functions

fixed_pressure(std::shared_ptr<const settings> s, std::shared_ptr<const ns_data> n)

virtual ~fixed_pressure()

virtual void get_param_info(std::vector<std::string> &names, std::vector<std::string> &units, ubvector &low, ubvector &high)

Set parameter information [pure virtual].

virtual void compute_eos(const ubvector &e, int &success, std::ofstream &scr_out, model_data &dat)

Compute the EOS corresponding to parameters in e and put output in tab_eos.

virtual void initial_point(ubvector &e)

Function to compute the initial guess.

generic_quarks

class generic_quarks : public bamr::two_polytropes

Generic quark model from Steiner13tn (9 parameters)

This class is referred to as Model D in [Steiner13tn].

This model uses o2scl::eos_had_schematic near saturation, a polytrope (with a uniform prior in the exponent like alt_polytropes) and then a generic quark matter EOS at high densities.

[Alford05hs] parameterizes quark matter with

\[ P = \frac{3 b_4}{4 \pi^2} \mu^4 - \frac{3 b_2}{4 \pi^2} \mu^2 -B \]

where \( \mu \) is the quark chemical potential. QCD corrections can be parameterized by expressing \( b_4 \equiv 1-c \) , and values of \( c \) up to 0.4 (or maybe even larger) are reasonable (see discussion after Eq. 4). Note that, in charge-neutral matter in beta equilibrium, \( \sum_{i=u,d,s,e} n_i \mu_i = \mu_B n_B = \mu n_Q \). where \( \mu_B \) and \( n_B \) are the baryon chemical potential and baryon density and \( n_Q \) is the number density of quarks.

The parameter \( b_2 = m_s^2 - 4 \Delta^2 \) for CFL quark matter, and can thus be positive or negative. A largest possible range might be somewhere between \( (400~\mathrm{MeV})^2 \), which corresponds to the situation where the gap is zero and the strange quarks receive significant contributions from chiral symmetry breaking, to \( (150~\mathrm{MeV})^2-4 (200~\mathrm{MeV})^2 \) which corresponds to a bare strange quark with a large gap. In units of \( \mathrm{fm}^{-1} \) , this corresponds to a range of about \( -3.5 \) to \( 4~\mathrm{fm}^{-2} \) . In Alford et al. (2010), they choose a significantly smaller range, from \( -1 \) to \( 1~\mathrm{fm}^{-2} \).

Simplifying the parameterization to

\[ P = a_4 \mu^4 +a_2 \mu^2 - B \]

gives the following ranges

\[ a_4 = 0.045~\mathrm{to}~0.08 \]

and

\[ a_2 = -0.3~\mathrm{to}~0.3~\mathrm{fm}^{-2} \]

for the “largest possible range” described above or

\[ a_2 = -0.08~\mathrm{to}~0.08~\mathrm{fm}^{-2} \]

for the range used by Alford et al. (2010).

The energy density is

\[ \varepsilon = B + a_2 \mu^2 + 3 a_4 \mu^4 \]

Note that

\[\begin{split}\begin{eqnarray*} \frac{dP}{d \mu} &=& 2 a_2 \mu + 4 a_4 \mu^3 = n_Q \nonumber \\ \frac{d\varepsilon}{d \mu} &=& 2 a_2 \mu + 12 a_4 \mu^3 \end{eqnarray*}\end{split}\]

Public Functions

generic_quarks(std::shared_ptr<const settings> s, std::shared_ptr<const ns_data> n)

virtual ~generic_quarks()

virtual void get_param_info(std::vector<std::string> &names, std::vector<std::string> &units, ubvector &low, ubvector &high)

Set parameter information [pure virtual].

virtual void compute_eos(const ubvector &e, int &success, std::ofstream &scr_out, model_data &dat)

Compute the EOS corresponding to parameters in e and put output in tab_eos.

virtual void initial_point(ubvector &e)

Function to compute the initial guess.

quark_star

class quark_star : public bamr::model

A strange quark star model from Steiner13tn (4 parameters)

This class is referred to as Model E in [Steiner13tn].

Public Functions

virtual void setup_params(o2scl::cli &cl)

Setup new parameters.

virtual void remove_params(o2scl::cli &cl)

Remove model-specific parameters.

quark_star(std::shared_ptr<const settings> s, std::shared_ptr<const ns_data> n)

virtual ~quark_star()

int pressure(size_t nv, const ubvector &x, ubvector &y)

Compute the pressure as a function of the chemical potential.

double pressure2(double mu)

Compute the pressure as a function of the chemical potential.

virtual void get_param_info(std::vector<std::string> &names, std::vector<std::string> &units, ubvector &low, ubvector &high)

Set parameter information [pure virtual].

virtual void compute_eos(const ubvector &e, int &success, std::ofstream &scr_out, model_data &dat)

Compute the EOS corresponding to parameters in e and put output in tab_eos.

virtual void initial_point(ubvector &e)

Function to compute the initial guess.

Public Members

double B

The bag constant.

double c

The paramter controlling non-perturbative corrections to \( \mu^4 \).

double Delta

The gap.

double ms

The strange quark mass.

o2scl::mroot_hybrids gmh

The solver to find the chemical potential for zero pressure.

o2scl::root_brent_gsl grb

An alternative root finder.

qmc_neut

class qmc_neut : public bamr::model

Use QMC computations of neutron matter from Steiner12cn (7 parameters)

[Steiner12cn] used a parameterization for neutron matter which is designed to fit results from quantum Monte Carlo (QMC) simulations in [Gandolfi12mm] .

The parameterization is

\[ E_{\mathrm{neut}} = a \left( \frac{n_B}{n_0} \right)^{\alpha} + b \left( \frac{n_B}{n_0} \right)^{\beta} \]

where \( E_{\mathrm{neut}} \) is the energy per particle in neutron matter, \( n_B \) is the baryon number density, and \( n_0 \equiv 0.16~\mathrm{fm}^{-3} \) is the saturation density. The parameter ranges are

\[\begin{split}\begin{eqnarray*} a &=& 13 \pm 0.3~\mathrm{MeV} \nonumber \\ \alpha &=& 0.50 \pm 0.02 \nonumber \\ b &=& 3 \pm 2~\mathrm{MeV} \nonumber \\ \beta &=& 2.3 \pm 0.2 \, . \end{eqnarray*}\end{split}\]

At high density polytropes are used in a way similar to that in bamr::two_polytropes. The transition between neutron matter and the first polytrope is at a baryon density specified in nb_trans. The remaining 3 parameters are index1, trans1, and index2.

In [Steiner12cn] the polytrope indices are between 0.2 and 2.0. The upper limit on polytropic indices has since been changed from 2.0 to 4.0. The transition between the first and second polytrope at the energy density in trans1 which is between 2.0 and 8.0 \(\mathrm{fm}^{-4}\).

Interpolation objects

ubvector ed_corr

ubvector pres_corr

ubvector pres_err

Public Functions

qmc_neut(std::shared_ptr<const settings> s, std::shared_ptr<const ns_data> n)

virtual ~qmc_neut()

virtual void get_param_info(std::vector<std::string> &names, std::vector<std::string> &units, ubvector &low, ubvector &high)

Set parameter information [pure virtual].

virtual void compute_eos(const ubvector &e, int &success, std::ofstream &scr_out, model_data &dat)

Compute the EOS corresponding to parameters in e and put output in tab_eos.

virtual void initial_point(ubvector &e)

Function to compute the initial guess.

Public Members

double nb0

Saturation density in \( \mathrm{fm}^{-3} \).

double nb_trans

Transition density (default 0.48)

o2scl::interp_vec si

Ratio interpolation object.

o2scl::interp_vec si_err

Ratio error interpolation object.

o2scl::prob_dens_gaussian pdg

Gaussian distribution for proton correction factor.

qmc_fixp

class qmc_fixp : public bamr::model

QMC + line segments model created for Steiner15un (8 parameters)

This model was also used in [Fryer15tf], [Nattila16eo], and [Steiner16ns] .

This EOS model is similar to bamr::qmc_threep, except that the high-density EOS is a set of line-segments, similar to bamr::fixed_pressure. The transition between neutron matter from the QMC parameterization and the first line segment is set to a baryon density of nb_trans . The energy density at this transition density is referred to as ed_trans, and the corresponding pressure is pr_trans. The four high-density parameters pres1 through pres4 are then defined by

\[ P(\mathrm{ed1}) - \mathrm{pr\_trans} = \mathrm{pres1}; \quad P(\mathrm{ed2}) - P(\mathrm{ed1}) = \mathrm{pres2}; \quad P(\mathrm{ed3}) - P(\mathrm{ed2}) = \mathrm{pres3}; \quad P(\mathrm{ed4}) - P(\mathrm{ed3}) = \mathrm{pres4} \]

where the energy density grid is set by the class members ed1, ed2, ed3, and ed4. The lower limits on parameters pres1 through pres4 are all zero. The upper limit on pres1 is \( 0.3~\mathrm{fm}^{-4} \). The upper limits on the remaining pressure parameters are set so that the EOS is not acausal (even though causality is separately double-checked by the code in bamr.cpp anyway).

The limits on the high-density EOS parameters are the same as those in bamr::fixed_pressure.

The energy densities which define the grid

double ed1

double ed2

double ed3

double ed4

Public Functions

qmc_fixp(std::shared_ptr<const settings> s, std::shared_ptr<const ns_data> n)

virtual ~qmc_fixp()

virtual void get_param_info(std::vector<std::string> &names, std::vector<std::string> &units, ubvector &low, ubvector &high)

Set parameter information [pure virtual].

virtual void compute_eos(const ubvector &e, int &success, std::ofstream &scr_out, model_data &dat)

Compute the EOS corresponding to parameters in e and put output in tab_eos.

virtual void initial_point(ubvector &e)

Function to compute the initial guess.

Public Members

double nb0

Saturation density in \( \mathrm{fm}^{-3} \).

double nb_trans

Transition baryon density (default 0.16, different than bamr::qmc_neut)

qmc_threep

class qmc_threep : public bamr::model

QMC + three polytropes created for Steiner15un (9 parameters)

This model was also used in [Fryer15tf], [Nattila16eo], and [Steiner16ns] . For neutron-rich matter near the saturation density, this class uses the QMC parameterization from [Steiner12cn] as in bamr::qmc_neut. The parameter ranges for for \(a\) and \(\alpha\) are expanded and \(b\) and \(\beta\) are recast in terms of \(S\) and \(L\).

\[\begin{split}\begin{eqnarray*} a &=& 12.5~\mathrm{to}~13.5~\mathrm{MeV} \nonumber \\ \alpha &=& 0.47~\mathrm{to}~0.53 \nonumber \\ S &=& 29.5~\mathrm{to}~36.1~\mathrm{MeV}\nonumber \\ L &=& 30~\mathrm{to}~70~\mathrm{MeV} \end{eqnarray*}\end{split}\]

The correlation between \( S \) and \( L \) defined by

\[ L < \left(\frac{9.17}{\mathrm{MeV}}\right) S - 266~\mathrm{MeV} \quad \mathrm{and} \quad L > \left(\frac{14.3}{\mathrm{MeV}}\right) S - 379~\mathrm{MeV} \]

from Lattimer14co is enforced. Alternatively, expressing these constraints in \( (S,L) \) space, are between the line through (29,0) and (35,55) and the line through (26.5,0) and (33.5,100) .

Three polytropes are added at high density similar to bamr::two_polytropes and bamr::qmc_neut based on five parameters index1, trans1, index2, trans2, and index3. The transition between neutron matter and the first polytrope is at a baryon density specified in nb_trans. The transition between the first and second polytrope is specified in trans1, and the transition between the second and third polytrope is specified in trans2. The polytropic indices are allowed to be between 0.2 and 8.0 and the transition densities are allowed to be between 0.75 and 8.0 \( \mathrm{fm}^{-4} \).

Public Functions

qmc_threep(std::shared_ptr<const settings> s, std::shared_ptr<const ns_data> n)

virtual ~qmc_threep()

virtual void get_param_info(std::vector<std::string> &names, std::vector<std::string> &units, ubvector &low, ubvector &high)

Set parameter information [pure virtual].

virtual void compute_eos(const ubvector &e, int &success, std::ofstream &scr_out, model_data &dat)

Compute the EOS corresponding to parameters in e and put output in tab_eos.

virtual void initial_point(ubvector &e)

Function to compute the initial guess.

Public Members

double nb0

Saturation density in \( \mathrm{fm}^{-3} \).

double nb_trans

Transition density (default 0.16, different than bamr::qmc_neut)

qmc_twolines

class qmc_twolines : public bamr::model

QMC plus two line segments with arbitrary energy densities (8 parameters)

Public Functions

qmc_twolines(std::shared_ptr<const settings> s, std::shared_ptr<const ns_data> n)

virtual ~qmc_twolines()

virtual void get_param_info(std::vector<std::string> &names, std::vector<std::string> &units, ubvector &low, ubvector &high)

Set parameter information [pure virtual].

virtual void compute_eos(const ubvector &e, int &success, std::ofstream &scr_out, model_data &dat)

Compute the EOS corresponding to parameters in e and put output in tab_eos.

virtual void initial_point(ubvector &e)

Function to compute the initial guess.

Public Members

double nb0

Saturation density in \( \mathrm{fm}^{-3} \).

double nb_trans

Transition density (default 0.16, different than bamr::qmc_neut)