Class tov_love (o2scl)

O2scl_eos : Class List

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 column cs2.)

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.