Physics Overview

This section outlines the mathematical framework used to develop the physics for this module. The first part describes the process of deriving the equations that characterize isolated and slowly rotating neutron stars, explaining the key concepts and steps involved. The second part discusses the numerical implementation of these solutions, detailing the methods used to solve the equations and extract physical quantities.

The Hartle-Thorne method

The Hartle-Thorne procedure is a method for obtaining an approximate solution for slowly rotating neutron stars. It begins with an exact solution to Einstein’s equations for a non-rotating, spherically symmetric star, which serves as the background spacetime. Small disturbances are then introduced to this background in a slow-rotation expansion, controlled by a dimensionless expansion parameter (\(\epsilon\)). This parameter represents the ratio of the star’s angular speed (\(\Omega\)) to the critical angular speed at which mass shedding occurs (\(\Omega_{\textrm{shed}}\)). The model is based on the following assumptions:

1. Equation of State: In the non-rotating configuration, matter follows a barotropic equation of state represented as \(p=p(\varepsilon)\), where \(p\) denotes the pressure and \(\varepsilon\) stands for the total energy density. In conditions of extremely high densities as in the interior of NS, the Fermi temperature (\(T_{\textrm{f}} \sim 3 \times 10^{11} \, \textrm{K}\)) greatly surpasses the temperature of the star, resulting in an EoS which is unaffected by the star’s temperature. Thus, the main contribution to the pressure is due to nuclear forces and corrections caused by thermal agitation can be neglected.

2. Matter distribution: The matter distribution of the star in the non-rotating configuration is described by a perfect-fluid stress energy tensor. Shear stresses and energy transport are negligible.

3. Axial and reflection symmetry: The perturbed spacetime that characterizes the rotating star exhibits axial symmetry with respect to an arbitrary axis, which is assumed to be aligned with the star’s rotation axis. Moreover, the geometry remains unchanged when mirrored across a plane perpendicular to the axis of rotation. These symmetry properties allow for setting the structure of the geometry and enable the decoupling of Einstein’s Equations, as elaborated upon below.

4. Uniform rotation: Although differential rotation can be considered, uniformly rotating configurations are assumed, as they represent stable states that minimize the total mass-energy.

5. Slow rotation: The non-rotating metric configuration is perturbed in a dimensionless small spin parameter \(\epsilon = \Omega / \Omega_{\textrm{shed}}\). In this context, \(\Omega\) represents the angular speed of the neutron star as observed from spatial infinity, while \(\Omega_{\mathrm{shed}}\) is a characteristic frequency defined by \(\Omega^{2}_{\mathrm{shed}} \sim M_{\star}/R_{\star}^{3}\), where \(R_{\star}\) and \(M_{\star}\) refer to the radius and the mass in the non-rotating configuration. The quantity \(\Omega_{\mathrm{shed}}\) corresponds to the order of the Keplerian orbital frequency of a test particle at a radius \(R_{\star}\) moving around a mass \(M_{\star}\). Therefore, the approximation is valid for \(\epsilon \ll 1\), but breaks down when \(\epsilon \sim 1\).

Slow-rotation expansion

The slow-rotation expansion relies on the perturbation of the exact spherically symmetric background solution. In mathematical terms this means that the perturbed metric can be written as a power series

\[\begin{equation} g_{\alpha \beta} = g^{(0)}_{\alpha \beta} \, + \, \epsilon h^{(1)}_{\alpha \beta} \, + \, \frac{1}{2!} \epsilon^{2} h^{(2)}_{\alpha \beta} \, + \, \cdots \end{equation}\]

The quantity \(g^{(0)}_{\alpha \beta}\) is the static and spherically symmetric background metric defined as

\[\begin{equation} ds_{(0)}^{2} \equiv g^{(0)}_{\alpha \beta} dx^{\alpha}dx^{\beta} = -e^{\nu(r)} dt^{2} + e^{\lambda(r)} dr^{2} + r^{2}(d\theta^{2}+\sin^{2}\theta d\phi^{2}) \end{equation}\]

