tov_love.h
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-2020, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef TOV_LOVE_H
24 #define TOV_LOVE_H
25 
26 #include <gsl/gsl_math.h>
27 
28 #include <boost/numeric/ublas/vector.hpp>
29 
30 #include <o2scl/astep_gsl.h>
31 #include <o2scl/table_units.h>
32 #include <o2scl/astep_nonadapt.h>
33 #include <o2scl/ode_rk8pd_gsl.h>
34 #include <o2scl/ode_iv_solve.h>
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /** \brief Determination of the neutron star Love number
41 
42  We use \f$ c=1 \f$ but keep factors of \f$ G \f$, which has
43  units \f$ \mathrm{km}/\mathrm{M_{\odot}} \f$.
44 
45  Following the notation in \ref Postnikov10, define
46  the function \f$ H(r) \f$, which is the solution of
47  \f[
48  H^{\prime \prime} (r) + H^{\prime}(r) \left\{
49  \frac{2}{r} + e^{\lambda(r)} \left[ \frac{2 G m(r)}{r^2} +
50  4 \pi G r P(r) - 4 \pi G r \varepsilon(r) \right]
51  \right\} + H(r) Q(r) = 0
52  \f]
53  where (now supressing the dependence on \f$ r \f$),
54  \f[
55  \nu^{\prime} \equiv 2 G e^{\lambda}
56  \left(\frac{m+4 \pi P r^3}{r^2}\right) \,
57  \f]
58  which has units of \f$ 1/\mathrm{km} \f$ ,
59  \f[
60  e^{\lambda} \equiv \left(1-\frac{2 G m}{r}\right)^{-1} \, ,
61  \f]
62  and
63  \f[
64  Q \equiv 4 \pi G e^{\lambda} \left( 5 \varepsilon + 9 P +
65  \frac{\varepsilon+P}{c_s^2}\right) - 6 \frac{e^{\lambda}}{r^2}
66  - \nu^{\prime 2}
67  \f]
68  which has units of \f$ 1/\mathrm{km}^2 \f$ .
69  The boundary conditions on \f$ H(r) \f$ are that \f$ H(r) = a_0
70  r^2 \f$ and \f$ H^{\prime} = 2 a_0 r \f$ for an arbitrary
71  constant \f$ a_0 \f$ (\f$ a_0 \f$ is chosen to be equal to 1).
72  Internally, \f$ P \f$ and \f$ \varepsilon \f$ are stored in
73  units of \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$ .
74 
75  From this we can define another (unitless) function
76  \f$ y(r) \equiv r H^{\prime}(r)/H(r) \f$, which obeys
77  \f[
78  r y^{\prime} + y^2 + y e^{\lambda} \left[
79  1+4 \pi G r^2 \left( P-\varepsilon \right)
80  \right] + r^2 Q = 0
81  \f]
82  with boundary condition is \f$ y(0) = 2 \f$ .
83  Solving for \f$ y^{\prime} \f$,
84  \f[
85  y^{\prime} = \frac{1}{r}
86  \left\{-r^2 Q - y e^{\lambda} \left[ 1+ 4 \pi G r^2
87  \left( P - \varepsilon \right) \right] -y^2 \right\}
88  \f]
89  Define \f$ y_R = y(r=R) \f$. This form for \f$ y^{\prime}(r) \f$
90  is specified in \ref y_derivs() .
91 
92  The unitless quantity \f$ k_2[\beta,y_R] \f$ (the Love number)
93  is defined by (this is the expression from \ref Postnikov10 )
94  \f{eqnarray*}
95  k_2[\beta,y(r=R)] &\equiv& \frac{8}{5} \beta^5
96  \left(1-2 \beta\right)^2
97  \left[ 2 - y_R + 2 \beta \left( y_R - 1 \right) \right]
98  \nonumber \\
99  && \times \left\{ 2 \beta \left( 6 - 3 y_R + 3 \beta ( 5 y_R - 8)
100  + 2 \beta^2 \left[ 13 - 11 y_R + \beta (3 y_R-2)
101  \right.\right.\right. \nonumber \\
102  && + \left.\left.\left. 2 \beta^2 (1+y_R) \right] \right) + 3
103  (1-2 \beta)^2 \left[ 2 - y_R + 2 \beta (y_R - 1) \right]
104  \log (1-2 \beta) \right\}^{-1}
105  \f}
106 
107  \ref Hinderer10 writes the differential equation for \f$ H(r) \f$
108  in a slightly different (but equivalent) form,
109  \f[
110  H^{\prime \prime}(r) = 2 \left( 1 - \frac{2 G m}{r}\right)^{-1}
111  H(r) \left\{ - 2 \pi G \left[ 5 \varepsilon + 9 P + \frac{\left(
112  \varepsilon+P\right)}{c_s^2} \right] + \frac{3}{r^2} +
113  2 \left( 1 - \frac{2 G m}{r}\right)^{-1}
114  \left(\frac{G m}{r^2}+4 \pi G r P\right)^2 \right\}
115  +\frac{2 H^{\prime}(r)}{r} \left( 1 - \frac{2 G m}{r}\right)^{-1}
116  \left[ -1+\frac{G m}{r} + 2 \pi G r^2 \left(\varepsilon-P\right)
117  \right] \, .
118  \f]
119  This is the form given in \ref H_derivs() .
120 
121  The tidal deformability is then
122  \f[
123  \lambda \equiv \frac{2}{3} k_2 R^5
124  \f]
125  and has units of \f$ \mathrm{km}^5 \f$ or can be converted to
126  \f$ \mathrm{g}~\mathrm{cm}^2~\mathrm{s}^2 \f$ .
127 
128  It is assumed that \ref tab stores a stellar profile (such as
129  one computed with \ref tov_solve::fixed(), \ref
130  tov_solve::fixed_pr(), or \ref tov_solve::max() ) has been
131  specified before-hand and contains (at least) the following
132  columns
133  - <tt>ed</tt> energy density in units of
134  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$
135  - <tt>pr</tt> pressure in units of
136  \f$ \mathrm{M}_{\odot}/\mathrm{km}^3 \f$
137  - <tt>cs2</tt> sound speed squared (unitless)
138  - <tt>gm</tt> gravitational mass in \f$ \mathrm{M}_{\odot} \f$
139  - <tt>r</tt> radius in \f$ \mathrm{km} \f$
140  (Note that the \ref o2scl::tov_solve class doesn't automatically compute
141  the column <tt>cs2</tt>.)
142 
143  This class handles the inner boundary by starting from the small
144  non-zero radius stored in \ref eps instead of at \f$ r=0 \f$. The
145  value of \ref eps defaults to 0.2 km.
146 
147  If there is a discontinuity in the EOS (i.e. a jump in
148  the energy density at some radius \f$ r_d \f$), then
149  the function \f$ y(r) \f$ must satisfy (see \ref Damour09
150  and \ref Postnikov10)
151  \f[
152  y(r_d+\delta) - y(r_d-\delta) =
153  \frac{
154  \varepsilon(r_d+\delta)-\varepsilon(r_d-\delta)}{m(r_d)/(4 \pi r_d^3) +
155  p}
156  \f]
157 
158  \note The function \ref calc_H() cannot yet handle
159  discontinuities (if there are any then the error handler
160  is called).
161 
162  \note This class does not yet handle strange quark
163  stars which have an energy discontinuity at the surface
164  (this currently causes the error handler to be called).
165  */
166  class tov_love {
167 
168  public:
169 
170  typedef std::function<int(double,size_t,
171  const std::vector<double> &,
172  std::vector<double> &)> ode_funct2;
173 
174 #ifndef DOXYGEN_INTERNAL
175 
176  protected:
177 
178  /// The ODE integrator
180 
181  /** \brief The derivative \f$ y^{\prime}(r) \f$
182  */
183  int y_derivs(double r, size_t nv, const std::vector<double> &vals,
184  std::vector<double> &ders);
185 
186  /** \brief The derivatives \f$ H^{\prime \prime}(r) \f$ and
187  \f$ H^{\prime}(r) \f$
188  */
189  int H_derivs(double r, size_t nv, const std::vector<double> &vals,
190  std::vector<double> &ders);
191 
192  /// Schwarzchild radius in km (set in constructor)
193  double schwarz_km;
194 
195  /** \brief Compute \f$ k_2(\beta,y_R) \f$ using the analytic
196  expression
197 
198  Used in both \ref tov_love::calc_y() and \ref
199  tov_love::calc_H().
200  */
201  double eval_k2(double beta, double yR);
202 
203  /// List of discontinuities
204  std::vector<double> disc;
205 
206 #endif
207 
208  public:
209 
210  tov_love();
211 
212  /** \brief If greater than zero, show the ODE output (default 0)
213  */
214  int show_ode;
215 
216  /** \brief Additional testing if the ODE solver fails
217  */
219 
220  /** \brief If true, call the error handler if the solution does
221  not converge (default true)
222  */
224 
225  /// A table containing the solution to the differential equation(s)
227 
228  /** \brief The radial step for resolving discontinuities in km
229  (default \f$ 10^{-4} \f$)
230  */
231  double delta;
232 
233  /// The first radial point in \f$ \mathrm{km} \f$ (default 0.02)
234  double eps;
235 
236  /// The default ODE integrator
238 
239  /// Pointer to the input profile
240  std::shared_ptr<o2scl::table_units<> > tab;
241 
242  /// Set ODE integrator
243  void set_ODE(o2scl::ode_iv_solve<ode_funct2,std::vector<double> >
244  &ois_new) {
245  oisp=&ois_new;
246  }
247 
248  /** \brief Compute the love number using y
249  */
250  int calc_y(double &yR, double &beta, double &k2, double &lambda_km5,
251  double &lambda_cgs, bool tabulate=false);
252 
253  /** \brief Add a discontinuity at radius \c rd (in km)
254  */
255  void add_disc(double rd) {
256  disc.push_back(rd);
257  }
258 
259  /** \brief Remove all discontinuities
260  */
261  void clear_discs() {
262  disc.clear();
263  }
264 
265  /** \brief Compute the love number using H
266  */
267  int calc_H(double &yR, double &beta, double &k2, double &lambda_km5,
268  double &lambda_cgs);
269 
270  };
271 
272 #ifndef DOXYGEN_NO_O2NS
273 }
274 #endif
275 
276 #endif
o2scl::ode_iv_solve
o2scl::tov_love::def_ois
o2scl::ode_iv_solve< ode_funct2, std::vector< double > > def_ois
The default ODE integrator.
Definition: tov_love.h:237
o2scl::tov_love::y_derivs
int y_derivs(double r, size_t nv, const std::vector< double > &vals, std::vector< double > &ders)
The derivative .
o2scl::tov_love::calc_y
int calc_y(double &yR, double &beta, double &k2, double &lambda_km5, double &lambda_cgs, bool tabulate=false)
Compute the love number using y.
o2scl::tov_love::delta
double delta
The radial step for resolving discontinuities in km (default )
Definition: tov_love.h:231
o2scl::tov_love::add_disc
void add_disc(double rd)
Add a discontinuity at radius rd (in km)
Definition: tov_love.h:255
o2scl::tov_love::eval_k2
double eval_k2(double beta, double yR)
Compute using the analytic expression.
o2scl::tov_love::show_ode
int show_ode
If greater than zero, show the ODE output (default 0)
Definition: tov_love.h:214
o2scl::tov_love::addl_testing
bool addl_testing
Additional testing if the ODE solver fails.
Definition: tov_love.h:218
o2scl::tov_love::disc
std::vector< double > disc
List of discontinuities.
Definition: tov_love.h:204
o2scl::table_units
o2scl::tov_love::set_ODE
void set_ODE(o2scl::ode_iv_solve< ode_funct2, std::vector< double > > &ois_new)
Set ODE integrator.
Definition: tov_love.h:243
o2scl::tov_love::H_derivs
int H_derivs(double r, size_t nv, const std::vector< double > &vals, std::vector< double > &ders)
The derivatives and .
o2scl::tov_love::oisp
o2scl::ode_iv_solve< ode_funct2, std::vector< double > > * oisp
The ODE integrator.
Definition: tov_love.h:179
o2scl::tov_love::results
o2scl::table_units results
A table containing the solution to the differential equation(s)
Definition: tov_love.h:226
o2scl::tov_love::eps
double eps
The first radial point in (default 0.02)
Definition: tov_love.h:234
o2scl::tov_love::err_nonconv
bool err_nonconv
If true, call the error handler if the solution does not converge (default true)
Definition: tov_love.h:223
o2scl::tov_love::calc_H
int calc_H(double &yR, double &beta, double &k2, double &lambda_km5, double &lambda_cgs)
Compute the love number using H.
o2scl::tov_love::tab
std::shared_ptr< o2scl::table_units<> > tab
Pointer to the input profile.
Definition: tov_love.h:240
o2scl::tov_love::clear_discs
void clear_discs()
Remove all discontinuities.
Definition: tov_love.h:261
o2scl::tov_love
Determination of the neutron star Love number.
Definition: tov_love.h:166
o2scl::tov_love::schwarz_km
double schwarz_km
Schwarzchild radius in km (set in constructor)
Definition: tov_love.h:193

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).