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 and matter distribution: 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. 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 assumed to be negligible.

2. 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\).

3. Uniform rotation: Although differential rotation is particularly relevant in newly formed neutron stars, it has been shown that normal stars with differential rotation tend to evolve toward a uniform equilibrium state [1]. The method originally assumes uniform rotation for simplicity, although it could be extended to include differential rotation.

4. Stationary and axial 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. Thus, the metric components must be independent of the time coordinate \(t\) and the azimuthal angle \(\phi\).

5. Reflection symmetry: The geometry remains unchanged when mirrored across a plane perpendicular to the axis of rotation. This implies that the geometry should be invariant under the transformation \(\theta\) \(\rightarrow\) \(\pi-\theta\).

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 [2]

\[\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-frequency 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 four-velocity of the fluid, 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] u_{\mu}(r,\theta) &= u_{\mu}^{(0)}(r) + \epsilon u_{\mu}^{(1)}(r,\theta) + \frac{1}{2} \epsilon^{2} u_{\mu}^{(2)} (r,\theta) + \cdots \end{align}\end{split}\]

\[\]

The fluid’s four velocity is defined as

\[\begin{align} u^{\mu} = (u^{t}, 0, 0, \epsilon \Omega u^{t}) \end{align}\]

\[\]

where the component \(u^{t}\) can be found from the normalization condition \(u_{\alpha}u^{\alpha}=-1\). Therefore, the stress-energy tensor can be expanded as

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

\[\]

The metric perturbation function depends on both the radial coordinate, \(r\), and the angular coordinate, \(\theta\). To decouple Einstein’s equations, a harmonic decomposition may be used in the Regge-Wheeler gauge (see [2] for details)

\[\begin{split}\begin{align} \label{eq:hexp} h^{(2)}(r,\theta) &= \sum_{\substack{\ell = 0 \, \\[0.2ex] ( \ell \textsf{ even} ) }}^{2} h^{(2)}_{\ell}(r) P_{\ell}(\cos \theta) \, , \\[1ex] \label{eq:kexp} k^{(2)}(r,\theta) &= \sum_{\substack{\ell = 2 \, \\[0.2ex] ( \ell \textsf{ even} ) }}^{2} k^{(2)}_{\ell}(r) P_{\ell}(\cos \theta) \, , \\[1ex] \label{eq:mexp} m^{(2)}(r,\theta) &= \sum_{\substack{\ell = 0 \\[0.2ex] ( \ell \textsf{ even} ) }}^{2} m^{(2)}_{\ell}(r) P_{\ell}(\cos \theta) \, , \\[1ex] \label{eq:wexp} \omega^{(1)}(r,\theta) &= \sum_{\substack{\ell = 1 \\[0.2ex] ( \ell \textsf{odd} ) }}^{k} \omega^{(1)}_{\ell}(r) \dfrac{dP_{\ell}(\cos \theta)}{d\cos\theta} \, , \end{align}\end{split}\]

\[\]

and similar expressions can be obtained for the perturbations of the pressure and energy density. 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. Since the function \(\xi^{(2)}\) behaves as a scalar under rotations it can be decomposed as,

\[\begin{split}\begin{align} \xi^{(s)}(R,\Theta) = \sum_{\substack{\ell = 0 \\[0.2ex] ( \ell \textsf{ even} ) }}^{2} \xi^{(2)}_{\ell}(R) P_{\ell}(\cos \Theta) \, . \end{align}\end{split}\]

\[\]

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

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

Up to first order in the star’s spin frequency, there is a single equation for the \((t, \phi)\) component of the perturbed Einstein’s equations, given by:

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

where

\[\begin{equation} \varpi^{(1)}_{1}(R,\Theta) \equiv \Omega - \omega^{(1)}_{1}(R,\Theta) \end{equation}\]

is the coordinate angular velocity of the fluid element at \((R, \Theta)\) as seen by a free-falling observer, and \(\Omega\) is the angular velocity of the fluid as observed by an observer at rest at some point within the fluid.

Second Order

Up to second order in the spin-rotation expansion, there are two systems of equations associated with the \(\ell=0\) and \(\ell=2\) modes. The first system describes the correction to the TOV mass due to rotation, while the second system allows us to obtain the quadrupole deformation of the star induced by rotation. These system of equations for each mode are given by

\(\ell=0\)