According to the symmetry assumptions, the perturbed metric can be expanded up to \(\mathcal{O}(\epsilon^{2})\) as []

\[\begin{split}\begin{align} ds^{2} =& -e^{\nu(r)} \big[ 1 + 2 \epsilon^{2}h^{(2)}(r,\theta) \big] dt^{2} + e^{\lambda(r)} \left[ 1 + \dfrac{2 \epsilon^{2} m^{(2)}(r,\theta)}{r - 2 M(r)} \right] dr^{2} \\[1ex] &+ r^{2} \big[ 1 + 2 \epsilon^{2}k^{(2)}(r,\theta) \big] \bigg( d\theta^{2} + \sin^{2}\theta \big[ d\phi - \epsilon \omega^{(1)}(r,\theta) dt \big]^{2} \bigg) \end{align}\end{split}\]

where \(\omega^{(1)}\) denotes a first order correction while \(h^{(2)}\), \(m^{(2)}\) and \(k^{(2)}\) correspond to second order corrections in the spin parameter (\(\epsilon\)). The form of the metric is justified by symmetry arguments. Note that for an axially symmetric spacetime a transformation of the form \(\Omega \rightarrow - \Omega\) is equivalent to \(t \rightarrow -t\). Therefore, \(g_{t\phi}\) should contain only odd powers of \(\epsilon\) and \(g_{tt}\), \(g_{rr}\), \(g_{\theta\theta}\), \(g_{\phi\phi}\) only even powers. From this metric ansatz, the perturbations \(h^{(1)}_{\alpha \beta}\) and \(h^{(2)}_{\alpha \beta}\) are therefore,

\[\begin{split}\begin{align} h^{(1)}_{\alpha \beta} dx^{\alpha}dx^{\beta} &= 2 r^{2} \omega^{(1)}(r,\theta) \sin^{2}\theta dt d\phi \, , \, \\[1ex] h^{(2)}_{\alpha \beta} dx^{\alpha}dx^{\beta} &= - \bigg[ 4 e^{\nu(r)} h^{(2)}(r,\theta) + 2 r^{2} \sin^{2}\theta \omega^{(1)}(r, \theta)^{2} \bigg] dt^{2} + 4 e^{\lambda(r,\theta)} \frac{m^{(2)}(r,\theta)}{r-2 M(r)}dr^{2} \\[1ex] &+ 4 r^{2} k^{(2)}(r,\theta) \left( d\theta^{2} +\sin^{2}\theta d\phi^{2} \right) \, . \end{align}\end{split}\]

The slow rotation perturbations will also produce corrections to the stress-energy tensor which is considered to be a perfect fluid given by

\[\begin{equation} T_{\alpha \beta} = (\varepsilon + p )u_{\alpha}u_{\beta} + pg_{\alpha \beta} \end{equation}\]

where \(p\) is the pressure, \(\varepsilon\) is the total energy density and \(u^{\alpha}\) is the four-velocity of the fluid. Due to the metric perturbations, the pressure and energy density will also develop perturbations,

\[\begin{split}\begin{align} p(\epsilon) &= p^{(0)}(r) + \epsilon p^{(1)}(r,\theta) + \frac{1}{2!} \epsilon^{2}p^{(2)}(r,\theta) \, , \\[1ex] \varepsilon(\epsilon) &= \varepsilon^{(0)}(r) + \epsilon \varepsilon^{(1)}(r,\theta) + \frac{1}{2!} \epsilon^{2}\varepsilon^{(2)}(r,\theta) \\[1ex] \end{align}\end{split}\]

\[\]

and therefore

\[\begin{equation} T_{\alpha \beta} = T^{(0)}_{\alpha \beta} \, + \, \epsilon T^{(1)}_{\alpha \beta} \, + \, \frac{1}{2!} \epsilon^{2} T^{(2)}_{\alpha \beta} \, + \, \cdots \end{equation}\]

