tov_love.h
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-2018, 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  \rho(r_d+\delta)-\rho(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  int show_ode;
213 
214  /** \brief Additional testing
215  */
217 
218  /** \brief If true, call the error handler if the solution does
219  not converge (default true)
220  */
222 
223  /// A table containing the solution to the differential equation(s)
225 
226  /** \brief The radial step for resolving discontinuities in km
227  (default \f$ 10^{-4} \f$)
228  */
229  double delta;
230 
231  /// The first radial point in \f$ \mathrm{km} \f$ (default 0.02)
232  double eps;
233 
234  /// The default ODE integrator
236 
237  /// Pointer to the input profile
238  std::shared_ptr<o2scl::table_units<> > tab;
239 
240  /// Set ODE integrator
241  void set_ODE(o2scl::ode_iv_solve<ode_funct2,std::vector<double> >
242  &ois_new) {
243  oisp=&ois_new;
244  }
245 
246  /** \brief Compute the love number using y
247  */
248  int calc_y(double &yR, double &beta, double &k2, double &lambda_km5,
249  double &lambda_cgs, bool tabulate=false);
250 
251  /** \brief Add a discontinuity at radius \c rd (in km)
252  */
253  void add_disc(double rd) {
254  disc.push_back(rd);
255  }
256 
257  /** \brief Remove all discontinuities
258  */
259  void clear_discs() {
260  disc.clear();
261  }
262 
263  /** \brief Compute the love number using H
264  */
265  int calc_H(double &yR, double &beta, double &k2, double &lambda_km5,
266  double &lambda_cgs);
267 
268  };
269 
270 #ifndef DOXYGEN_NO_O2NS
271 }
272 #endif
273 
274 #endif
int y_derivs(double r, size_t nv, const std::vector< double > &vals, std::vector< double > &ders)
The derivative .
double delta
The radial step for resolving discontinuities in km (default )
Definition: tov_love.h:229
std::vector< double > disc
List of discontinuities.
Definition: tov_love.h:204
double schwarz_km
Schwarzchild radius in km (set in constructor)
Definition: tov_love.h:193
bool err_nonconv
If true, call the error handler if the solution does not converge (default true)
Definition: tov_love.h:221
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::ode_iv_solve< ode_funct2, std::vector< double > > * oisp
The ODE integrator.
Definition: tov_love.h:179
void set_ODE(o2scl::ode_iv_solve< ode_funct2, std::vector< double > > &ois_new)
Set ODE integrator.
Definition: tov_love.h:241
double eval_k2(double beta, double yR)
Compute using the analytic expression.
Determination of the neutron star Love number.
Definition: tov_love.h:166
int calc_H(double &yR, double &beta, double &k2, double &lambda_km5, double &lambda_cgs)
Compute the love number using H.
void add_disc(double rd)
Add a discontinuity at radius rd (in km)
Definition: tov_love.h:253
int H_derivs(double r, size_t nv, const std::vector< double > &vals, std::vector< double > &ders)
The derivatives and .
double eps
The first radial point in (default 0.02)
Definition: tov_love.h:232
o2scl::table_units results
A table containing the solution to the differential equation(s)
Definition: tov_love.h:224
std::shared_ptr< o2scl::table_units<> > tab
Pointer to the input profile.
Definition: tov_love.h:238
bool addl_testing
Additional testing.
Definition: tov_love.h:216
void clear_discs()
Remove all discontinuities.
Definition: tov_love.h:259
o2scl::ode_iv_solve< ode_funct2, std::vector< double > > def_ois
The default ODE integrator.
Definition: tov_love.h:235

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