\[\begin{split}\begin{align} \label{eq:m0} \dfrac{dm^{(2)}_{0}}{dR} &= \, - 4 \pi R^{2} \dfrac{d\varepsilon}{dR}\xi^{(2)}_{0} \, + \, S_{m^{(2)}_{0}} \ , \\[1ex] \label{eq:h0} \dfrac{dh^{(2)}_{0}}{dR} &= \, \dfrac{1 + 8 \pi p R^{2}}{R^{2}}e^{2 \lambda} m^{(2)}_{0} \, + \, \dfrac{4 \pi R (M + 4 \pi p R^{3})(p +\varepsilon)}{R^{2}}e^{2\lambda} \xi^{(2)}_{0} \, + \, S_{h^{(2)}_{0}} \, , \\[1ex] \dfrac{d\xi^{(2)}_{0}}{dR} &= \, \dfrac{2[M(R+8\pi p R^{3}) - M^{2} - 2 \pi R^{4} (p + \varepsilon + 8 \pi p R^{2} \varepsilon )]}{R^{2} (M+ 4\pi p R^{3} )}e^{\lambda} \xi^{(2)}_{0} \\ \, &- \, \dfrac{1+8 \pi p R^{2}}{M + 4 \pi p R^{3} }e^{\lambda} m^{(2)}_{0} \, + \, S_{\xi^{(2)}_{0}} \, , \label{eq:xi0} \end{align}\end{split}\]

\[\]

where

\[\begin{split}\begin{align} S_{m^{(2)}_{0}} &= \frac{1}{12} R^3 e^{-\nu (R)} \left\{ \left( R-2 M \right) \left(\dfrac{d\varpi^{(1)}_{1}}{dR} \right)^{2} +32 \pi R \left( \varpi^{(1)}_{1} \right)^{2} (p +\varepsilon)\right\} \, ,\\ S_{h^{(2)}_{0}} &=-\dfrac{1}{12} R^{3} e^{-\nu(R)} \left( \dfrac{d\varpi^{(1)}_{1}}{dR}\right)^{2} \, ,\\ S_{\xi^{(2)}_{0}} &= \dfrac{R^2 e^{-\nu} }{12 \left(M +4 \pi R^3 p\right)} \bigg\{ R^{2} (R-2 M) \left( \dfrac{d\varpi^{(1)}_{1}}{dR}\right)^{2}\\ & -8 \left( \varpi^{(1)}_{1} \right)^{2} \left(3 M + 4 \pi R^3 p - R \right) + 8 R (R-2 M) \varpi^{(1)}_{1} \dfrac{d\varpi^{(1)}_{1}}{dR} \bigg\} \, . \end{align}\end{split}\]

\[\]

\(\ell=2\)

\[\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 OF NON-ROTATING STARS

Another important macroscopic observable of non-rotating neutron stars is the tidal deformability, which measures how easily a star can be deformed under the influence of a static external field. Although tidal deformations have been formally explored in ([6] , [7]), an interesting approach to deriving the equation governing tidal deformation is to use the Hartle-Thorne expansion up to second order by setting the function \(\varpi^{(1)}_{1}=0\) in the differential equation for \(h^{(2)}_{2}\) [5]. This leads to an equation of the following form,

\[\begin{split}\begin{align} &\dfrac{d^{2}h^{(2)}_{2}}{dR^{2}} + \left\{ \dfrac{2}{R} + \left[ \dfrac{2M}{R^{2}} + 4 \pi R(p - \varepsilon) \right] e^{\lambda} \right\} \dfrac{dh^{(2)}_{2}}{dR} \\ &- \left\{ \dfrac{ 6e^{\lambda} }{R^{2}} - 4\pi \left[ 5\varepsilon + 9p + (\varepsilon + p ) \dfrac{d\varepsilon}{dp} \right] e^{\lambda} + \left( \dfrac{d\nu}{dR} \right)^{2} \right\}h^{(2)}_{2} = 0 \, . \end{align}\end{split}\]

From the previous equation, the tidal Love number, or tidal deformability of the star can be extracted. For more details, see the subsection on observables extraction.

Asymptotic Analysis

Since the interior equations at each spin-frequency order diverge at \(R = 0\), it is necessary to generate boundary conditions within a very small core around the center of the star. Regularity is imposed at the center, and the unknown functions are Taylor-expanded around \(R = 0\). By substituting the power series expansion into the differential equations, an algebraic system is obtained for the coefficients of the series expansion of each unknown function. The asymptotic expansions used for generating the boundary conditions are presented in subsection Numerical Implementation.

Observables extraction

The process of extracting the observables or multipole moments of the star involves solving the interior equations and matching them to the exterior solutions at each order. This procedure must be carried out order by order in the slow-rotation expansion to extract each macroscopic observable.

Zeroth Order

At zeroth order, since the spacetime is static and spherically symmetric, the exterior solution is the Schwarzschild solution. Thus,

\[\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] \end{align}\end{split}\]