\[\]

However, as pointed out by Hartle and Thorne, careful is needed with the coordinate frame used when expanding those quantites in a slow rotation expansion because they won’t be valid throughout the entire star. That would be valid if the fractional changes in energy density and pressure at each point were small. However, this is not true close to the boundary of the star where the density may be finite where the non-rotating configuration vanishes. In order to overcome this issue, we choose a new radial coordinate \(R\) in the non-rotating configuration such that the isodensity and isobaric contours of the perturbed configuration are the same as the isodensity and isobaric contours defined at \(R\) in the non-rotating configuration. In matematical terms, this means

\[\begin{split}\begin{align} \varepsilon[r(R,\Theta), \theta] &= \varepsilon (R) = \varepsilon^{(0)}(R) \\[2ex] p[r(R,\Theta), \theta] &= p(R) = p^{(0)}(R) \end{align}\end{split}\]

\[\]

where the coordinate transformation between the old and new frame is given by

\[\begin{equation} r(R,\Theta) = R + \epsilon^{2} \xi^{(2)}(R,\Theta) \hspace{1cm} ; \hspace{1cm} \theta(R,\Theta) = \Theta \end{equation}\]

\[\]

Let’s call this new frame the Hartle-Thorne coordinate frame and this will be the frame to formally perform the perturbation expansions. By construction in this new frame the pressure and total energy density only contains the static part, i.e. no spin perturbations. Instead, the small displacement \(\xi^{(2)}\) plays the role of the density and pressure perturbations in the Hartle-Thorne frame. In this frame the line element reads up to \(\mathcal{O}(\epsilon^{2})\) as

\[\begin{split}\begin{align} ds^{2}_{\mathsf{HT}} &= - \left[ e^{\nu} \left( 1 + 2\epsilon^{2} h^{(2)} + \epsilon^{2}\xi^{(2)} \dfrac{d\nu}{dR} \right) - \epsilon^{2}R^{2} \omega^{(1)\, 2} \sin^{2}\Theta \right]dt^{2} - 2 \epsilon R^{2} \omega^{(1)}\sin^{2}\Theta dt d\phi \\[2ex] & \hspace{0.5cm} + \left[ R^{2}(1+2\epsilon^{2} k^{(2)}) + 2 \epsilon^{2} R \xi \right] \sin^2\Theta d\phi^{2} + 2\epsilon^{2} e^{\lambda} \dfrac{\partial \xi^{(2)}}{\partial \Theta} dR d\Theta \\[2ex] & \hspace{0.5cm} + e^{\lambda}\left( 1 + \epsilon^{2} \dfrac{2 m^{(2)}}{R-2 M} + \epsilon^{2} \xi^{(2)} \dfrac{d\lambda}{dR} + 2 \epsilon^{2} \dfrac{\partial \xi^{(2)}}{\partial R} \right) dR^{2} \\[2ex] & \hspace{0.5cm} + \left[ R^{2} \left( 1 + 2 \epsilon^{2} k^{(2)} \right) + 2\epsilon^{2} R \xi^{(2)} \right] d\Theta^{2} \end{align}\end{split}\]

\[\]

For a stationary and axisymmetric star, the four-velocity is given by

\[u^{\alpha} = (u^{0},0,0,\epsilon \Omega u^{0}).\]

\[\]

with \(\Omega\) being the angular speed of the star measured by an observer at rest at some point in the fluid and \(u^{0}\) is determined from the normalization condition \(u_{\alpha}u^{\alpha}=-1\). Using the metric in the Hartle-Thorne frame then,

