Class nstar_rot (o2scl)¶
-
class nstar_rot¶
Rotating neutron star class based on RNS v1.1d from N. Stergioulas et al.
Several changes have been made to the original code. The code using Numerical Recipes has been removed and replaced with an equivalent based on GSL and . The overall interface has been changed and some code has been updated with C++ equivalents.
Initial guess
The original RNS code suggests that the initial guess is typically a star with a smaller angular momentum.
References
The original RNS v1.1d can be obtained from http://www.gravity.phys.uwm.edu/rns/, and you may find Nick Stergioulas’s web page http://www.astro.auth.gr/~niksterg/, or Sharon Morsink’s page http://fermi.phys.ualberta.ca/~morsink/ useful.
See [Bonazzola73], [Bonazzola94], [Cook92], [Cook94], [Friedman88], [Gourgoulhon94], [Komatsu89], [Laarakkers99], [Nozawa98], [Stergioulas95], and [Stergioulas03].
Todo
In class nstar_rot:
Better documentation is needed everywhere.
Test the resize() function
It appears that KAPPA and KSCALE contant an arbitrary constant, try changing it and see if we get identical results. Try to ensure that the values are consistent between the eos_nstar_rot class and the nstar_rot class.
Variables r_is_gp, p_center, h_center, and others only occur in spherical_star(), integrate(), and make_center(), and can be moved to function parameters or otherwise reorganized.
Directly compare spherical_star() output with tov_solve results
- Idea for Future:
Consider moving the int_z() algorithm to vector_derint.h
Remove the unit-indexed arrays everywhere.
Try moving some of the storage to the heap?
Some of the arrays seem larger than necessary.
The function o2scl::nstar_rot::new_search() is inefficient because it has to handle the boundary conditions separately. This could be improved.
Make the solvers more robust. The ang_vel() and ang_vel_alt() functions appear particularly unstable.
Equation of State
The thermodynamic enthalpy is
\[ H = E + P V \]and the specific enthalpy is often written \( h = H/M \) where \( M \) is some mass scale. In GR this is often chosen to be a generic “baryon mass”, \( m_B \) so that\[ h = \frac{\varepsilon + P}{\rho} \]where \( \varepsilon \equiv E/V \) and \( \rho \equiv m_B/V \). The pseudo-enthalpy is typically defined by\[ d\hat{h} = dP/(P+\varepsilon) \](sometimes denoted as \( H \) or \( h \) and calledh
in the code, but referred to here as \( \hat{h} \) to avoid confusion with the enthalpy and specific enthalpy defined above). Additionally, this quantity is sometimes defined with an additional factor of \( c^2 \). At \( T=0 \) and presuming only one conserved charge, one can use the Gibbs-Duhem relation, \( dP = (\varepsilon + P)/\mu d \mu \) to write\[ d \hat{h} = \frac{d \mu}{\mu} \]\[ \hat{h} = \ln \mu + C \]The constant \( C \) is often chosen to be \( C \equiv - \ln \mu(P=0) \), so that\[ \hat{h} = \ln \frac{\mu}{\mu(P=0)} \]Typically, \( \mu(P=0) \) is around 931 MeV. The eos_nstar_rot_interp class takes an EOS tabulated in powers of \( \mathrm{fm}^{-1} \) and recasts it into a form which can be used by nstar_rot.Coordinate system
The space is mapped to coordinates \( (s,\mu) \) where
\[ r = r_{\mathrm{eq}} \left( \frac{s}{1-s} \right) \]and \( r_{\mathrm{eq}} \) is the equatorial radius, thus\[ s = \frac{r}{r+r_{\mathrm{eq}}} \](This is Eq. 48 in Stergioulas’ thesis)Draft documentation
For spherical stars, the isotropic radius \( r_{\mathrm{is}} \) is defined by
\[ \frac{d r}{d r_{\mathrm{is}}} = \left(\frac{r}{r_{\mathrm{is}}}\right) \left(1 - 2 \frac{m}{r}\right)^{1/2} \]Quadrupole moments
Quadrupole moments computed using the method in [Laarakkers99].
Axisymmetric Instability
[Friedman88] shows that a secular axisymmetric instability sets in when the mass becomes maximum along a sequence of constant angular momentum. Equivalently, [Cook92] shows that the instability occurs when the angular momentum becomes minimum along a sequence of constant rest mass.
A GR virial theorem for a stationary and axisymmetric system was found in [Bonazzola73]. A more general two-dimensional virial identity was found in [Bonazzola94]. The three-dimensional virial identity found in [Gourgoulhon94] is a generalization of the Newtonial virial theorem.
Using the stationary and axisymmetric metric ( \( G = c = 1 \) )
\[ ds^2 = - e^{\gamma+\rho} dt^2 + e^{2 \alpha} \left( dr^2 + r^2 d\theta^2 \right) + e^{\gamma-\rho} r^2 \sin^2 \theta ( d \phi - \omega dt) ^2 \]one solves for the four metric functions \( \rho(r,\theta) \), \( \gamma(r,\theta) \), \( \alpha(r,\theta) \) and \( \omega(r,\theta) \) .It is assumed that matter is a perfect fluid, and the stress-energy tensor is
\[ T^{\mu \nu} = \left( \rho_0 + \rho_i + P \right) u^{\mu} u^{\nu} + P g^{\mu \nu} \]Einstein’s field equations imply four field equations for a specified rotation law,
\[ u^{t} u_{\phi} = F(\Omega) \]for some function \( F(\omega) \) .Using Eq. (27) in [Cook92], one can write
\[ \rho(s,\mu) = - e^{-\gamma/2} \sum_{n=0}^{\infty} P_{2n}(\mu) \int_0^{1}~ds^{\prime} \int_0^1~d \mu f_{\rho}(n,s,s^{\prime}) P_{2n}{\mu^{\prime}} \tilde{S}(s^{\prime},\mu^{\prime}) \]where the function \( f_{\rho} \) is defined by\[ f_{\rho} \equiv \Theta(s^{\prime}-s) \left(\frac{1-s}{s}\right)^{2 n+1} \left[\frac{s^{\prime 2n}}{(1-s^{\prime})^{2n+2}}\right] + \Theta(s^{\prime}-s) \left(\frac{1-s}{s}\right)^{2 n+1} \left[\frac{s^{\prime 2n}}{(1-s^{\prime})^{2n+2}}\right] \]This function is stored in f_rho . Similar definitions are made for f_gamma and f_omega .The Keplerial orbit at the equator is
\[ \Omega_K = \frac{\omega^{\prime}}{2 \psi^{\prime}} ... \](eq. 31 in [Stergioulas03])
Note
This class is still experimental.
EOS member variables
-
double eq_radius_tol_rel¶
Relative accuracy for the equatorial radius, \( r_e \) (default \( 10^{-5} \))
Used in iterate() .
-
double alt_tol_rel¶
Accuracy for equatorial radius using alternate solvers (default \( 10^{-9} \))
-
int verbose¶
Verbosity parameter.
-
bool eos_set¶
If true, then an EOS has been set.
-
bool scaled_polytrope¶
If true, then use a polytrope and rescale.
-
eos_nstar_rot *eosp¶
Pointer to the user-specified EOS.
-
nstar_rot()¶
Output
-
double e_center¶
Central energy density (in units of \( 10^{15} \mathrm{g}/\mathrm{cm}^3 \))
-
double r_ratio¶
Ratio of polar to equatorial radius.
-
double r_e¶
Coordinate equatorial radius.
Quantities computed by nstar_rot::comp() (in order)
-
double r_p¶
Radius at pole.
-
double s_p¶
The value of the s-coordinate at the pole.
-
double s_e¶
The value of the s-coordinate at the equator.
-
double velocity_equator¶
The velocity at the equator.
-
double R_e¶
Circumferential radius in cm (i.e. the radius defined such that \( 2 \pi R_e \) is the proper circumference)
-
double Mass_p¶
Proper mass (in \( \mathrm{g} \) )
-
double Mass¶
Gravitational mass (in \( \mathrm{g} \) )
-
double Mass_0¶
Baryonic mass (in \( \mathrm{g} \) )
-
double J¶
Angular momentum (in \( \mathrm{g}~\mathrm{cm}^2 / \mathrm{s} \) )
-
double Omega¶
Angular velocity (in \( \mathrm{radians}/\mathrm{s} \) )
-
double T¶
Total rotational kinetic energy.
-
double I¶
Moment of inertia.
-
double W¶
Gravitational binding energy.
-
double Z_p¶
Polar redshift.
-
double Z_f¶
Forward equatorial redshift.
-
double Z_b¶
Backward equatorial redshift.
-
double Omega_K¶
Kepler rotation frequency (in 1/s)
-
double eccentricity¶
The eccentricity.
-
double vel_plus¶
Desc.
-
double vel_minus¶
Desc.
-
double h_plus¶
Height from surface of last stable co-rotating circular orbit in equatorial plane.
If this is zero then all orbits are stable.
-
double h_minus¶
Height from surface of last stable counter-rotating circular orbit in equatorial plane.
If this is zero then all orbits are stable.
-
double Omega_plus¶
Desc.
-
double u_phi¶
Desc.
-
double Omega_p¶
Angular velocity of a particle in a circular orbit at the equator.
-
double grv2¶
Desc.
-
double grv2_new¶
Desc.
-
double grv3¶
Desc.
-
double om_over_Om¶
Ratio of potential \( \omega \) to angular velocity \( \Omega \).
-
double mass_quadrupole¶
Mass quadrupole moment.
Settings
-
double cf¶
The convergence factor (default 1.0)
Internal constants
-
double C¶
Speed of light in vacuum (in CGS units)
-
double G¶
Gravitational constant (in CGS units)
-
double MSUN¶
Mass of sun (in g)
-
double KAPPA¶
Square of length scale in CGS units, \( \kappa \equiv 10^{-15} c^2/G \).
-
double MB¶
The mass of one baryon (in g)
-
double KSCALE¶
The value \( \kappa G c^{-4} \).
-
double PI¶
The constant \( \pi \).
-
void constants_rns()¶
Use the values of the constants from the original RNS code.
-
void constants_o2scl()¶
Use the O2scl constants (the default)
Grid quantities set in make_grid()
Grid values computed in integrate() for spherical_star()
Metric functions
Initial guess computed by the comp() function
-
double r_e_guess¶
Guess for the equatorial radius.
EOS quantities
Other quantities defined over the full two-dimensional grid
Quantities defined for fixed values of mu
-
double gamma_pole_h¶
The value of \( \hat{\gamma} \) at the pole.
-
double gamma_center_h¶
The value of \( \hat{\gamma} \) at the center.
-
double gamma_equator_h¶
The value of \( \hat{\gamma} \) at the equator.
-
double rho_pole_h¶
The value of \( \hat{\rho} \) at the pole.
-
double rho_center_h¶
The value of \( \hat{\rho} \) at the center.
-
double rho_equator_h¶
The value of \( \hat{\rho} \) at the equator.
-
double omega_equator_h¶
The value of \( \hat{\omega} \) at the equator.
-
double Omega_h¶
Angular velocity, \( \hat{\omega} \).
-
double p_center¶
Central pressure.
-
double h_center¶
Central enthalpy.
Helper functions for Green’s function expansion
Legendre polynomials
-
double tol_abs¶
The tolerance for the functions with the prefix “fix” (default \( 10^{-4} \) )
Thermodyanmic quantities near the surface
-
double p_surface¶
Pressure at the surface.
-
double e_surface¶
Energy density at the surface.
-
double enthalpy_min¶
Minimum specific enthalpy.
Interpolation functions
-
int n_nearest¶
Cache for interpolation.
-
double interp(ubvector &xp, ubvector &yp, int np, double xb)¶
Driver for the interpolation routine.
First we find the tab. point nearest to xb, then we interpolate using four points around xb.
Used by int_z(), e_at_p(), p_at_e(), p_at_h(), h_at_p(), n0_at_e(), comp_omega(), comp_M_J(), comp(), spherical_star(), iterate().
Basic Usage
-
inline void set_eos(eos_nstar_rot &eos)¶
Set the EOS.
-
inline void polytrope_eos(double index)¶
Use a polytropic EOS with a specified index.
-
int fix_cent_eden_axis_rat(double cent_eden, double axis_rat, bool use_guess = false)¶
Construct a configuration with a fixed central energy density and a fixed axis ratio.
The central energy density should be in \( \mathrm{g}/\mathrm{cm}^3 \) and the axis ratio is unitless. This is fastest of the high-level interface functions as it doesn’t require an additional solver.
-
int fix_cent_eden_grav_mass(double cent_eden, double grav_mass)¶
Construct a configuration with a fixed central energy density and a fixed gravitational mass.
The central energy density should be in \( \mathrm{g}/\mathrm{cm}^3 \) and the gravitational mass should be in solar masses.
-
int fix_cent_eden_bar_mass(double cent_eden, double bar_mass)¶
Construct a configuration with a fixed central energy density and a fixed baryonic mass.
The central energy density should be in \( \mathrm{g}/\mathrm{cm}^3 \) and the baryonic mass should be in solar masses.
-
int fix_cent_eden_with_kepler(double cent_eden)¶
Construct a configuration with a fixed central energy density and the Keplerian rotation rate.
The central energy density should be in \( \mathrm{g}/\mathrm{cm}^3 \) .
-
int fix_cent_eden_with_kepler_alt(double cent_eden, bool use_guess = false)¶
Experimental alternate form for fix_cent_eden_with_kepler()
-
int fix_cent_eden_grav_mass_alt(double cent_eden, double grav_mass, bool use_guess = false)¶
Experimental alternate form for fix_cent_eden_grav_mass()
-
int fix_cent_eden_bar_mass_alt(double cent_eden, double bar_mass, bool use_guess = false)¶
Experimental alternate form for fix_cent_eden_bar_mass()
-
int fix_cent_eden_ang_vel_alt(double cent_eden, double ang_vel, bool use_guess = false)¶
Experimental alternate form for fix_cent_eden_ang_vel()
-
int fix_cent_eden_ang_mom_alt(double cent_eden, double ang_mom, bool use_guess = false)¶
Experimental alternate form for fix_cent_eden_ang_mom()
-
int fix_cent_eden_non_rot(double cent_eden)¶
Construct a non-rotating configuration with a fixed central energy density.
The central energy density should be in \( \mathrm{g}/\mathrm{cm}^3 \) .
-
int fix_cent_eden_ang_vel(double cent_eden, double ang_vel)¶
Construct a configuration with a fixed central energy density and a fixed angular velocity.
The central energy density should be in \( \mathrm{g}/\mathrm{cm}^3 \) and the angular velocity should be in \( \mathrm{rad}/\mathrm{s} \). The final angular velocity (possibly slightly different than
ang_vel
is stored in Omega .Note
In the original RNS code, the
ang_vel
argument is different because it was rescaled by a factor of \( 10^{4} \).
-
int fix_cent_eden_ang_mom(double cent_eden, double ang_mom)¶
Construct a configuration with a fixed central energy density and a fixed angular momentum.
The central energy density should be in \( \mathrm{g}/\mathrm{cm}^3 \). The angular momentum should be in units of \( G M_{\odot}^2/C \) .
Testing functions
These functions compare this class with hard-coded results obtained with the RNS code.
-
void test1(o2scl::test_mgr &t)¶
Test determining configuration with fixed central energy density and fixed radius ratio with EOS C.
-
void test2(o2scl::test_mgr &t)¶
Test configuration rotating and Keplerian frequency with a fixed central energy density and EOS C.
-
void test3(o2scl::test_mgr &t)¶
Test fixed central energy density and fixed gravitational mass with EOS C.
-
void test4(o2scl::test_mgr &t)¶
Test fixed central energy density and fixed baryonic mass with EOS C.
-
void test5(o2scl::test_mgr &t)¶
Test fixed central energy density and fixed angular velocity with EOS C.
-
void test6(o2scl::test_mgr &t)¶
Test fixed central energy density and fixed angular momentum with EOS C.
EOS functions
-
double e_at_p(double pp)¶
Compute \( \varepsilon(P) \)
Used in dm_dr_is(), dp_dr_is(), integrate() and iterate().
-
double p_at_e(double ee)¶
Compute \( P(\varepsilon) \)
Used in make_center() and integrate().
-
double h_at_p(double pp)¶
Enthalpy at fixed pressure.
Used in make_center() and integrate().
-
double n0_at_e(double ee)¶
Baryon density at fixed energy density.
Used in comp_M_J() and comp() .
Derivatives on the grid
Initialization functions
-
double legendre(int n, double x)¶
Returns the Legendre polynomial of degree n, evaluated at x.
This uses the recurrence relation and is used in comp_f_P() which is called by the constructor.
-
void comp_f_P()¶
Compute two-point functions.
This function computes the 2-point functions \( f^m_{2n}(r,r') \) used to integrate the potentials \( \rho, \gamma \) and \( \omega \) (See Komatsu et al. 1989 for details). Since the grid points are fixed, we can compute the functions f_rho, f_gamma, f_omega, P_2n, and P1_2n_1 once at the beginning.
See Eqs. 27-29 of [Cook92] and Eqs. 33-35 of [Komatsu89]. This function is called by the constructor.
-
void make_grid()¶
Create computational mesh.
Create the computational mesh for \( s=r/(r+r_e) \) (where \( r_e \) is the coordinate equatorial radius) and \( \mu = \cos \theta \) using
\[ s[i]=\mathrm{SMAX}\left(\frac{i-1}{\mathrm{SDIV}-1}\right) \]\[ \mu[j]=\left(\frac{i-1}{\mathrm{MDIV}-1}\right) \]When \( r=0 \), \( s=0 \), when \( r=r_e \), \( s=1/2 \), and when \( r = \infty \), \( s=1 \) . Inverting the relationship between \( r \) and \( s \) gives \( r = r_e s / (1-s) \) .Points in the mu-direction are stored in the array
mu[i]
. Points in the s-direction are stored in the arrays_gp[j]
.This function sets s_gp, s_1_s, one_s, mu, one_m2, theta and sin_theta . All of these arrays are unit-indexed. It is called by the constructor.
-
void make_center(double e_center)¶
Compute central pressure and enthalpy from central energy density.
For polytropic EOSs, this also computes
rho0_center
.
Post-processing functions
-
void comp_omega()¶
Compute Omega and Omega_K.
-
void comp_M_J()¶
Compute rest mass and angular momentum.
-
void comp()¶
Compute various quantities.
The main post-processing function
For computing spherical stars
-
void spherical_star()¶
Computes a spherically symmetric star.
The metric is
\[ ds^2 = -e^{2\nu}dt^2 + e^{2\lambda} dr^2 + r^2 d\theta^2 + r^2 sin^2\theta d\phi^2 \]where \( r \) is an isotropic radial coordinate (corresponding tor_is
in the code).- Todo:
AWS: 7/20/20: Better document out how this metric definition leads to \( \gamma=\nu+\lambda \) and \( \rho = \nu - \lambda \) and the relationship between r and r_is .
-
double dm_dr_is(double r_is, double r, double m, double p)¶
Derivative of gravitational mass with respect to isotropic radius.
-
double dp_dr_is(double r_is, double r, double m, double p)¶
Derivative of pressure with respect to isotropic radius.
-
double dr_dr_is(double r_is, double r, double m)¶
Derivative of radius with respect to isotropic radius.
-
void integrate(int i_check, double &r_final, double &m_final, double &r_is_final)¶
Integrate one of the differential equations for spherical stars.
Desc
-
int iterate(double r_ratio, double tol_rel)¶
Main iteration function.
Public Types
-
typedef boost::numeric::ublas::vector<double> ubvector¶
-
typedef boost::numeric::ublas::range ub_range¶
-
typedef boost::numeric::ublas::vector_range<boost::numeric::ublas::vector<double>> ubvector_range¶
-
typedef boost::numeric::ublas::matrix<double> ubmatrix¶
Public Functions
-
void resize(int MDIV_new, int SDIV_new, int LMAX_new, int RDIV_new)¶
Resize the grid.
-
inline int set_solver(o2scl::mroot<> &m)¶
Set new solver.
Public Members
-
int MDIV¶
The number of grid points in the \( \mu \) direction.
-
int SDIV¶
The number of grid points in the \( s \) direction.
-
int LMAX¶
The number of Legendre polynomials.
-
o2scl::mroot_hybrids def_mroot¶
Default solver.
Protected Functions
-
int solve_grav_mass(size_t nv, const ubvector &x, ubvector &y, double grav_mass)¶
Solve for the gravitational mass.
-
int solve_bar_mass(size_t nv, const ubvector &x, ubvector &y, double bar_mass)¶
Solve for the gravitational mass.
Protected Attributes
-
o2scl::mroot *mrootp¶
Solver.
-
o2scl::root_bkt_cern<polytrope_solve> rbc¶
The polytrope solver.
-
o2scl::search_vec<ubvector_range> sv_ub¶
Array search object.
-
int RDIV¶
The number of grid points in integration of TOV equations for spherical stars.
-
double SMAX¶
Maximum value of s-coordinate (default 0.9999)
-
double DS¶
Spacing in \( s \) direction, \( \mathrm{SMAX}/(\mathrm{SDIV}-1) \).
-
double DM¶
Spacing in \( \mu \) direction, \( 1/(\mathrm{MDIV}-1) \).
-
double RMIN¶
Minimum radius for spherical stars (default \( 10^{-15} \))