If we take the \((t,t)\) component of the exterior solution and expand it at spatial infinity, we obtain \(g^{(0)}_{tt} = -1 + \frac{2M_{\star}}{R} + \cdots\). Thus, the quantity \(M_{*}\) is interpreted as the TOV mass monopole of the source distribution of a non-rotating star. The TOV radius, \(R_{\star}\), of the star is obtained by matching the pressure at the surface, which implies finding the radial coordinate, \(R\), at which \(p(R=R_{\star}) = 0\).

\[\]

First Order

At linear order, the spin angular momentum of the slowly rotating star can be extracted by matching the interior and exterior solution at the neutron star surface \(R_{\star}\) and also imposing continuity of the function \(\omega^{(1)}_{1}(R)\). This implies, that the matching conditions are given by,

\[\begin{split}\begin{align} \varpi^{(1)}_{1, \textsf{int}}(R_{\star}) &= \varpi^{(1)}_{1, \textsf{ext}}(R_{\star}) \\ \dfrac{d\varpi^{(1)}_{1, \textsf{int}}(R_{\star})}{dR} &= \dfrac{d\varpi^{(1)}_{1, \textsf{ext}}(R_{\star})}{dR} \end{align}\end{split}\]

where the exterior solution gives,

\[\begin{split}\begin{align} \varpi^{(1)}_{1, \, \mathsf{ext}}(R) &= \Omega - \frac{2 S}{R^{3}} = \Omega \left( 1 - \frac{2 I}{R^{3}} \right) \\[2ex] \end{align}\end{split}\]

with \(\Omega\) being the angular speed of the neutron star, \(S\) it’s spin angular momentum and \(I\) the moment of inertia. Since there are two equations for two unknowns, \((\Omega, S)\) then they are giving by

\[\begin{align} S = \dfrac{1}{6} R_{\star}^{4} \left( \dfrac{d \varpi^{(1)}_{1, \textsf{int}}}{dR} \right)_{R=R_{\star}} \hspace{0.5cm} ; \hspace{0.5cm} \Omega = \varpi^{(1)}_{1, \textsf{int}}(R_{\star}) + \dfrac{2S}{R_{\star}^{3}} \, . \end{align}\]

Once \(S\) and \(\Omega\) are obtained, the moment of inertia is defined as \(I = S/\Omega\). This module computes the dimensionless moment of inertia which is defined as

\[\begin{align} \bar{I} \equiv \dfrac{I}{M_{\star}^{3}} \, . \end{align}\]

As mentioned by Hartle [4], if a different value of \(\Omega\) is desired (e.g., another value interested by the user), one rescales the function \(\varpi^{(1)}_{1}\) as

\[\begin{align} \varpi^{(1)}_{1, \textsf{user}} = \dfrac{\Omega_{\textsf{user}}}{\Omega} \varpi^{(1)}_{1} \, . \end{align}\]

Second Order

At second order, we can extract the first correction to the mass monopole, denoted as \(\delta M_{*}\), the equatorial radius, \(R_{\textrm{eq}}\), the first correction to the mean radius of the star, \(\langle \delta R_{\star} \rangle\), the stellar eccentricity, \(e_{s}\), and the rotational quadrupole moment, \(Q^{\textrm{(rot)}}\). The procedure is similar to the previous one, where the interior and exterior solutions must be matched at the surface of the neutron star.

\(\ell=0\)

Mass correction

To extract the first correction to the TOV mass \(\delta M_{\star}\), we solve the \(\ell=0\) equations and math the interior and exterior solutions at the surface of the neutron star.

\[\begin{align} m^{(2)}_{0, \, \mathsf{int}}(R) &=m^{(2)}_{0, \, \mathsf{ext}}(R) \end{align}\]

where

\[\begin{align} m^{(2)}_{0, \, \mathsf{ext}}(R) &=-\frac{S^{2}}{R^{3}} + \delta M_{\star} \, . \end{align}\]

From the matching condition the quantity \(\delta M_{\star}\) is obtained as