\[u^{0} = e^{-\nu/2} + \dfrac{1}{2}\epsilon^{2} e^{-3\nu/2} \left[ (\Omega - \omega^{(1)})^{2}R^{2}\sin^{2} \Theta - e^{\nu} \left( 2h^{(2)} + \dfrac{d\nu}{dR} \xi^{(2)} \right) \right] + \mathcal{O}(\epsilon^{4})\]

\[\]

Once the background metric has been perturbed to the desired order in the Hartle-Thorne frame, one can systematically construct the perturbation equations order by order in the perturbation parameter from Einstein’s equations,

\[\begin{equation} E_{\mu\nu} := G_{\mu\nu} - 8\pi T_{\mu \nu} = 0 \ . \end{equation}\]

Since Einstein’s equations can be written as

\[\begin{equation} E_{\mu\nu}(\epsilon) = E_{\mu\nu}^{(0)} + E_{\mu\nu}^{(1)} \epsilon + \dfrac{1}{2} E_{\mu\nu}^{(2)} \epsilon^{2} + \cdots \end{equation}\]

every term must vanish independently. Therefore, the perturbation equations at the perturbation order \(n\) are given by

\[\begin{equation} E_{\mu\nu}^{(n)} = 0 \ . \end{equation}\]

These are partial differential equations in \(R\) and \(\Theta\). To decouple them, a spherical harmonic decomposition is necessary. This process is similar to the decoupling of the Schrödinger or Poisson equations using scalar spherical harmonics. In the context of general relativity, where Einstein’s equations are tensorial in nature, the decomposition must be approached with greater care. The components of the perturbation tensors transform as scalar, vector, or tensor spherical harmonics under rotations. Then, by choosing an appropriate coordinate gauge, it can be shown that the metric perturbation functions can be expanded as

\[\begin{split}\begin{align} \varpi^{(1)}(R,\Theta) &= \sum_{\ell} \varpi^{(1)}_{\ell}(R) \dfrac{d P_{\ell}(\cos\Theta)}{d\cos\Theta} \\[2ex] h^{(2)}(R,\Theta) &= \sum_{\ell} h^{(2)}_{\ell} P_{\ell}(\cos\Theta) \\[2ex] m^{(2)}(R,\Theta) &= \sum_{\ell} m^{(2)}_{\ell} P_{\ell}(\cos\Theta) \\[2ex] k^{(2)}(R,\Theta) &= \sum_{\ell} k^{(2)}_{\ell} P_{\ell}(\cos\Theta) \\[2ex] \xi^{(2)}(R,\Theta) &= \sum_{\ell} \xi^{(2)}_{\ell} P_{\ell}(\cos\Theta) \end{align}\end{split}\]

\[\]

where \(\varpi^{(1)}(R,\Theta) \equiv \Omega - \omega^{(1)}(R,\Theta)\).

Interior equations

Zeroth Order: The TOV equations

At order \(\mathcal{O}(1)\), the metric reduces to the background spacetime which is static and spherically symmetric. By defining the following quantity

\[\begin{equation} e^{\lambda(R)} = \left[ 1 - \dfrac{2M(R)}{R} \right]^{-1} \end{equation}\]

the \((t,t)\) and \((R, R)\) components of the Einstein’s equations yield

\[\begin{split}\begin{align} \dfrac{dM}{dR} &= 4 \pi R^{2} \varepsilon \\[3ex] \dfrac{d\nu}{dR} &= 2 \dfrac{M + 4\pi R^{3} p}{R(R-2 M)} \\[3ex] \end{align}\end{split}\]

respectively. In addition, if the equation of motion is used \(\nabla^{\mu}T_{\mu R}=0\) with the equation for \(\nu(R)\), the Tolman-Oppenheimer-Volkoff equation (TOV) is obtained,

\[\begin{equation} \dfrac{dp}{dR} = -(\varepsilon + p) \dfrac{M+4 \pi R^{3} p }{R(R-2M)} \ . \end{equation}\]

These equations along with a defined equation of state (EoS) \(p=p(\varepsilon)\) close the system of differential equations.

