Class tov_solve (o2scl)¶

class tov_solve¶

Class to solve the Tolman-Oppenheimer-Volkov equations.

See the Solution of the Tolman-Oppenheimer-Volkov Equations section of the User’s Guide for the mathematical background.

This class uses adaptive integration to compute the gravitational mass, the radius, the baryonic mass (if the EOS supplies the baryon density), and the gravitational potential (if requested). The equation of state may be changed at any time, by specifying the appropriate eos_tov object

Basic Usage

After specifying the EOS through tov_solve::set_eos(), one calls either tov_solve::mvsr(), tov_solve::fixed(), tov_solve::max() or tov_solve::fixed_pr() to compute the mass versus radius curve, the profile of a star of a given fixed mass, the profile of the maximum mass star, or the profile of a star with a fixed central pressure. These functions all generate results in the form of a table_units object, which can be obtained from tov_solve::get_results().

Screen output:

verbose=0 - Nothing
verbose=1 - Basic information
verbose=2 - For each profile computation, report solution information at every kilometer.
verbose=3 - Report profile information at every 20 grid points and output the final interpolation to zero pressure. A keypress is required after each profile.

Mass versus radius curve

The neutron star mass versus radius curve is constructed by computing profiles of all neutron stars with a range of central pressures. This range is from prbegin to prend, and the ratio between successive central pressures is specified in princ (making a logarithmic central pressure grid).

Profiles of fixed mass

Profiles for a fixed gravitational mass are computed by solving for the correct central pressure. In order to ensure that the solver does not accidentally select a central pressure beyond the maximum mass neutron star, the profile with the maximum mass is computed beforehand automatically (using max()). The intial guess to the solver is always the value of fixed_pr_guess, which defaults to \( 5.2 \times 10^{-5}~\mathrm{Msun}/\mathrm{km}^3 \) . Alternatively, the user can specify the central pressure of the maximum mass star so that it does not have to be computed.

In order to handle multiply-branched mass-radius relations, the value of fixed_pr_guess must be specified in order to ensure that the correct profile is generated. Even if fixed_pr_guess is specified, there is no guarantee that the Newton-Raphson will select the correct profile (though it is very likely if the guess for the central pressure is sufficiently accurate).

The fixed() function does not support computing profiles with central pressures with gravitational masses larger than the maximum value. For this, use fixed_pr() .

Profile for the maximum mass star

This is provided by the function max() . Internally, this uses the minimizer specified by set_min() or the default minimizer, def_min, to minimize the function max_fun() . In order to generate a good initial guess, the max() function begins by looping over central pressures from max_begin to max_end and choosing the best guess from that set of configurations.

Profile for a fixed central pressure

This is provided by the function fixed_pr() and is relatively fast because it does not require any solving or minimizing. This function allows central pressures larger than that of the maximum mass star.

Output tables

The functions tov_solve::fixed(), tov_solve::fixed_pr(), and tov_solve::max() produce output tables which represent the profile of the neutron star of the requested mass. The columns are

gm, the enclosed gravitational mass in \( \mathrm{M}_{\odot} \)
r, the radial coordinate in \( \mathrm{km} \)
gp, the gravitational potential (unitless) when calc_gpot is true
rjw, the value of \( g \) (see definition below; present if ang_vel is true)
omega_rat, the value of \( f \) (see definition below; present if ang_vel is true)
bm, the baryonic mass in \( \mathrm{M}_{\odot} \) (when eos_tov::baryon_column is true).
pr, the pressure in user-specified units
ed, the energy density in user-specified units
nb, the baryon density in user-specified units (if eos_tov::baryon_column is true)
sg, the local surface gravity (in \( \mathrm{g}/\mathrm{cm}^{2} \) )
rs, the local redshift (unitless),
dmdr, the derivative of the enclosed gravitational mass in \( \mathrm{M}_{\odot}/\mathrm{km} \)
dlogpdr, the derivative of the natural logarithm of the pressure
dgpdr, the derivative of the gravitational potential in \( 1/\mathrm{km} \) (if calc_gpot is true)
dbmdr, the derivative of the enclosed baryonic mass (if eos_tov::baryon_column is true).

The function tov_solve::mvsr() produces a different kind of output table corresponding to the mass versus radius curve. Some points on the curve may correspond to unstable configurations.