\[\begin{align} \delta M_{\star} = m^{(2)}_{0, \textsf{int}}(R_{\star}) + \dfrac{S}{R_{\star}^{3}} \, . \end{align}\]

In principle the quantity \(\delta M_{\star}\) is a integration constant. However, if we take the exterior metric solution and expand the \((t,t)\) components at spatial infinity it is obtained

\[\begin{align} g_{tt} = -1 + \dfrac{2\left( M_{\star} + \delta M_{\star} \right)}{R} + \cdots \end{align}\]

and therefore, the mass monopole of the distribution is given by

\[\begin{align} M = M_{\star} + \delta M_{\star} \, . \end{align}\]

The quantity \(\delta M_{\star}\) depends in general on the angular spin-frequency of the star \(\Omega\). To obtain the same quantity for another value of \(\Omega\) specified by the user, \(\Omega_{\textsf{user}}\). Then,

\[\begin{align} \delta M_{\star, \textsf{user}} = \delta M_{\star} \dfrac{\Omega^{2}_{\textsf{user}}}{\Omega^{2}} \, . \end{align}\]

Mean radius correction

When the star is rotating, the surface of the neutron star gets deformed. To find an invariant parametrization of the stellar surface, we look for the surface in 3D space with spherical coordinates \((\tilde{r}, \theta, \phi)\) with the same intrinsic geometry as the surface of constant density of the star at \(R=R_{\star}\). The surface contour is then described by

\[\begin{align} \tilde{r} = R_{\star} + \epsilon^{2}\left[ \xi^{(2)}_{0} + \left( R_{\star} k^{(2)}_{2} + \xi^{(2)}_{2} \right) P_{2}(\cos\theta) \right]_{R_{\star}} \end{align}\]

By taking an integral average over the 2-sphere then the total mean radius of the star is obtained,

\[\begin{align} \langle \tilde{r} \rangle = R_{*} + \epsilon^{2} \langle \delta R \rangle \end{align}\]

where \(\langle \delta R \rangle \equiv \xi^{(2)}_{0}(R_{\star})\). We refer to the quantity \(\langle \delta R \rangle\) as the first mean contribution the TOV radius \(R_{\star}\). Similarly, the quantity \(\langle \delta R \rangle\) can be rescaled as

\[\begin{align} \langle \delta R \rangle_{\textsf{user}} = \langle \delta R \rangle \dfrac{\Omega^{2}_{\textsf{user}}}{\Omega^{2}} \, . \end{align}\]

\(\ell=2\)

Rotational quadrupole moment

From the \(\ell = 0\) set of equations, the rotational quadrupole moment of the star, denoted as \(Q^{\textsf{(rot)}}\), can be extracted. First, the two metric perturbation functions at second order in the slow-rotation expansion, \(h^{(2)}_{2}(R)\) and \(k^{(2)}_{2}(R)\), are solved numerically inside the star. Then, by matching the internal solution to the external one, it is possible to extract two integration constants: one inside the star and the other outside. In particular, the integration constant outside the star is the one required to extract the quadrupole moment. To see this, recall that the general solution of an ODE consists of the particular solution plus the homogeneous solution. Thus,

\[\begin{split}\begin{align} h^{\, (2)}_{2, \, \mathsf{int}}(R) &= h^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{int}}(R) \, + \, C^{\, (2)}_{\ell, \, \mathsf{int}} \, h^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{int}}(R) \, , \\ h^{\, (2)}_{2, \, \mathsf{ext}}(R) &= h^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{ext}}(R) \, + \, C^{\, (2)}_{\ell, \, \mathsf{ext}} \, h^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{ext}}(R) \\ \end{align}\end{split}\]

and similarly for \(k^{(2)}_{2}(R)\)

\[\begin{split}\begin{align} k^{\, (2)}_{2, \, \mathsf{int}}(R) &= k^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{int}}(R) \, + \, C^{\, (2)}_{\ell, \, \mathsf{int}} \, k^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{int}}(R) \, , \\ k^{\, (2)}_{2, \, \mathsf{ext}}(R) &= k^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{ext}}(R) \, + \, C^{\, (2)}_{\ell, \, \mathsf{ext}} \, k^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{ext}}(R) \\ \, . \end{align}\end{split}\]

Then, by matching both functions, \(h^{(2)}_{2}(R)\) and \(k^{(2)}_{2}(R)\) at the neutron star boundary, the integration constants can be found as,