First Order

\[\begin{equation} \dfrac{d^{2}\varpi^{(1)}_{\ell}}{dR^{2}} + 4 \dfrac{1-\pi R^{2}(\varepsilon + p)e^{\lambda}}{R} \dfrac{d\varpi^{(1)}_{\ell}}{dR} - \left[ \dfrac{(\ell+2)(\ell-1)}{R^{2}} + 16\pi (\varepsilon + p) \right] e^{\lambda} \varpi^{(1)}_{\ell} = 0 \end{equation}\]

\[\begin{equation} \dfrac{d^{2}\varpi^{(1)}_{1}}{dR^{2}} + 4 \dfrac{1-\pi R^{2}(\varepsilon + p)e^{\lambda}}{R} \dfrac{d\varpi^{(1)}_{1}}{dR} - 16\pi (\varepsilon + p) e^{\lambda} \varpi^{(1)}_{1} = 0 \end{equation}\]

Second Order

\[\begin{split}\begin{align} \varpi^{(1)}(R,\Theta) &= \varpi^{(1)}_{1}(R) \\[2ex] h^{(2)}(R,\Theta) &= h^{(2)}_{0}(R) + h^{(2)}_{2}(R) P_{2}(\cos\Theta) \\[2ex] m^{(2)}(R,\Theta) &= m^{(2)}_{0}(R) + m^{(2)}_{2}(R) P_{2}(\cos\Theta) \\[2ex] \xi^{(2)}(R,\Theta) &= \xi^{(2)}_{0}(R) + \xi^{(2)}_{2}(R) P_{2}(\cos\Theta) \\[2ex] k^{(2)}(R,\Theta) &= k^{(2)}_{2}(R) P_{2}(\cos\Theta) \end{align}\end{split}\]

\[\]

\[\begin{split}\begin{align} \dfrac{dm^{(2)}_{0}}{dR} = & \, - 4 \pi R^{2} \dfrac{d\rho}{dR}\xi^{(2)}_{0} \, + \, S_{m^{(2)}_{0}} \\[2ex] \dfrac{d\xi^{(2)}_{0}}{dR} = & \, \dfrac{2[M(R+8\pi p R^{3}) - M^{2} - 2 \pi R^{4} (p + \rho + 8 \pi p R^{2} \rho )]}{R^{2} (M+ 4\pi p R^{3} )}e^{\lambda} \xi^{(2)}_{0} \\[2ex] & \ \ \, - \, \dfrac{1+8 \pi p R^{2}}{M + 4 \pi p R^{3} }e^{\lambda} m^{(2)}_{0} \, + \, S_{\xi^{(2)}_{0}} \end{align}\end{split}\]

\[\]

\[\begin{split}\begin{align} \dfrac{dk^{(2)}_{2}}{dR} = & - \dfrac{dh^{(2)}_{2}}{dR} + \dfrac{R-3M-4\pi p R^{3}}{R^{2}} e^{\lambda} h^{(2)}_{2} + \dfrac{R - M + 4\pi p R^{3}}{R^{3}}e^{2 \lambda} m^{(2)}_{2} \\[2ex] \dfrac{dh^{(2)}_{2}}{dR} = & - \dfrac{R-M+ 4\pi p R^{3}}{R} e^{\lambda} \dfrac{dk^{(2)}_{2}}{dR} + \dfrac{3-4\pi(\varepsilon + p)R^{2}}{R} e^{\lambda} h^{(2)}_{2} + \dfrac{2}{R} e^{\lambda} k^{(2)}_{2} \\[1ex] & + \dfrac{1+8\pi p R^{2}}{R^{2}} e^{2\lambda} m^{(2)}_{2} + \dfrac{R^{3}}{12}e^{-\nu} \left( \dfrac{d\varpi^{(1)}_{1}}{dR} \right)^{2} \end{align}\end{split}\]

\[\]