gm, the total gravitational mass in \( \mathrm{M}_{\odot} \)
r, the radius in \( \mathrm{km} \)
gp, the gravitational potential in the center (unitless) when calc_gpot is true
rjw, the value of \( g \) at the surface (see definition below; present if ang_vel is true)
omega_rat, the value of \( f \) at the surface (see definition below; present if ang_vel is true)
bm, total the baryonic mass in \( \mathrm{M}_{\odot} \) (when eos_tov::baryon_column is true).
pr, the central pressure in user-specified units
ed, the central energy density in user-specified units
nb, the central baryon density in user-specified units (if eos_tov::baryon_column is true)
sg, the surface gravity (in \( \mathrm{g}/\mathrm{cm}^{2} \) )
rs, the redshift at the surface,
dmdr, the derivative of the gravitational mass
dlogpdr, the derivative of the natural logarithm of the pressure
dgpdr, the derivative of the gravitational potential in \( 1/\mathrm{km} \) (if calc_gpot is true)
dbmdr, the derivative of the enclosed baryonic mass (if eos_tov::baryon_column is true).

Unit systems

By default, this class operates with energy density and pressure in units of \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) and baryon density in units of \( \mathrm{fm}^{-3} \).

The function set_units(std::string,std::string,std::string) allows one to use different unit systems for energy density, pressure, and baryon density. The following list of units of energy density and pressure are hard-coded into the library and always work:

"g/cm^3",
"erg/cm^3",
"dyne/cm^2",
"MeV/fm^3",
"1/fm^4", and
"Msun/km^3" (i.e. \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) )

The list of hard-coded units for baryon density are:

"1/m^3",
"1/cm^3", and
"1/fm^3".

Other units choices will work if the conversion is either already added to the global unit conversion object (from o2scl_settings.get_convert_units() ) or the global unit conversion object is able to compute them by a system() call to GNU units (see documentation in convert_units). Note that the choice of what units the tables are produced in is independent of the unit system specified in the associated eos_tov object, i.e. the input EOS and output EOS units need not be the same.

Alternatively, using set_units(double,double,double) allows one to specify the conversion factors directly without using the global unit conversion object.

Accuracy

The present code, as demonstrated in the tests, gives the correct central pressure and energy density of the analytical solution by Buchdahl to within less than 1 part in \( 10^8 \).

Rotation

Rotation is considered if tov_solve::ang_vel is set to true. This adds two more differential equations to solve for each central pressure.

See also the Moment of Inertia in the Slowly-Rotating Approximation section of the User’s Guide.

The differential equation for \( \bar{\omega} \) is independent of the relative scale for \( \bar{\omega} \) and \( j \) . (Note that \( j \) is a metric function not simply related to the angular momentum, \( J \) .) First, one rescales \( \bar{\omega} \) and rewrites everything in terms of \( f\equiv \bar{\omega}/\Omega \) and \( g \equiv r^4 j~df/dr \) . The quantity \( f \) is unitless and \( g \) has units of \( \mathrm{km}^3 \) . Second, pick a central pressure, \( m(r=0)=g(r=0)=0 \), arbitrary values for \( \Phi(r=0) \) and \( f(r=0) \), and integrate

\[\begin{split}\begin{eqnarray*} \frac{d P}{d r} &=& - \frac{G \varepsilon m}{r^2} \left( 1+\frac{P}{\varepsilon}\right) \left( 1+\frac{4 \pi P r^3}{m} \right) \left( 1-\frac{2 G m}{r}\right)^{-1} \nonumber \\ \frac{d m}{d r} &=& 4 \pi r^2 \varepsilon \nonumber \\ \frac{d \Phi}{d r} &=& - \frac{1}{\varepsilon} \frac{ d P}{d r} \left(1+\frac{P}{\varepsilon}\right)^{-1} \nonumber \\ \frac{d g}{dr} &=& -4 r^3 \frac{d j}{dr} f \nonumber \\ \frac{d f}{dr} &=& \frac{g}{r^4 j} \end{eqnarray*}\end{split}\]

Afterwards, shift \( \Phi \) by a constant to ensure the correct value at \( r=R \), and multiply \( g \) by a constant to ensure that \( g=r^4 j (df/dr) \) holds for the new potential \( \Phi \). Then, multiply \( f \) by a constant to ensure that

\[ f(r=R) + \frac{R}{3} \left(\frac{df}{dr}\right)_{r=R} = 1 \]