\[\begin{split}\begin{align} C^{\, (2)}_{2, \, \mathsf{int}} =& \dfrac{ \left[ h^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{ext}}(R_{*}) - h^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{int}}(R_{*})\right] k^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{ext}}(R_{*}) + h^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{ext}}(R_{*}) \left[ k^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{int}}(R_{*}) - k^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{ext}}(R_{*}) \right] }{ h^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{int}}(R_{*}) k^{\, (2) \, \mathsf{H}}_{\ell, \, \mathsf{ext}}(R_{*}) - h^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{ext}}(R_{*}) k^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{int}}(R_{*}) } \, , \\[1ex] C^{\, (2)}_{2, \, \mathsf{ext}} =& \dfrac{ \left[ h^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{ext}}(R_{*}) - h^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{int}}(R_{*})\right] k^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{int}}(R_{*}) + h^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{int}}(R_{*}) \left[ k^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{int}}(R_{*}) - k^{\, (2) \, \mathsf{P}}_{2, \, \mathsf{ext}}(R_{*}) \right] }{ h^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{int}}(R_{*}) k^{\, (2) \, \mathsf{H}}_{\ell, \, \mathsf{ext}}(R_{*}) - h^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{ext}}(R_{*}) k^{\, (2) \, \mathsf{H}}_{2, \, \mathsf{int}}(R_{*}) } \, . \end{align}\end{split}\]

The homogeneous and particular exterior solutions can be read out from

\[\begin{split}\begin{align} 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}\]

where \(Q^{m}_{\ell}\) is the associated Legendre function of the second kind,

\[\begin{split}\begin{align} Q^{1}_{1}(\zeta) &= (\zeta^{2}-1)^{1/2} \left[ \dfrac{3\zeta^{2}-2}{\zeta^{2}-1} -\dfrac{3}{2}\log \left( \dfrac{\zeta +1}{\zeta -1} \right) \right] \, , \\[1ex] Q^{2}_{2}(\zeta) &= \dfrac{3}{2}(\zeta^{2}-1) \log \left( \dfrac{\zeta + 1}{\zeta -1} - \dfrac{3\zeta^{3}-5\zeta}{\zeta^{2}-1} \right) \, . \end{align}\end{split}\]

Using the analytical exterior solution and expanding the \(g_{tt}\) component of the metric at spatial infinity it is obtained:

\[\begin{align} g_{tt}(R) = - 1 + \frac{2 \left( M_{\star } + \delta M_{\star} \right)}{R} - \dfrac{ 2 \left( 8 M_{\star}^4 C^{(2)}_{2,\mathsf{ext}} + 5 S^{2} \right) }{5 R^{3} M_{\star}} P_{2}(\cos\Theta) + \cdots \end{align}\]

The quadrupole moment is the term proportional to \(2P_{2}/R^{3}\). Therefore,

\[\begin{align} Q^{(\textrm{rot})} = -\frac{S^2}{M_{\star }} -\frac{8}{5} M_{\star }^3 C^{(2)}_{2, \mathsf{ext}} \, . \end{align}\]

As before, such quantity depends in general on the speed spin-frequency of the star. Then, this quantity can be rescaled as

\[\begin{align} Q^{\mathrm{(rot)}}_{\mathsf{user}} = Q^{\mathrm{(rot)}} \dfrac{\Omega^{2}_{\textsf{user}}}{\Omega^{2}} \, . \end{align}\]

Equatorial radius and stellar eccentricity

To compute the equatorial radius correction due to the rotation of the neutron star, we can use the contour invariant parametrization of the surface. The equatorial and polar radius of the star are defined at \(\Theta=\pi/2\) and \(\Theta=0\), respectively. Thus,

\[\begin{split}\begin{align} \label{eq:Req} R_{\textrm{eq}} &= R_{*} + \delta R_{\textrm{eq}} \, , \\ R_{\textrm{pol}} &= R_{*} + \delta R_{\textrm{pol}} \, . \end{align}\end{split}\]

where

\[\begin{split}\begin{align} \delta R_{\textrm{eq}} &= \left[ \xi^{(2)}_{0} -\dfrac{1}{2} \left( R_{*} k^{(2)}_{2} + \xi^{(2)}_{2} \right) \right]_{R_{*}} \, , \\ \delta R_{\textrm{pol}} &= \left[ \xi^{(2)}_{0} + R_{*} k^{(2)}_{2} + \xi^{(2)}_{2} \right]_{R_{*}} \, . \end{align}\end{split}\]

The stellar eccentricity can be defined from the previous expressions as:

