Class tov_love (o2scl)

O2scl : Class List

class 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 Postnikov et al. (2010) )

\[\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,

\[\begin{split}\begin{eqnarray*} 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} \right. \nonumber \\ && \left. + 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} \nonumber \\ && \times \left[ -1+\frac{G m}{r} + 2 \pi G r^2 \left(\varepsilon-P\right) \right] \, . \end{eqnarray*}\end{split}\]

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], [Postnikov10], and [Hinderer10]).

\[ 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} \]

Idea for Future:

Improve calc_H() to handle discontinuities and to tabulate the EOS.

Improve the handling at small r using an expansion, similar to that used in e.g. Detweiler and Lindblom (1985)?

Allow specification of, e.g., an eos_tov like object rather than an EOS table.

Warning

The Love number can be sensitive to inaccuracy in the input EOS table, thus it’s important to make sure the EOS table has a sufficiently fine grid to give accurate results.

Note

The function calc_H() cannot yet handle discontinuities (if there are any then the error handler is called).

Public Types

typedef std::function<int(double, size_t, const std::vector<double>&, std::vector<double>&)> ode_funct2

The ODE function type.

Public Functions

tov_love()

inline 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.

The initial values of yR, beta, k2, lambda_km5, and lambda_cgs are ignored.

If tabulate is true, then the results of the ODE integration are stored in results . The o2scl::table::clear() function is called beforehand, so any data stored in in results is destroyed.

Note

The final results for the tidal deformability may differ when tabulate is true versus when tabulate is false. This difference is likely within the uncertainty of the ODE integration.

inline void add_disc(double rd)

Add a discontinuity at radius rd (in km)

inline 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.

The initial values of yR, beta, k2, lambda_km5, and lambda_cgs are ignored.

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.

This table is filled when calc_y() is called with tabulate=true.

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

A pointer to the ODE integrator.

double schwarz_km

Schwarzchild radius in km (set in constructor)

std::vector<double> disc

List of discontinuities.