The functions \( f \) and \( g \) are stored in columns called "omega_rat" and "rjw", respectively. One can compute the baryonic mass by integration or by adding one additional differential equation, bringing the total to six.

The moment of inertia is,

\[ I = \frac{R^4}{6 G} \left.\frac{df}{dr}\right|_{r=R} \]

where the last fraction is stored in domega_rat . For an object named ts of type tov_solve after a call to fixed(), max(), or fixed_pr(), the moment of inertia can be computed with, e.g.

tov_solve ts;
ts.max();
double schwarz_km=o2scl_cgs::schwarzchild_radius/1.0e5;
double I=ts.domega_rat*pow(ts.rad,4.0)/3.0/schwarz_km;

After a call to mvsr(), the values of \( f(r=R) \) and \( g(r=R) \) are stored in columns labeled "omega_rat" and "rjw". The moment of inertia of a 1.4 solar mass neutron star can be computed with, e.g.

tov_solve ts;
ts.mvsr();
std::shared_ptr<table_units<> > tab=ts.get_results();
double schwarz_km=o2scl_cgs::schwarzchild_radius/1.0e5;
double I_14=tab->interp("gm",1.4,"rjw")/3.0/schwarz_km;

Convergence details

By default, if the TOV solver fails to converge, the error handler is called and an exception is thrown. If o2scl::tov_solve::err_nonconv is false, then o2scl::tov_solve::mvsr(), o2scl::tov_solve::fixed(), and o2scl::tov_solve::max(), return an integer which gives some information about why the solver failed to converge.

If err_nonconv is set to false, then the fixed() function temporarily sets the value of both o2scl::mroot::err_nonconv for the current solver and o2scl::jacobian::err_nonconv for the jacobian object, o2scl::mroot_hybrids::def_jac, of def_solver equal to false.

Other details

The ODE solution is stored in a buffer which can be directly accessed using o2scl::tov_solve::get_rkx(), o2scl::tov_solve::get_rky(), and o2scl::tov_solve::get_rkdydx(). In the case that the ODE steps are small enough that there isn’t enough space in the buffers, the ODE solution data is thinned by a factor of two to allow for the remaining ODE steps to be stored. The size of the buffers can be controlled in buffer_size .

If o2scl::tov_solve::reformat_results is true (the default), then the results are placed in a shared pointer to the table_units object which can be accessed using o2scl::tov_solve::get_results(). The maximum size of the output table can be controlled with max_table_size. The output table may be smaller than this, as it cannot be larger than the number of steps stored in the buffer.

Idea for Future:: Convert to o2scl::ode_iv_solve?

Note

If prend is too close to (or smaller than) the central pressure of the maximum mass star, then the mass-radius curve generated by mvsr() might not include the maximum mass star.

Note

The table generated by mvsr() may include unstable configurations (including those with central pressures larger than that of the maximum mass star).

Note

The function o2scl::tov_solve::integ_star() returns gsl_efailed without calling the error handler in the case that the solver can recover gracefully from, for example, a negative pressure.

User EOS

eos_tov *te¶: The EOS.

bool eos_set¶: True if the EOS has been set.

Units for output table

std::string eunits¶: Units for energy density.

std::string punits¶: Units for pressure.

std::string nunits¶: Units for baryon density.

double efactor¶: unit conversion factor for energy density (default 1.0)

double pfactor¶: unit conversion factor for pressure (default 1.0)

double nfactor¶: unit conversion factor for baryon density (default 1.0)

double min_log_pres¶

Smallest allowed pressure for integration (default: -100)

This quantity can’t be much smaller than -100 since we need to compute numbers near \( e^{-\mathrm{min\_log\_pres}} \)

Idea for Future:: Replace this with the proper value from std::limits?

Integration storage

ubvector rkx¶: Radial coordinate (in kilometers)

std::vector<ubvector> rky¶

ODE functions.

If rky is viewed as a matrix, then the first column of each row is the gravitational mass in solar masses, and the second column is the natural logarithm of the pressure in \( \mathrm{M}_{\odot}/km^3 \) . When calc_gpot is true, the next column is the gravitational potential (which is unitless), and when eos_tov::baryon_column is true, the next column is the baryonic mass in \( \mathrm{M}_{\odot} \).

std::vector<ubvector> rkdydx¶: The derivatives of the ODE functions.

std::shared_ptr<table_units<>> out_table¶: The output table.

Numerical methods

