Class tov_love (o2scl)¶
-
class
o2scl
::
tov_love
¶ Determination of the neutron star Love number.
We use \( c=1 \) but keep factors of \( G \), which has units \( \mathrm{km}/\mathrm{M_{\odot}} \).
Following the notation in Postnikov10, define the function \( H(r) \), which is the solution of
\[ H^{\prime \prime} (r) + H^{\prime}(r) \left\{ \frac{2}{r} + e^{\lambda(r)} \left[ \frac{2 G m(r)}{r^2} + 4 \pi G r P(r) - 4 \pi G r \varepsilon(r) \right] \right\} + H(r) Q(r) = 0 \]where (now supressing the dependence on \( r \)),\[ \nu^{\prime} \equiv 2 G e^{\lambda} \left(\frac{m+4 \pi P r^3}{r^2}\right) \, \]which has units of \( 1/\mathrm{km} \) ,\[ e^{\lambda} \equiv \left(1-\frac{2 G m}{r}\right)^{-1} \, , \]and\[ Q \equiv 4 \pi G e^{\lambda} \left( 5 \varepsilon + 9 P + \frac{\varepsilon+P}{c_s^2}\right) - 6 \frac{e^{\lambda}}{r^2} - \nu^{\prime 2} \]which has units of \( 1/\mathrm{km}^2 \) . The boundary conditions on \( H(r) \) are that \( H(r) = a_0 r^2 \) and \( H^{\prime} = 2 a_0 r \) for an arbitrary constant \( a_0 \) ( \( a_0 \) is chosen to be equal to 1). Internally, \( P \) and \( \varepsilon \) are stored in units of \( \mathrm{M}_{\odot}/\mathrm{km}^3 \) .From this we can define another (unitless) function \( y(r) \equiv r H^{\prime}(r)/H(r) \), which obeys
\[ r y^{\prime} + y^2 + y e^{\lambda} \left[ 1+4 \pi G r^2 \left( P-\varepsilon \right) \right] + r^2 Q = 0 \]with boundary condition is \( y(0) = 2 \) . Solving for \( y^{\prime} \),\[ y^{\prime} = \frac{1}{r} \left\{-r^2 Q - y e^{\lambda} \left[ 1+ 4 \pi G r^2 \left( P - \varepsilon \right) \right] -y^2 \right\} \]Define \( y_R = y(r=R) \). This form for \( y^{\prime}(r) \) is specified in y_derivs() .The unitless quantity \( k_2[\beta,y_R] \) (the Love number) is defined by (this is the expression from Postnikov10 )
\[\begin{split}\begin{eqnarray*} k_2[\beta,y(r=R)] &\equiv& \frac{8}{5} \beta^5 \left(1-2 \beta\right)^2 \left[ 2 - y_R + 2 \beta \left( y_R - 1 \right) \right] \nonumber \\ && \times \left\{ 2 \beta \left( 6 - 3 y_R + 3 \beta ( 5 y_R - 8) + 2 \beta^2 \left[ 13 - 11 y_R + \beta (3 y_R-2) \right.\right.\right. \nonumber \\ && + \left.\left.\left. 2 \beta^2 (1+y_R) \right] \right) + 3 (1-2 \beta)^2 \left[ 2 - y_R + 2 \beta (y_R - 1) \right] \log (1-2 \beta) \right\}^{-1} \end{eqnarray*}\end{split}\]Hinderer10 writes the differential equation for \( H(r) \) in a slightly different (but equivalent) form,
\[ H^{\prime \prime}(r) = 2 \left( 1 - \frac{2 G m}{r}\right)^{-1} H(r) \left\{ - 2 \pi G \left[ 5 \varepsilon + 9 P + \frac{\left( \varepsilon+P\right)}{c_s^2} \right] + \frac{3}{r^2} + 2 \left( 1 - \frac{2 G m}{r}\right)^{-1} \left(\frac{G m}{r^2}+4 \pi G r P\right)^2 \right\} +\frac{2 H^{\prime}(r)}{r} \left( 1 - \frac{2 G m}{r}\right)^{-1} \left[ -1+\frac{G m}{r} + 2 \pi G r^2 \left(\varepsilon-P\right) \right] \, . \]This is the form given in H_derivs() .The tidal deformability is then
\[ \lambda \equiv \frac{2}{3} k_2 R^5 \]and has units of \( \mathrm{km}^5 \) or can be converted to \( \mathrm{g}~\mathrm{cm}^2~\mathrm{s}^2 \) .It is assumed that tab stores a stellar profile (such as one computed with tov_solve::fixed(), tov_solve::fixed_pr(), or tov_solve::max() ) has been specified before-hand and contains (at least) the following columns
ed
energy density in units of \( \mathrm{M}_{\odot}/\mathrm{km}^3 \)pr
pressure in units of \( \mathrm{M}_{\odot}/\mathrm{km}^3 \)cs2
sound speed squared (unitless)gm
gravitational mass in \( \mathrm{M}_{\odot} \)r
radius in \( \mathrm{km} \) (Note that the o2scl::tov_solve class doesn’t automatically compute the columncs2
.)
This class handles the inner boundary by starting from the small non-zero radius stored in eps instead of at \( r=0 \). The value of eps defaults to 0.2 km.
If there is a discontinuity in the EOS (i.e. a jump in the energy density at some radius \( r_d \)), then the function \( y(r) \) must satisfy (see Damour09 and Postnikov10)
\[ y(r_d+\delta) - y(r_d-\delta) = \frac{ \varepsilon(r_d+\delta)-\varepsilon(r_d-\delta)}{m(r_d)/(4 \pi r_d^3) + p} \](See [Damour09] and [Postnikov10])
- Note
The function calc_H() cannot yet handle discontinuities (if there are any then the error handler is called).
- Note
This class does not yet handle strange quark stars which have an energy discontinuity at the surface (this currently causes the error handler to be called).
Public Types
-
typedef std::function<int(double, size_t, const std::vector<double>&, std::vector<double>&)>
ode_funct2
¶
Public Functions
-
tov_love
()¶
-
void
set_ODE
(o2scl::ode_iv_solve<ode_funct2, std::vector<double>> &ois_new)¶ Set ODE integrator.
-
int
calc_y
(double &yR, double &beta, double &k2, double &lambda_km5, double &lambda_cgs, bool tabulate = false)¶ Compute the love number using y.
-
void
add_disc
(double rd)¶ Add a discontinuity at radius
rd
(in km)
-
void
clear_discs
()¶ Remove all discontinuities.
-
int
calc_H
(double &yR, double &beta, double &k2, double &lambda_km5, double &lambda_cgs)¶ Compute the love number using H.
Public Members
-
int
show_ode
¶ If greater than zero, show the ODE output (default 0)
-
bool
addl_testing
¶ Additional testing if the ODE solver fails.
-
bool
err_nonconv
¶ If true, call the error handler if the solution does not converge (default true)
-
o2scl::table_units
results
¶ A table containing the solution to the differential equation(s)
-
double
delta
¶ The radial step for resolving discontinuities in km (default \( 10^{-4} \))
-
double
eps
¶ The first radial point in \( \mathrm{km} \) (default 0.02)
-
o2scl::ode_iv_solve<ode_funct2, std::vector<double>>
def_ois
¶ The default ODE integrator.
-
std::shared_ptr<o2scl::table_units<>>
tab
¶ Pointer to the input profile.
Protected Functions
-
int
y_derivs
(double r, size_t nv, const std::vector<double> &vals, std::vector<double> &ders)¶ The derivative \( y^{\prime}(r) \).
-
int
H_derivs
(double r, size_t nv, const std::vector<double> &vals, std::vector<double> &ders)¶ The derivatives \( H^{\prime \prime}(r) \) and \( H^{\prime}(r) \).
-
double
eval_k2
(double beta, double yR)¶ Compute \( k_2(\beta,y_R) \) using the analytic expression.
Used in both tov_love::calc_y() and tov_love::calc_H().
Protected Attributes
-
o2scl::ode_iv_solve<ode_funct2, std::vector<double>> *
oisp
¶ The ODE integrator.
-
double
schwarz_km
¶ Schwarzchild radius in km (set in constructor)
-
std::vector<double>
disc
¶ List of discontinuities.