tov_love.h
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-2017, 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 surpressing 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 has been specified before-hand
129  and has (at least) the following columns
130  - <tt>ed</tt> energy density in units of \f$ 1/\mathrm{fm}^4 \f$
131  - <tt>pr</tt> pressure in units of \f$ 1/\mathrm{fm}^4 \f$
132  - <tt>cs2</tt> sound speed squared (unitless)
133  - <tt>gm</tt> gravitational mass in \f$ \mathrm{M}_{\odot} \f$
134  - <tt>r</tt> radius in \f$ \mathrm{km} \f$
135 
136  \future Use \ref o2scl::ode_iv_solve instead of several steps of
137  type \ref o2scl::astep_gsl .
138 
139  */
140  class tov_love {
141 
142  public:
143 
145 
146 #ifndef DOXYGEN_INTERNAL
147 
148  protected:
149 
150  /// The default ODE integrator
152 
153  /// The ODE integrator
155 
156  /** \brief The derivative \f$ y^{\prime}(r) \f$
157  */
158  int y_derivs(double r, size_t nv, const ubvector &vals,
159  ubvector &ders);
160 
161  /** \brief The derivatives \f$ H^{\prime \prime}(r) \f$ and
162  \f$ H^{\prime}(r) \f$
163  */
164  int H_derivs(double r, size_t nv, const ubvector &vals,
165  ubvector &ders);
166 
167  /// Schwarzchild radius in km (set in constructor)
168  double schwarz_km;
169 
170  /** \brief Compute \f$ k_2(\beta,y_R) \f$ using the analytic
171  expression
172 
173  Used in both \ref tov_love::calc_y() and \ref
174  tov_love::calc_H().
175  */
176  double eval_k2(double beta, double yR);
177 
178 #endif
179 
180  public:
181 
182  tov_love();
183 
184  /// A table containing the solution to the differential equation(s)
186 
187  /// The first radial point in \f$ \mathrm{km} \f$ (default 0.02)
188  double eps;
189 
190  /// Pointer to the input profile
191  std::shared_ptr<o2scl::table_units<> > tab;
192 
193  /// Set ODE integrator
195  oisp=&ois_new;
196  }
197 
198  /** \brief Compute the love number using y
199  */
200  void calc_y(double &yR, double &beta, double &k2, double &lambda_km5,
201  double &lambda_cgs, bool tabulate=false);
202 
203  /** \brief Compute the love number using H
204  */
205  void calc_H(double &yR, double &beta, double &k2, double &lambda_km5,
206  double &lambda_cgs);
207 
208  };
209 
210 #ifndef DOXYGEN_NO_O2NS
211 }
212 #endif
213 
214 #endif
void calc_H(double &yR, double &beta, double &k2, double &lambda_km5, double &lambda_cgs)
Compute the love number using H.
double schwarz_km
Schwarzchild radius in km (set in constructor)
Definition: tov_love.h:168
int H_derivs(double r, size_t nv, const ubvector &vals, ubvector &ders)
The derivatives and .
void set_ODE(o2scl::ode_iv_solve<> &ois_new)
Set ODE integrator.
Definition: tov_love.h:194
void calc_y(double &yR, double &beta, double &k2, double &lambda_km5, double &lambda_cgs, bool tabulate=false)
Compute the love number using y.
double eval_k2(double beta, double yR)
Compute using the analytic expression.
int y_derivs(double r, size_t nv, const ubvector &vals, ubvector &ders)
The derivative .
Determination of the neutron star Love number.
Definition: tov_love.h:140
o2scl::ode_iv_solve def_ois
The default ODE integrator.
Definition: tov_love.h:151
double eps
The first radial point in (default 0.02)
Definition: tov_love.h:188
o2scl::ode_iv_solve * oisp
The ODE integrator.
Definition: tov_love.h:154
o2scl::table_units results
A table containing the solution to the differential equation(s)
Definition: tov_love.h:185
std::shared_ptr< o2scl::table_units<> > tab
Pointer to the input profile.
Definition: tov_love.h:191

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