mroot<mm_funct, ubvector, jac_funct> *mroot_ptr¶: The solver for fixed gravitational masses.

min_base *min_ptr¶: The minimizer for maximum masses.

astep_base<ubvector, ubvector, ubvector, ode_funct> *as_ptr¶: The adaptive stepper.

size_t buffer_size¶: Size of the ODE solution buffer (default \( 10^5 \))

size_t max_table_size¶: Maximum number of lines in the output table (default 400)

virtual int derivs(double x, size_t nv, const ubvector &y, ubvector &dydx)¶: The ODE step function.

virtual double max_fun(double maxx)¶: The minimizer function to compute the maximum mass.

virtual int integ_star(size_t ndvar, const ubvector &ndx, ubvector &ndy)¶: The solver function to compute the stellar profile.

tov_solve()¶

virtual ~tov_solve()¶

Basic properties

double mass¶: Gravitational mass (in \( \mathrm{M}_{\odot} \))

double rad¶: Radius (in km)

double bmass¶: Baryonic mass (when computed; in \( \mathrm{M}_{\odot} \))

double gpot¶: Gravitational potential (when computed; unitless)

double last_rjw¶: The value of \( r^4 j df / dr \) at the surface.

double last_f¶: The value of \( \bar{\omega} / \Omega \) at the surface.

double domega_rat¶: The value of \( d (\bar{\omega}/\Omega)/dr \) at the surface (when ang_vel is true)

double pcent_max¶

Maximum value for central pressure in \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) (default \( 10^{20} \) )

This variable is set by the mvsr(), max(), fixed() and fixed_pr() functions and used in integ_star() .

Solution parameters

These parameters can be changed at any time.

bool reformat_results¶

If true, then fixed() and max() reformat results into a o2scl::table_units object.

Note that mvsr() always places its results in the output table independent of the value of this variable.

double baryon_mass¶

The mass of one baryon.

The mass of one baryon in kg for the total baryon mass calculation (defaults to the proton mass).

bool ang_vel¶: If true, solve for the angular velocity (default false)

bool gen_rel¶: Apply general relativistic corrections (default true)

bool calc_gpot¶: calculate the gravitational potential (default false)

double step_min¶: smallest allowed radial stepsize in km (default 1.0e-4)

double step_max¶: largest allowed radial stepsize in km (default 0.05)

double step_start¶: initial radial stepsize in km (default 4.0e-3)

int verbose¶: control for output (default 1)

size_t max_integ_steps¶: Maximum number of integration steps (default 100000)

bool err_nonconv¶: If true, call the error handler if the solution does not converge (default true)

double pmax_default¶: Default value of maximum pressure for maximum mass star in \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) (default \( 10^{20} \))

Mass versus radius parameters

double prbegin¶: Beginning pressure in \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) (default 7.0e-7)

double prend¶: Ending pressure in \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) (default 8.0e-3)

double princ¶: Increment factor for pressure (default 1.1)

std::vector<double> pr_list¶

List of pressures at which more information should be recorded.

If pressures (in the user-specified units) are added to this vector, then in mvsr(), the radial location, enclosed gravitational mass, and (if o2scl::eos_tov::baryon_column is true) enclosed baryon mass are stored in the table for each central pressure. The associated columns are named r0, gm0, bm0, r1, gm1, bm1, etc.

Fixed mass parameter

double fixed_pr_guess¶

Guess for central pressure in \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) (default \( 5.2 \times 10^{-5} \))

This guess is used in the functions fixed().

Maximum mass profile parameters

double max_begin¶: Beginning pressure for maximum mass guess in \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) (default 7.0e-5)

double max_end¶: Ending pressure for maximum mass guess in \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) (default 5.0e-3)

double max_inc¶: Increment for pressure for maximum mass guess (unitless; default 1.3)

Default numerical objects

min_brent_gsl def_min¶: The default minimizer.

mroot_hybrids<mm_funct, ubvector, ubmatrix, jac_funct> def_solver¶: The default solver.

astep_gsl<ubvector, ubvector, ubvector, ode_funct> def_stepper¶: The default adaptive stepper.

Convergence information flags

static const int fixed_solver_failed = 128¶

static const int fixed_integ_star_failed = 256¶

static const int fixed_wrong_mass = 512¶

static const int max_minimizer_failed = 1024¶

static const int max_integ_star_failed = 2048¶

static const int mvsr_integ_star_failed = 4096¶