\[\begin{split}\begin{align} e_{s} &\equiv \left[ \left( R_{\textrm{eq}}/R_{\textrm{pol}} \right)^{2} -1 \right]^{1/2} \\ &= \left[ -\dfrac{3 \epsilon^{2}}{R_{*}}\left( R_{*} k^{(2)}_{2} + \xi^{(2)}_{2} \right) \right]^{1/2}_{R_{*}} \, . \label{eq:ecc} \end{align}\end{split}\]

Notice that the quantities \(\delta R_{eq}\), \(\delta R_{pol}\) and \(e_{s}\) depend on the angular spin-frequency of the star. To rescale the quantities according to a desired speed angular frequency then,

\[\begin{align} e^{\mathsf{user}}_{s} = \dfrac{e_{s}}{\Omega} \hspace{0.5cm} ; \hspace{0.5cm} \left( \delta R_{\mathrm{eq}} \right)_{\textsf{user}} = \delta R_{\mathrm{eq}} \dfrac{\Omega^{2}_{\textsf{user}}}{\Omega^{2}} \hspace{0.5cm} ; \hspace{0.5cm} \left( \delta R_{\mathrm{pol}} \right)_{\textsf{user}} = \delta R_{\mathrm{pol}} \dfrac{\Omega^{2}_{\textsf{user}}}{\Omega^{2}} \, . \end{align}\]

Tidal deformability

The exterior solution to the equation presented earlier for the tidal deformation is given by

\[\begin{split}\begin{align} h^{(2), \mathrm{tid}}_{2,\mathsf{ext}} &= c_{1} \left( \dfrac{R}{M_{\star}} \right) \left( 1 - \dfrac{2 M_{\star}}{R} \right) \bigg[ - \dfrac{2M_{\star} \left( R-M_{\star} \right) \left( 3R^{2} - 6M_{\star} R -2 M_{\star}^{2} \right)}{R^{2}(R-2M_{\star})^{2}} \\ &+ 3 \ln\left( \dfrac{R}{R-2M_{\star}} \right) \bigg] + c_{2} \left( \dfrac{R}{M_{\star}} \right)^{2} \left( 1 - \dfrac{2M_{\star}}{R} \right) \end{align}\end{split}\]

Expanding this expression in the buffer zone we obtain [5]

\[\begin{align} h^{(2), \mathrm{tid}}_{2,\mathsf{ext}} = \dfrac{16}{5} c_{1}\dfrac{M_{\star}^{3}}{R^{3}} + c_{2} \dfrac{R^{2}}{M_{\star}^{2}} + \cdots \end{align}\]

By comparing the general expansion of the metric of a tidally distorted object of mass \(M_{\star}\) in ACMC coordinates, it can be seen that \(c_{1}\) is related to the tidal quadrupole moment \(Q^{\mathrm{(tid)}}\) and the constant \(c_{2}\) to the tidal field \(\mathcal{E}^{\mathsf{(tid)}}\). Then the apsidal constant \(k_{2}\) which is related to the tidal deformability is given by \((8/5)C^{5}(c_{1}/c_{2})\). By matching the internal and exterior solution at the neutron star surface,

\[\begin{split}\begin{align} h^{(2),\textrm{tid}}_{2, \textrm{int}}(R_{\star}) &= h^{(2),\textrm{tid}}_{2, \textrm{ext}}(R_{\star}) \\ h'^{(2),\textrm{tid}}_{2, \textrm{int}}(R_{\star}) &= h'^{(2),\textrm{tid}}_{2, \textrm{ext}}(R_{\star}) \\ \end{align}\end{split}\]

we can express \(c_{1}\) and \(c_{2}\) in terms of \(h^{(2),\textrm{tid}}_{2, \textrm{int}}(R_{\star})\) and it’s derivative \(h'^{(2),\textrm{tid}}_{2, \textrm{int}}(R_{\star})\). Thus, the apsidal constant gives,

\[\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_{\star}/R_{\star}\) is the compactness of the neutron star and

\[\begin{align} Y\equiv y(R_{\star}) = R_{\star} \dfrac{h'^{(2),\textrm{tid}}_{2, \textrm{int}}(R_{\star})}{ h^{(2),\textrm{tid}}_{2, \textrm{int}}(R_{\star})} \, . \end{align}\]

Once the apsidal constant is computed, the dimensionless static tidal Love number or tidal deformability is given by

\[\begin{align} \bar{\lambda}^{(\mathrm{tid})} = \dfrac{2}{3} k_{2} C^{-5}. \end{align}\]

Numerical Implementation

Zeroth Order

Enthalpy formulation

From a numerical standpoint, as we integrate from the center to the surface, it is essential to check at each step whether the pressure has already reached 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}\)). To avoid computing the pressure at each step, the TOV equations can be reformulated by applying a change of variables [3].