\[\begin{split}\begin{align} \xi^{(2)}_{2} &= - \dfrac{ R^2 e^{-\lambda} }{ 3(M + 4\pi p R^{3}) } \left[ 3 h^{(2)}_{2} + e^{-\nu}R^{2} \left( \varpi^{(1)}_{1} \right)^{2} \right] \\[2ex] m^{(2)}_{2} &= - Re^{-\lambda} h^{(2)}_{2} + \dfrac{1}{6}R^{4}e^{-(\nu + \lambda)} \left[ R e^{-\lambda} \left( \dfrac{d\varpi^{(1)}_{1}}{dR} \right)^{2} + 16\pi R \left( \varpi^{(1)}_{1} \right)^{2} \left( \varepsilon + p \right) \right] \\[2ex] \end{align}\end{split}\]

Tidal Deformability and Love number

The ODE for the \(\ell\)-th multipole is given by

\[\begin{split}\begin{align} &R \dfrac{dy}{dR} + y(y-1)+\dfrac{2}{f} \left[ 1 - \dfrac{3M}{R} - 2 \pi R^{2} (\varepsilon + 3p ) \right] y \\[2ex] &- \dfrac{1}{f}\left[\ell(\ell+1)-4 \pi R^{2}(\varepsilon+p)\left(3 + \dfrac{d\varepsilon}{dp} \right)\right] = 0 \end{align}\end{split}\]

Initial condition

\[\begin{equation} y(R=0)=\ell \end{equation}\]

Love number \(k_{2}\):

\[\begin{split}\begin{align} \begin{split} k_{2} =& \dfrac{8}{5}(1-2C)^{2} C^{5}\left[ 2C(Y-1) - Y + 2 \right] \\ &\times \bigg\{ 2C \bigl[ 4(Y+1) C^{4} + (6Y-4) C^{3} + (26-22 Y)C^{2} \\[1ex] & + 3(5Y-8) C- 3Y + 6 \bigr] - 3(1-2C)^{2} \bigl[ 2C(Y-1) \\ &- Y + 2 \bigr] \ln \left( \dfrac{1}{1-2C} \right) \bigg\}^{-1} \end{split} \end{align}\end{split}\]

where \(C= M/R\) is the compactness of the neutron star.

Exterior solutions

\[\begin{split}\begin{align} \varepsilon_{\mathsf{ext}}(R) &= p_{\mathsf{ext}}(R) = 0 \\[2ex] M_{\mathsf{ext}}(R) &= M_{\star} \\[2ex] \nu_{\mathsf{ext}}(R) &= -\lambda_{\mathsf{ext}}(R) = \ln \left( 1 - 2 M_{\star}/R \right) \\[2ex] \varpi^{(1)}_{1, \, \mathsf{ext}}(R) &= \Omega - \frac{2 S}{R^{3}} = \Omega \left( 1 - \frac{2 I}{R^{3}} \right) \\[2ex] m^{(2)}_{0, \, \mathsf{ext}}(R) &=-\frac{S^{2}}{R^{3}} + \delta M_{\star} \\[2ex] h^{(2)}_{0, \, \mathsf{ext}}(R) &= \dfrac{S^{2}}{R^{3}(R-2M_{\star})} - \dfrac{\delta M_{\star}}{(R-2M_{\star})} \\[2ex] h^{(2)}_{2, \, \mathsf{ext}}(R) &= \frac{1}{M_{\star} R^{3}}\left( 1+\frac{M_{\star}}{R} \right)S^{2} + A Q^{2}_{2}\left( \frac{R}{M_{\star}} - 1 \right) \\[2ex] k^{(2)}_{2, \, \mathsf{ext}}(R) &= - \frac{1}{M_{\star} R^{3}}\left( 1+\frac{2M_{\star}}{R} \right)S^{2} + \frac{2 A M_{\star}}{\sqrt{R(R-2M_{\star})}}Q^{1}_{2}\left( \frac{R}{M_{\star}} - 1 \right) - A Q^{2}_{2} \left( \frac{R}{M_{\star}} - 1 \right) \\[2ex] \end{align}\end{split}\]