static const int ang_vel_failed = 8192¶

static const int cent_press_large = 16384¶

static const int cent_press_neg = 32768¶

static const int over_max_steps = 65536¶

static const int last_step_large = 131072¶

Basic operation

inline void set_eos(eos_tov &ter)¶: Set the EOS to use.

void set_units(double s_efactor = 1.0, double s_pfactor = 1.0, double s_nbfactor = 1.0)¶: Set output units for the table.

void set_units(std::string eunits = "", std::string punits = "", std::string nunits = "")¶: Set output units for the table.

void get_units(std::string &eunits, std::string &punits, std::string &nunits)¶: Get output units for the table.

virtual int mvsr()¶: Calculate the mass vs. radius curve.

virtual int fixed(double target_mass, double pmax = 1.0e20)¶

Calculate the profile of a star with fixed mass.

This function computes the profile for a star with a fixed mass. If the target mass is negative, it is interpreted as subtracting from the maximum mass configuration. For a 2.15 solar mass neutron star, target_mass=-0.2 corresponds to 1.95 solar mass configuration.

The variable pmax is the maximum allowable central pressure in \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) (This often, but not always, equal to the central pressure of the maximum mass star.) This ensures that the function does not unintentionally select an unstable configuration. If pmax is greater than or equal to the default value (pmax_default), then the maximum mass star will be computed with max() first in order to determine the maximum allowable central pressure.

Note that this function will likely fail when the mass-radius curve has two central pressures with the same gravitational mass.

virtual int fixed_pr(double pcent, double pmax = 1.0e20)¶

Calculate the profile of a star with a specified central pressure.

This function computes the profile of a star with a fixed central pressure. The central pressure, pcent, should be in the unit system specified by the user which defaults to solar masses per cubic kilometer “Msun/km^3” but can be changed with a call to one of the set_units() functions.

The variable pmax is the maximum allowable central pressure in \( \mathrm{M}_{\odot}/\mathrm{km}^3 \), and must be larger than the value of pcent converted to to \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) .

virtual int max()¶

Calculate the profile of the maximum mass star.

Note that this function maximizes the gravitational mass. If the M-R curve has two stable branches, then this function does not necessarily give the configuration with the largest central pressure.

This function may also depend on the accuracy of the initial interval determined by max_begin and max_end.

virtual void make_table()¶

Construct a table from the results.

This function constructs a table_units object from the information in rkx, rky, and rkdydx . Note that the table is always constructed by default so this function need not be called unless tov_solve::reformat_results is false.>

inline std::shared_ptr<table_units<>> get_results()¶: Return the results data table.

inline void set_table(std::shared_ptr<table_units<>> t)¶

Return the results data table.

This function immediately adds four constants to the table, schwarz, Msun, pi and mproton.

Get internally formatted results (in internal unit system)

inline const ubvector &get_rkx() const¶: Get the vector for the radial grid.

inline const std::vector<ubvector> &get_rky() const¶: Get a list of vectors for the ODEs.

inline const std::vector<ubvector> &get_rkdydx() const¶: Get a list of vectors for the corresponding derivatives.

Control numerical methods

inline void set_mroot(mroot<mm_funct, ubvector, jac_funct> &s_mrp)¶: Set solver.

inline void set_min(min_base<> &s_mp)¶: Set minimizer.

inline void set_stepper(astep_base<ubvector, ubvector, ubvector, ode_funct> &sap)¶: Set the adaptive stepper.

Public Types

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

typedef boost::numeric::ublas::matrix<double> ubmatrix¶

typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column¶

typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row¶

Protected Functions

void column_setup(bool mvsr_mode = false)¶

Set up column names and units.

When this function is used by mvsr(), mvsr_mode is set to true.

void make_unique_name(std::string &col, std::vector<std::string> &cnames)¶

Ensure col does not match strings in cnames.

Underscores are added to col until it matches none of the strings in cnames.

Protected Attributes

ode_funct ofm¶: ODE function object.

interp_vec<ubvector> iop¶: Interpolation object for listed radii in mvsr()

bool integ_star_final¶: If true, integ_star() computes all the profile info, otherwise it only computes the gravitational mass.

size_t ix_last¶: The last row index in rky.

double schwarz_km¶: The schwarzchild radius in km.

double tmass¶

Target mass for integ_star()

Has a value of zero, unless set to a non-zero value by fixed().