Important

TOV Equations in pseudo-enthalpy form

Let \(h=\ln\textsf{h}\) and then \(dh=\dfrac{d \textsf{h}}{\textsf{h}} = \dfrac{dp}{\varepsilon + p}\). Then, the TOV equation can be written as

\[\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)\). The extra equation for the metric function \(\nu\) gives,

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

and the solution is just given by

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

First Order

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

Second Order

Important

• \(\ell=0\) equations:

\[\begin{split}\begin{align} \dfrac{dm^{(2)}_{0}}{dR} &= \, - 4 \pi R^{2} \dfrac{d\varepsilon}{dR}\xi^{(2)}_{0} \, + \, S_{m^{(2)}_{0}} \ , \\[1ex] \dfrac{d\xi^{(2)}_{0}}{dR} &= \, \dfrac{2[M(R+8\pi p R^{3}) - M^{2} - 2 \pi R^{4} (p + \varepsilon + 8 \pi p R^{2} \varepsilon )]}{R^{2} (M+ 4\pi p R^{3} )}e^{\lambda} \xi^{(2)}_{0} \\ \, &- \, \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}\]

\[\]

where

\[\begin{split}\begin{align} S_{m^{(2)}_{0}} &= \frac{1}{12} R^3 e^{-\nu (R)} \left\{ \left( R-2 M \right) \left(\dfrac{d\varpi^{(1)}_{1}}{dR} \right)^{2} +32 \pi R \left( \varpi^{(1)}_{1} \right)^{2} (p +\varepsilon)\right\} \, ,\\ S_{\xi^{(2)}_{0}} &= \dfrac{R^2 e^{-\nu} }{12 \left(M +4 \pi R^3 p\right)} \bigg\{ R^{2} (R-2 M) \left( \dfrac{d\varpi^{(1)}_{1}}{dR}\right)^{2}\\ & -8 \left( \varpi^{(1)}_{1} \right)^{2} \left(3 M + 4 \pi R^3 p - R \right) + 8 R (R-2 M) \varpi^{(1)}_{1} \dfrac{d\varpi^{(1)}_{1}}{dR} \bigg\} \, . \end{align}\end{split}\]

\[\]

Boundary conditions

\[\begin{split}\begin{align} m^{(2)}_{0}(R) &= \frac{4 \pi e^{-\nu_c}}{15 p'_{c}} \left(p_c+\varepsilon _c\right) \left(1 + 2 p'_{c}\right)\left( \varpi^{(1)}_{1,c} \right)^{2}R^{5} + \mathcal{O}(R^{7}) \\[1ex] \xi^{(2)}_{0}(R) &= \frac{ e^{-\nu_c} }{4 \pi \left( 3 p_c+ \varepsilon_c \right)}\left( \varpi^{(1)}_{1,c}\right)^{2}R + \mathcal{O}(R^{3}) \end{align}\end{split}\]

• \(\ell=2\) equations:

Inhomogeneous system

\[\begin{split}\begin{align} \dfrac{dh^{(2)}_{2}}{dR} &= a_{1} h^{(2)}_{2} + b_{1} k^{(2)}_{2} + \tilde{S}_{h^{(2)}_{2}} \\ \dfrac{dk^{(2)}_{2}}{dR} &= a_{2} k^{(2)}_{2} + b_{2} h^{(2)}_{2} + \tilde{S}_{k^{(2)}_{2}} \\ \end{align}\end{split}\]

Homogeneous system

\[\begin{split}\begin{align} \dfrac{dh^{(2)}_{2}}{dR} &= a_{1} h^{(2)}_{2} + b_{1} k^{(2)}_{2} \\ \dfrac{dk^{(2)}_{2}}{dR} &= a_{2} k^{(2)}_{2} + b_{2} h^{(2)}_{2} \\ \end{align}\end{split}\]

where \(a_{1}\), \(b_{1}\), \(a_{2}\), and \(b_{2}\) are radial functions given by