Observables

Numerical Implementation

Enthalpy formulation

From the numerical perspective, since we are integrating from the center up to the surfce, we need to check at each step if the pressure already reaches zero. In practice this is done by assuming that the surface is reached when \(p=\delta p_{c}\) where \(\delta\) is a very small value (e.g., \(\delta=10^{-9}\)). In order to avoid finding the pressure at each step, the TOV equations can be reformulated in different way under a change of variables.

Let \(h=\ln\textsf{h}\) and then \(dh=\dfrac{d \textsf{h}}{\textsf{h}} = \dfrac{dp}{\varepsilon + p}\)

Important

TOV Equations in pseudo-enthalpy form

\[\begin{split}\begin{align} \begin{cases} \hspace{0.2cm} \dfrac{dR}{dh} =- \dfrac{R(R-2M)}{M + 4\pi R^{3} p} \\[3ex] \hspace{0.2cm} \dfrac{dM}{dh} =- \dfrac{4 \pi \varepsilon R^{3} (R-2M)}{M + 4 \pi R^{3} p} \\[3ex] \end{cases} \end{align}\end{split}\]

Boundary conditions

\[\begin{split}\begin{align} R(h_{\epsilon}) &= R_{\epsilon} \\[2ex] M(R_{\epsilon}) &= \dfrac{4 \pi }{3} \varepsilon_{c} R^{3}_{\epsilon} + \mathcal{O}(R^{5}_{\epsilon}) \end{align}\end{split}\]

Domain of integration: \([h_{\epsilon}, 0]\) where \(h_{\epsilon} = h_{c} - \dfrac{2 \pi}{3} R_{\epsilon}^{2}(\varepsilon_{c} + 3p_{c})\).

where \(p = p(h)\) and \(\varepsilon = \varepsilon(h)\).

\[\begin{equation} \dfrac{d\nu}{dh} = -2 \end{equation}\]

\[\begin{equation} \nu(h) = - 2h + \ln \left( 1-\dfrac{2M_{\star}}{R_{\star}} \right) \end{equation}\]

Important

\[\begin{split}\begin{align} \begin{cases} \hspace{0.2cm} \dfrac{d \varpi_{1}^{(1)}}{dR} = \alpha_{1}^{(1)} \\[2ex] \hspace{0.2cm} \dfrac{d\alpha^{(1)}_{1}}{dR} + 4 \dfrac{1-\pi R^{2}(\varepsilon + p)e^{\lambda}}{R} \dfrac{d\varpi^{(1)}_{1}}{dR} - 16\pi (\varepsilon + p) e^{\lambda} \varpi^{(1)}_{1} = 0 \end{cases} \end{align}\end{split}\]

\[\begin{split}\varpi^{(1)}_{1}(R_{\epsilon}) &= \varpi^{(1)}_{1, c} + \dfrac{8 \pi}{5} (\varepsilon_{c} + p_{c}) \varpi^{(1)}_{1,c} R_{\epsilon}^{2} + \mathcal{O}(R_{\epsilon}^{4}) \\[2ex] \alpha^{(1)}_{1}(R_{\epsilon}) &= \dfrac{16 \pi}{5} (\varepsilon_{c} + p_{c}) \varpi^{(1)}_{1,c} R_{\epsilon}^{2} + \mathcal{O}(R_{\epsilon}^{4})\end{split}\]

Important

\[\begin{equation} \begin{cases} R \dfrac{dy}{dR} + y^{2} + y e^{\lambda} \left[ 1 + 4 \pi R^{2} \left( p - \varepsilon \right) \right] + R^{2} Q = 0 \end{cases} \end{equation}\]

where

\[\begin{equation} Q = 4 \pi e^{\lambda} \left( 5 \varepsilon + 9 p + \dfrac{\varepsilon + p}{c_{s}^{2} } \right) - 6\dfrac{e^{\lambda}}{R^{2}} - \left( \dfrac{d\nu}{dR} \right)^{2} \end{equation}\]

Boundary condition: \(y(R_{\epsilon}) = 2\)

Asymptotic Analysis

Since these equations diverge exactly at \(R=0\) we need to generate boundary conditions at a very small core around the center of the star. Thus, we impose regularity at the center and Taylor-expand the unknown functions around \(R=0\). By plugging the power series expansion into the differential equations we obtain an algebraic system of equations for the coefficients of the series expansion for every unknown function. The final result is,

\[\begin{split}\begin{align} R(h_{\epsilon}) &= R_{\epsilon} \\[2ex] M(R_{\epsilon}) &= \dfrac{4 \pi }{3} \varepsilon_{c} R^{3}_{\epsilon} + \mathcal{O}(R^{5}_{\epsilon}) \\[2ex] \nu(R_{\epsilon}) &= \nu_{c} + \dfrac{4 \pi}{3}(\varepsilon_{c} + 3p_{c})R^{2} + \mathcal{O}(R^{4}) \\[2ex] p(R_{\epsilon}) &= p_{c} + \dfrac{2 \pi}{3}(\varepsilon_{c} + p_{c})(\varepsilon_{c} + 3 p_{c})R_{\epsilon}^{2} + \mathcal{O}(R^{4}_{\epsilon}) \\[2ex] \varepsilon(R_{\epsilon}) &= \varepsilon_{c} + \dfrac{2 \pi}{3}\dfrac{(\varepsilon_{c} + p_{c})(\varepsilon_{c} + 3 p_{c})R_{\epsilon}^{2}}{c^{2}_{s,c}} + \mathcal{O}(R^{4}_{\epsilon}) \\[2ex] \varpi^{(1)}_{1}(R_{\epsilon}) &= \varpi^{(1)}_{1, c} + \dfrac{8 \pi}{5} (\varepsilon_{c} + p_{c}) \varpi^{(1)}_{1,c} R_{\epsilon}^{2} + \mathcal{O}(R_{\epsilon}^{4}) \\[2ex] \varpi'^{(1)}_{1}(R_{\epsilon}) &= \dfrac{16 \pi}{5} (\varepsilon_{c} + p_{c}) \varpi^{(1)}_{1,c} R_{\epsilon}^{2} + \mathcal{O}(R_{\epsilon}^{4}) \end{align}\end{split}\]

with \(h_{\epsilon} = h_{c} - \dfrac{2 \pi}{3} R_{\epsilon}^{2}(\varepsilon_{c} + 3p_{c})\).

\[\begin{align} \end{align}\]

Numerical methods

• Interpolation: Stephan’s method

• ODE GSL solver: RKF45

References

[1] Hartle, J. B. (1967). Slowly rotating relativistic stars. I. Equations of structure. The Astrophysical Journal, 150, 1005.

[2] Lindblom, L. (1992). Determining the nuclear equation of state from neutron-star masses and radii. The Astrophysical Journal, 398, 569-573.

[3] Hartle, J. B., & Thorne, K. S. (1968). Slowly rotating relativistic stars. II. Models for neutron stars and supermassive stars. The Astrophysical Journal, 153, 807.

[4] Yagi, K., & Yunes, N. (2013). I-Love-Q relations in neutron stars and their applications to astrophysics, gravitational waves, and fundamental physics. Physical Review D, 88(2), 023009.

[5] Binnington, T., & Poisson, E. (2009). Relativistic theory of tidal Love numbers. Physical Review D, 80(8), 084018.

[6] Hinderer, T. (2008). Tidal Love numbers of neutron stars. The Astrophysical Journal, 677(2), 1216.