\[\begin{split}\begin{align} a_{1} =& \frac{2 \left\{ M \left[ R-4 \pi R^3 (3 p + \varepsilon)\right] + M^2+2 \pi R^4 \left(-8 \pi R^2 p^2+p+\varepsilon\right)-R^2\right\}}{R (R-2 M) \left(M+4 \pi R^3 p\right)} \\[1ex] b_{1} =& -\frac{2}{M+4 \pi R^3 p} \\[1ex] a_{2} =& \frac{2}{M(R)+4 \pi R^3 p(R)} \\[1ex] b_{2} =& \frac{2 \left(M-2 \pi R^3 (p+\varepsilon)+R\right)}{R \left(M+4 \pi R^3 p \right)} \end{align}\end{split}\]

and the source terms are

\[\begin{split}\begin{align} S_{k^{(2)}_{2}} &= \dfrac{R^2 e^{-\nu}}{12 (R-2 M) \left(M + 4 \pi R^3 p \right)} \Bigg\{ 512 \pi ^3 R^7 p^3 \left( \varpi^{(1)}_{1} \right)^{2} \\ & - \left(-4 R^2 M + 2 R M^{2} + 4 M^{3} + R^{3} \right) \left( \dfrac{d\varpi^{(1)}_{1}}{dR} \right)^{2} \\ & + 16 \varepsilon \pi R \left(-2 R M +2 M^2+R^2\right) \left( \varpi^{(1)}_{1} \right)^{2} + 16 \pi R p \Big[ R^2 M (R-2 M) \left( \dfrac{d\varpi^{(1)}_{1}}{dR} \right)^{2} \\ & + \left( \varpi^{(1)}_{1}\right)^{2} \left(2 R M \left(8 \pi R^2 \varepsilon-1\right)+2 M^2+R^2\right)\Big] \\ & + 32 \pi ^2 R^4 p^2 \Bigg[ M \left(8 \left( \varpi^{(1)}_{1}\right)^2-2 R^2 \left( \dfrac{d\varpi^{(1)}_{1}}{dR} \right)^{2} \right) \\ &+R^3 \left(\left( \dfrac{d\varpi^{(1)}_{1}}{dR} \right)^{2} +16 \pi \varepsilon \left( \varpi^{(1)}_{1}\right)^2 \right) \Bigg] \Bigg\} \, , \\[2ex] S_{h^{(2)}_{2}} &= \dfrac{e^{-\nu}R^{2}}{12 \left(M+4 \pi R^3 p\right)} \Bigg\{ 128 \pi ^2 R^4 p^2 \left( \varpi^{(1)}_{1} \right)^2 \\ & + 8 \pi R p \left[ 2 \left( \varpi^{(1)}_{1} \right)^{2} \left(2 M + 8 \pi R^3 \varepsilon - R\right) \\+ R^2 (R-2 M) \left( \dfrac{d\varpi^{(1)}_{1}}{dR} \right)^{2}\right] \\ & -(R-2 M) \left(16 \pi R \varepsilon \left( \varpi^{(1)}_{1} \right)^{2} - (2 M(R)+R)\left( \dfrac{d\varpi^{(1)}_{1}}{dR} \right)^{2}\right) \Bigg\} \, . \end{align}\end{split}\]

Boundary conditions

The following asymptotic expansion is used for both systems to generate boundary conditions at the center,

\[\begin{split}\begin{align} h^{(2)}_{2}(R) &= B R^{2} + \mathcal{O}(R^{4}) \, , \\ k^{(2)}_{2}(R) &= -B R^{2} + \mathcal{O}(R^{4}) \, . \end{align}\end{split}\]

Tidal deformability

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\)

Numerical methods

• The module utilizes Steffen’s algorithm as a monotonic increasing interpolator to process the Equation of State (EOS) data file. This approach ensures that the EoS remains monotonic, which is crucial for maintaining the stability of the star’s equilibrium configuration. The GNU Scientific Library (GSL) is used for this interpolation. For more information on how to use the GSL Steffen’s interpolator, please refer to the GSL website.

• To solve the interior equations, QLIMR employs the well-optimized and fast RF45 ODE solver from the GNU Scientific Library (GSL). This integrator is specifically designed for efficient performance. For more details on how to use the GSL integrator, please refer to GSL integrators.

References

[1] M. D. Duez, Y. T. Liu, S. L. Shapiro, and M. Shibata (2006). Evolution of magnetized, differentially rotating neutron stars: Simulations in full general relativity, Phys. Rev. D 73, 104015

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

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

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

[5] 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.

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

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