# Relativistic correction to the r-mode frequency in light of multi-messenger
constraints

Suprovo Ghosh Inter-University Centre for Astronomy and Astrophysics,
Pune University Campus, Pune - 411007, India
Dhruv Pathak Inter-University Centre for Astronomy and Astrophysics,
Pune University Campus, Pune - 411007, India
Debarati Chatterjee Inter-University Centre for Astronomy and Astrophysics,
Pune University Campus, Pune - 411007, India

###### Abstract

R-mode oscillations of rotating neutron stars are promising candidates for
continuous gravitational wave (GW) observations. The r-mode frequencies for
slowly rotating Newtonian stars are well-known and independent of the equation
of state (EOS) but for neutron stars, several mechanisms can alter the r-mode
frequency of which the relativistic correction is dominant and relevant for
most of the neutron stars. The most sensitive searches for continuous GWs are
those for known pulsars for which GW frequencies are in targeted narrow
frequency bands of few Hz. In this study, we investigate the effect of several
state-of-the-art multi-messenger constraints on the r-mode frequency for
relativistic, slowly rotating, barotropic stars. Imposing these recent
constraints on the EOS, we find that the r-mode frequency range is slightly
higher from the previous study and the narrow band frequency range can
increase upto 8-25% for the most promising candidate PSR J0537-6910 depending
on the range of compactness. We also derive universal relations between r-mode
frequency and dimensionless tidal deformability which can be used to estimate
the dynamical tide of the r-mode resonant excitation during the inspiral
signal. These results can be used to construct the parameter space for r-mode
searches in gravitational wave data and also constrain the nuclear equation of
state following a successful r-mode detection.

neutron stars, gravitational waves, r-mode, multi-messenger observations

##  1 Introduction

Gravitational waves (GW) can drive various modes of oscillation in a rotating
Neutron Star (NS) unstable (Chandrasekhar, 1970; Friedman & Schutz, 1978). The
r-mode is a toroidal mode of fluid oscillation for which the restoring force
is the Coriolis force. For any rotating star, r-modes can become unstable to
gravitational wave emission via the Chandrasekhar-Friedman-Schutz (CFS)
mechanism (Andersson, 1998, 2003). This instability can explain the spin down
of hot and young neutron stars (Andersson et al., 1999; Lindblom et al., 1998;
Alford & Schwenzer, 2014) as well as old, accreting neutron stars in low mass
x-ray binaries (LMXBs) (Ho et al., 2011), and provides a plausible explanation
for the absence of very fast rotating neutron stars in nature. Although shear
and bulk viscosity of the NS matter can damp these oscillation modes (Lindblom
& Owen, 2002; Chatterjee & Bandyopadhyay, 2006, 2007), while spinning down,
neutron stars can enter a region of temperature and rotational frequency in
which viscosity cannot damp the r-mode and the mode grows to a large amplitude
leading to GW emission (Lindblom et al., 1998).

Because of its astrophysical significance (Abbott et al., 2022), there have
been recent searches for continuous GW emission specifically from r-modes
using the LIGO-Virgo-Kagra (LVK) global network of gravitational wave
detectors (Aasi et al., 2015; Abbott et al., 2016; Acernese et al., 2014;
Akutsu et al., 2021). A search for GWs from r-modes in the Crab pulsar was
carried out by Rajbhandari et al. (2021) and from PSR J0537 by Fesik & Papa
(2020a, b) and the LVK collaboration (Abbott et al., 2021). No GWs were
detected in these searches, but upper limits on GW amplitude were obtained.

As the r-modes are rotationally restored stellar oscillations, their frequency
is proportional to the rotation frequency (Papaloizou & Pringle, 1978) for
slowly rotating stars. If the targeted pulsar rotational frequency is known
from electromagnetic data, the GW searches are done in a relatively narrow
frequency band obtained from theoretical estimates (Caride et al., 2019). The
r-mode frequency for slowly and uniformly rotating Newtonian star in the
perfect fluid approximation is found to be independent of the Equation of
state (EOS) (Provost et al., 1981)

$\kappa\equiv\frac{\sigma_{R}}{\Omega}=\frac{2m}{l(l+1)},$ (1)

where $\sigma_{R}$ is the r-mode frequency in the co-rotating frame, $\Omega$
is the rotational angular velocity of the star, $l$ and $m$ are the spherical
harmonic indices. In the inertial frame, the frequency is given by
$f=|(\kappa-m)\Omega|$. For neutron stars, dominant factors influencing the
r-mode frequency are relativistic effects (Lockitch et al., 2000, 2003; Idrisy
et al., 2015) and rapid rotation (Lindblom et al., 1999; Yoshida et al.,
2005). There are several other factors like the presence of solid crust (Levin
& Ushomirsky, 2001), stratification (Yoshida & Lee, 2000; Passamonti et al.,
2009), magnetic field (Ho & Lai, 2000; Morsink & Rezania, 2002) or
superfluidity in the core (Lindblom & Mendell, 2000; Andersson & Comer, 2001)
that might affect the r-mode frequency but their effect was found to be
negligible for most stars (Idrisy et al., 2015). Rapid rotation can increase
the value of $\kappa$ by $\sim 6\%$ for fastest rotating stars (Idrisy et al.,
2015). But since this correction is of the order $(\Omega/\Omega_{K})^{2}$
where $\Omega_{K}$ is the Keplerian frequency, the rotational corrections
become negligible for slowly spinning stars. Lockitch et al. (2000, 2003)
first derived the perturbation equations to solve for r-mode frequency for
relativistic barotropic stars and calculated them using only polytropic EOSs.
Later Idrisy et al. (2015) calculated the same for some tabulated EOSs. They
sampled 14 EOSs that support a maximum mass of at least 1.85$M_{\odot}$ and
obtained the universal relation between r-mode frequency and compactness (Fig
3 in (Idrisy et al., 2015)) as

$\kappa=0.627+0.079(M/R)-2.25(M/R)^{2}~{}.$ (2)

Since this work, there have been several recent multi-messenger observations
of neutron stars that have put constraints on the Neutron star equation of
state. The detection of gravitational waves from the binary neutron star
merger GW170817 (Abbott et al., 2017, 2018, 2019, 2020a), the NICER
measurement of mass and radius of pulsars (Miller et al., 2019; Riley et al.,
2019; Miller et al., 2021; Riley et al., 2021) has been extensively used to
constrain the nuclear equation of state (Pang et al., 2021; Traversi et al.,
2020; Dietrich et al., 2020; Legred et al., 2021; Annala et al., 2018; Biswas
et al., 2021; Ghosh et al., 2022). The maximum observed mass of pulsar is also
updated to around $2.08M_{\odot}$ (Fonseca et al., 2021). These multi-
messenger constraints (Abbott et al., 2020a) have ruled out with good
confidence several tabulated EOSs that were used in Idrisy et al. (2015).
Also, in the universal relation Eqn. 2, they did not impose two important
physical constraints on the r-mode frequency that(i) in the Newtonian limit
when $M/R\longrightarrow 0$, then $\kappa=2/3(0.667)$ for $l=m=2$ mode and
(ii) the linear coefficient in the fit relation should be negative, otherwise
it implies that the r-mode frequency increases with increase in compactness
upto a certain value (Lockitch et al., 2000).

In this study, we calculate the r-mode frequency for 15 tabulated equations of
state sampled from (Idrisy et al., 2015; Abbott et al., 2020a) that are
consistent with recent multi-messenger observations of neutron stars. Along
with these tabulated EOSs, we also consider a posterior for the EOS from
Legred et al. (2021), employing a nonparametric EOS model based on Gaussian
processes and combining information from pulsar masses, NICER observation of
mass radius and GW observations of binary neutron star mergers. The advantages
of this nonparametric EOS model is that this allows more model freedom in the
EOS representation than any direct parametrization with a small number of
parameters; can account for different degrees of freedom, including hyperonic,
quark models, phase transition and is not subject to any systematic errors
that arise with parametrized EOS families (Legred et al., 2021). A couple of
recent studies (Gupta et al., 2022; Ma et al., 2021) have also looked into the
possible resonant r-modes detection in the inspiral phase of the binary
mergers using the third generation detector Einstein telescope (Punturo et
al., 2010; Hild et al., 2011). The r-mode introduces an additional phase to
the waveform model which is estimated from the r-mode frequency $\kappa$ as
function of the dimensionless tidal deformability of the neutron star. So, for
both of these sets of EOS models, we calculate the r-mode frequency as a
function of dimensionless tidal deformability. We also give universal
relations of the r-mode frequency with both the compactness and dimensionless
tidal deformability which will be useful to constrain the parameter space for
searches for r-modes from both isolated pulsars and excitation during inspiral
phase in binary. In case, gravitational wave is detected from the r-mode
oscillations, these universal relations can also be used to constrain the
neutron star EOS.

The structure of the article is as follows: in Sec. 2, we describe the
formalism of the structure of equilibrium configuration of slowly and
uniformly rotating star, relevant perturbation and boundary equations. The
details of the numerical scheme and their convergence are discussed in Sec. 3.
We first test the scheme by comparing with the results of Idrisy et al.
(2015); Lockitch et al. (2003) in Sec 4. We then extend the analysis for our
EOS models and demonstrate the results of the investigation, including the
universal relations in Sec. 4. Finally, in Sec. 5 we discuss the main
implications of this work.

##  2 Formalism

###  2.1 Equilibrium configuration of slowly and uniformly rotating star

We consider a slowly and uniformly rotating perfect fluid star with angular
velocity $\Omega$. The slow rotation approximation requires that
$\Omega\ll\Omega_{K}$ where $\Omega_{K}$ is the Keplerian
frequency($\propto\sqrt{M/R^{3}}$). In this slow-rotation approximation, the
neutron star retains its spherical geometry as the centrifugal deformations
are an order $\Omega^{2}$ effect (Hartle, 1967). The equilibrium solution to
the slowly and uniformly rotating star is obtained by solving the Einstein
equations $G_{\alpha\beta}=8\pi T_{\alpha\beta}$ using the line element given
by (Lockitch et al., 2000)

$\begin{split}ds^{2}&=-e^{2\nu(r)}dt^{2}+e^{2\lambda(r)}dr^{2}+r^{2}d\theta^{2}\\\
&\quad+r^{2}\sin^{2}(\theta)d\phi^{2}-2\omega(r)r^{2}\sin^{2}(\theta)dtd\phi,\end{split}$
(3)

where $\omega(r)$ is the rotational frame-dragging inside the star. The energy
momentum tensor for the perfect fluid is given by

$T_{\alpha\beta}=(\varepsilon+p)u_{\alpha}u_{\beta}+pg_{\alpha\beta}.$ (4)

Here, $\varepsilon,p$ are the total energy density and pressure of the fluid
respectively measured by an observer co-moving with the 4-velocity

$u^{\alpha}=e^{-\nu}(t^{\alpha}+\Omega\phi^{\alpha}),$ (5)

where $t^{\alpha}$ and $\phi^{\alpha}$ are the timelike and rotational killing
vectors of the spacetime. Solving the Einstein equations for the metric and
fluid variables, they reduce to Tolman-Oppenheimer-Volkov (TOV) equations
given in Eqn. (LABEL:eq:tov)

$\displaystyle\frac{dm(r)}{dr}$ $\displaystyle=$ $\displaystyle
4\pi\varepsilon(r)r^{2}~{},$ $\displaystyle\frac{dp(r)}{dr}$ $\displaystyle=$
$\displaystyle-\frac{[p(r)+\varepsilon(r)][m(r)+4\pi
r^{3}p(r)]}{r(r-2m(r))}.~{}$

We also obtain the equations for the metric functions $\nu(r)$ and
$\lambda(r)$ as

$\displaystyle\frac{d\nu(r)}{dr}$ $\displaystyle=$
$\displaystyle-\frac{1}{\varepsilon+p}\frac{dp(r)}{dr},$ $\displaystyle
e^{-2\lambda}$ $\displaystyle=$ $\displaystyle 1-\frac{2m(r)}{r}~{}.$

For a given EOS $p=p(\varepsilon)$, the TOV equations (LABEL:eq:tov) are
integrated from the centre of the star to the surface using the boundary
conditions of vanishing mass, $m|_{r=0}=0$, at the centre of the star, and a
vanishing pressure, $p|_{r=R}=0$, at the surface ($r=R$). For the metric
function $\nu(r)$, we start the integration from $r=0$ with $\nu(0)=0$ and the
solution must match the exterior solution at the surface. We implement it by
making the following variable change (Glendenning, 2012)

$\nu(r)\longrightarrow\nu(r)-\nu(R)+\frac{1}{2}ln\left(1-\frac{2M}{R}\right),$
(8)

where $M=m(R)$ is the total mass of the star.

For the slowly rotating equilibrium configuration, in addition to the TOV
equations, we need to solve an equation for the other metric function
$\omega(r)$ given by the Hartle equation (Hartle, 1967)

$\frac{1}{r^{4}}\frac{d}{dr}\left(r^{4}j\frac{d\bar{\omega}}{dr}\right)+\frac{4}{r}\frac{dj}{dr}\bar{\omega}=0,$
(9)

where

$\bar{\omega}(r)=\Omega-\omega(r).$ (10)

j(r) is defined in terms of the metric functions as

$\displaystyle j(r)$ $\displaystyle=e^{-(\nu+\lambda)}\text{\hskip 17.07164pt
for $r\leq R$,}$ $\displaystyle=1\text{\hskip 45.5244pt for $r>R$.}$ (11)

The differential equation (9) for $\omega(r)$ can be integrated from r = 0
with an arbitrary choice of the central value $\bar{\omega}(0)$ and a
vanishing slope (Glendenning, 2012). At the surface, it should match the
exterior solution. From equation (9), the exterior solution ($r>R$) is given
by

$\omega(r)=\frac{2J}{r^{3}},$ (12)

where J is the angular momentum of the star. At the surface, the corresponding
boundary conditions should be matched

$\displaystyle J$ $\displaystyle=$
$\displaystyle\frac{1}{6}R^{4}\left(\frac{d\bar{\omega}}{dr}\right)_{R},$
$\displaystyle\Omega$ $\displaystyle=$
$\displaystyle\bar{\omega}(R)+\frac{2J}{R^{3}}.$

Since, $\bar{\omega}(r)$ depends on the rotational frequency $\Omega$, we
normalise $\bar{\omega}(r)$ by $\Omega$,
$\bar{\omega}(r)\equiv\frac{\bar{\omega}(r)}{\Omega}$ and make $\Omega=1$.

The tidal deformability parameter quantifies the degree of the tidal
deformation effects due to the companion in coalescing binary NS systems
during the early stages of an inspiral. It is defined as

$\lambda=-\frac{Q_{ij}}{\varepsilon_{ij}},$ (14)

where $Q_{ij}$ is the induced mass quadrupole moment of the NS and
$\varepsilon_{ij}$ is the gravitational tidal field of the companion. The
dimensionless tidal deformability ($\Lambda$) can be obtained by solving a set
of differential equations coupled with the TOV equations and it is related to
the dimensionless $l=2$ tidal Love number $k_{2}$ (Hinderer, 2008; Flanagan &
Hinderer, 2008) as

$\Lambda=\frac{2}{3}k_{2}\left(\frac{R}{M}\right)^{5}.~{}$ (15)

###  2.2 Perturbation equations

Here, we consider the non-radial perturbations of these slowly rotating
equilibrium models to linear order in $\Omega$. Since the equilibrium
spacetime is stationary and axisymmetric, we decompose our perturbations using
the Lagrangian formalism into modes of the form $e^{i(\sigma t+m\phi)}$
(Lockitch et al., 2000). We express the perturbed configuration in terms of
the set $(h_{\alpha\beta},\zeta^{\alpha},\delta\varepsilon,\delta p)$. Since,
the perturbed energy density and pressure are scalar, they have polar parity
and given as

$\delta\varepsilon=\delta\varepsilon(r)Y_{l}^{m},\delta p=\delta
p(r)Y_{l}^{m}.$ (16)

The Lagrangian displacement vector is defined as

$\begin{split}\zeta^{\alpha}&=\frac{1}{i\kappa\Omega}\sum_{l=m}^{\infty}\left(\frac{1}{r}W_{l}(r)Y_{l}^{m}r^{\alpha}+V_{l}(r)\nabla^{\alpha}Y_{l}^{m}\right.\\\
&\left.\quad-
iU_{l}(r)P^{\alpha}_{\nu}\epsilon^{\nu\beta\gamma\delta}\nabla_{\beta}Y^{l}_{m}\nabla_{\gamma}t\nabla_{\delta}r\right)e^{i\sigma
t},\end{split}$ (17)

where

$P^{\alpha}_{\nu}=e^{\nu+\lambda}(\delta_{\nu}^{\alpha}-t_{\nu}\nabla^{\alpha}t),$
(18)

and the co-moving frequency is given by

$\kappa\Omega\equiv\sigma+m\Omega.$ (19)

The perturbation variables $W_{l}\&V_{l}$ have polar parity and $U_{l}$ has
axial parity. In the Regge-Wheeler gauge the metric perturbation is given by

$\begin{split}h_{\mu\nu}=Y_{l}^{m}e^{i\sigma t}\sum_{l=m}^{\infty}\hskip
170.71652pt&\\\
\begin{bmatrix}H_{0,l}(r)e^{2\nu}&H_{1,l}(r)&h_{0,l}(r)(\frac{m}{\sin(\theta)})&ih_{0,l}(r)\sin(\theta)\partial_{\theta}\\\
H_{1,l}(r)e^{2\nu}&H_{2,l}(r)&h_{1,l}(r)(\frac{m}{\sin(\theta)})&ih_{0,l}(r)\sin(\theta)\partial_{\theta}\\\
symm&symm&r^{2}K_{l}(r)&0\\\
symm&symm&0&r^{2}\sin^{2}(\theta)K_{l}(r)\end{bmatrix},\end{split}$ (20)

which contains both axial ($h_{0,l},h_{1,l}$) and polar
($H_{0,l},H_{1,l},H_{2,l},K_{l}$) parity components. We obtain the
perturbation equations by requiring $\delta G_{\alpha\beta}=8\pi\delta
T_{\alpha\beta}$ term by term. Since we are considering upto order $\Omega$
variables and the displacement vector in Eqn. (17) has already a
$\kappa\Omega$ term, we only keep the zeroth order perturbation variables. The
variables can be grouped depending upon their order of dependence on the
rotational frequency $\Omega$. The relevant variables which are zeroth order
in $\Omega$ are $W_{l},V_{l},U_{l},H_{1,l},h_{0,l}$. Because $h_{1,l}$ is an
order $\Omega$ variable, we drop the “0” subscript and write $h_{0,l}$ as
$h_{l}$. The relevant $O(1)$ equations are given by

$H_{1,l}+\frac{16\pi(\varepsilon+p)}{l(l+1)}e^{2\lambda}rW_{l}=0~{},$ (21)
$V_{l}[l(l+1)(\epsilon+p)]-e^{-(\nu+\lambda)}[(\epsilon+p)e^{\nu+\lambda}W_{l}]^{\prime}=0,$
(22)
$\begin{split}r^{2}h_{l}^{\prime\prime}-r^{2}(\nu^{\prime}+\lambda^{\prime})h_{l}^{\prime}&+[(2-l^{2}-l)e^{2\lambda}\\\
-r(\nu^{\prime}+\lambda^{\prime})-2]h_{l}&-4r(\nu^{\prime}+\lambda^{\prime})U_{l}=0,\end{split}$
(23)

where a prime denotes a derivative with respect to r. We have used Eq. (21) to
eliminate the variable $H_{1,l}$ in favour of $W_{l}$.

To close the system of equations we obtain two other independent equations
that arise at $O(\Omega)$ which enforces the conservation of vorticity in
constant entropy surfaces

$\displaystyle 0=[l(l+1)\kappa\Omega(h_{l}+U_{l})-2m\bar{\omega}U_{l}]$
$\displaystyle+(l+1)Q_{l}[(2\bar{\omega}+\mu)W_{l-1}-2(l-1)\bar{\omega}V_{l-1}]$
$\displaystyle-Q_{l+1}[(2\bar{\omega}+\mu)W_{l+1}+2(l+2)\bar{\omega}V_{l+1}],$
(24)
$\begin{split}0&=(l-2)l(l+1)Q_{l-1}Q_{l}\left[-2\bar{\omega}rU^{\prime}_{l-2}\right.\\\
&\left.+2(l-1)\bar{\omega}U_{l-2}+(l-3)\mu U_{l-2}\right]\\\
&+(l+1)Q_{l}\left[(l-1)l\kappa\Omega
rV^{\prime}_{l-1}-2ml\bar{\omega}rV^{\prime}_{l-1}\right.\\\
&\left.-2l(l-1)\kappa\Omega
r\nu^{\prime}V_{l-1}+2ml(l-1)\bar{\omega}V_{l-1}\right.\\\ &\left.+m(l-3)l\mu
V_{l-1}\right]-(l+1)Q_{l}\\\ &\left[(l-1)l\kappa\Omega
e^{2\lambda}W_{l-1}-4\kappa\Omega
r(\nu^{\prime}+\lambda^{\prime})W_{l-1}\right]\\\ &+l(l+1)\left[m\kappa\Omega
r(h^{\prime}_{l}+U^{\prime}_{l})-2m\kappa\Omega
r\nu^{\prime}(h_{l}+U_{l})\right.\\\
&\left.+((l+1)Q^{2}_{l}-lQ^{2}_{l+1})(2\bar{\omega}rU^{\prime}_{l}+2\mu
U_{l})\right]\\\ &+l(l+1)\left[m^{2}+l(l+1)(Q^{2}_{l+1}+Q^{2}_{l}-1)\right]\\\
&\times(2\bar{\omega}+\mu)U_{l}-lQ_{l+1}\left[(l+1)(l+2)\kappa\Omega
rV^{\prime}_{l+1}\right.\\\
&\left.+2m(l+1)\bar{\omega}rV^{\prime}_{l+1}-2(l+1)(l+2)\kappa\Omega
r\nu^{\prime}V_{l+1}\right]\\\
&-lQ_{l+1}\left[2m(l+1)(l+2)\bar{\omega}V_{l+1}\right.\\\
&\left.+m(l+1)(l+4)\mu V_{l+1}-(l+1)(l+2)\kappa\Omega
e^{2\lambda}W_{l+1}\right.\\\ &\left.+4\kappa\Omega
r(\nu^{\prime}+\lambda^{\prime})W_{l+1}\right]\\\
&+l(l+1)(l+3)Q_{l+1}Q_{l+2}\left[2\bar{\omega}rU^{\prime}_{l+2}\right.\\\
&\left.+2(l+2)\bar{\omega}U_{l+2}+(l+4)\mu U_{l+2}\right],\end{split}$

where we have defined $\mu=re^{2\nu}(\bar{\omega}e^{-2\nu})^{\prime}$. Please
note that the equation (2.2) is a simplified form of Eq. (48) from Idrisy et
al. (2015)(or Eq. (22) of Lockitch et al. (2003)) which we have taken from Eq.
(346) of Lockitch (1999). For the barotropic stars (which is the case for cold
neutron stars), the conservation of vorticity gives rise to a mixing of axial
and polar modes at zeroth order in $\Omega$ (Lockitch et al., 2000).This
suggests that the modes of barotropic stars will generically be of a hybrid
nature.

###  2.3 Boundary Conditions

In order to solve the equations (22) (23) (2.2) (2.2), we need to apply the
appropriate boundary conditions. The first boundary condition is also called
the regularity condition which says that the perturbation equations must be
regular at the centre of the star. To implement this boundary condition, we
introduce a new variable $\tilde{F_{l}}$ corresponding to each perturbation
variable $F_{l}$ as

$F_{l}(r)=\left(\frac{r}{R}\right)^{l+q}\tilde{F_{l}}(r),$ (25)

where $F_{l}$ is any one of the perturbation variables
$W_{l},V_{l},U_{l},h_{l}$. The axial parity variables ($U_{l},h_{l}$) have
$q=1$ and polar parity variables ($W_{l},V_{l}$) have $q=0$ (Lockitch et al.,
2003).

The next boundary condition comes from the fact that the Lagrangian
perturbation of the pressure is zero at the surface of the star which
translates to

$W_{l}(R)=0~{}.$ (26)

Now $h_{l}$ is the only variable defined outside the star also where it
follows

$\left(1-\frac{2M}{r}\right)\frac{d^{2}h_{l}}{dr^{2}}-\left[\frac{l(l+1)}{r^{2}}-\frac{4M}{r^{3}}\right]h_{l}=0.$
(27)

The solution of this equation is given by the hypergeometric function
$h_{l}(r)=_{2}F_{1}(l-1,l+2;2l+2;2M/r)$ (Idrisy et al., 2015). The interior
and exterior solutions are matched via the following boundary conditions

$\lim_{\epsilon\to 0}[h_{l}(R-\epsilon)-h_{l}(R+\epsilon)]=0,$ (28)
$\lim_{\epsilon\to
0}[h_{l}(R-\epsilon)h^{\prime}_{l}(R+\epsilon)-h^{\prime}_{l}(R-\epsilon)h_{l}(R+\epsilon)]=0.$
(29)

##  3 Numerical Method

To get the r mode frequencies, we need to solve the perturbation equation (22)
(23) (2.2) (2.2) along with the boundary conditions (25) (26) (28) (29). For a
given EOS, solving the TOV equations (LABEL:eq:tov) and the Hartle equation
(9) gives the necessary equilibrium variables
$p(r),\varepsilon(r),\nu(r),\lambda(r)$ and $\omega(r)$. Since the
perturbation equations are coupled in terms of $l$, we need to set a upper
limit to our $l$ value ($l_{max}$) up to which we will solve these equations.
As we are focusing on axial-led hybrid modes, we fix $l_{max}$ to an odd value
to get a closed system of equations (Idrisy et al., 2015). For this choice of
axial-led hybrid modes, the perturbation variables has contribution only from
the terms with (Lockitch et al., 2000)

$\displaystyle\text{axial parity($h_{l},U_{l}$) with }l$ $\displaystyle=$
$\displaystyle m,m+2,m+4,$ $\displaystyle\text{polar parity($W_{l},V_{l}$)
with }l$ $\displaystyle=$ $\displaystyle m+1,m+3,m+5~{}.$

So, we solve for the eigenfunctions $h_{l},U_{l},W_{l+1},V_{l+1}$ for
$l=m,m+2,...$ and set others to zero.
Instead of integrating the coupled perturbation equations, we adopt a spectral
method using the Chebyshev polynomials to find the eigenfrequency $\kappa$
similar to Lockitch et al. (2003); Idrisy et al. (2015). We express our system
of ordinary differential equations in terms of sum of basis functions,
Chebyshev polynomials in our case and using some useful identities, we reduce
these differential equations in a system of algebric equations which is then
solved using root-finding techniques to get the eigenfrequency.

###  3.1 Chebyshev Polynomials

The Chebyshev polynomials of first kind are defined in $[-1,1]$ range by

$T_{i}(y)=\cos(i\arccos(y))~{}.$ (30)

Any function S(y) expanded in terms Chebyshev polynomials are given by

$S(y)=\sum_{i=0}^{i_{max}}s_{i}T_{i}(y)-\frac{1}{2}s_{0},$ (31)

where

$\begin{split}s_{i}&=\frac{2}{i_{max}}\sum_{j=0}^{i_{max}+1}\left[S\left[\cos\left(\frac{\pi(j+\frac{1}{2})}{i_{max}}\right)\right]\right.\\\
&\left.\cos\left(\frac{\pi
i(j+\frac{1}{2})}{i_{max}}\right)\right].\end{split}$ (32)

Since our variable r is in domain [0,R], we define a new variable
$y=2\left(\frac{r}{R}\right)-1$. Now, choosing $i_{max}$ is a trade-off as
increasing $i_{max}$ gives better convergence but it also increases the number
of equations to solve. Along with this, we also make use of two other
identities for Chebyshev polynomials involving derivatives of a function and
product of functions.

If $f_{l}^{\prime}$ and $f_{l}$ are the Chebyshev coefficients of the
derivative of a function and the function itself respectively, then they are
related by

$f^{\prime}_{l,i}-f^{\prime}_{l,i+2}=2(i+1)f_{l,i+1}.$ (33)

If $b_{i}$ and $f_{l,i}$ are the Chebyshev coefficients of a background
function B(y) and perturbation variable $F_{l}(y)$ respectively, then the
Chebyshev coefficients for their products are given by

$\pi_{l,i}=\sum_{j=0}^{i_{max}}[b_{i+j}+\Theta(j-1)b_{|i-j|}]f_{l,j},$ (34)

where $\Theta$ is the step function.

After imposing the regularity condition (25),each term in the perturbation
equation can be written in terms of a background function($B(r)$) which
depends on the star equilibrium profiles and a foreground
function($\tilde{F_{l}}(r)$) which are the perturbation variables. Since,
$f^{\prime}_{l}$ and $f_{l}$ are connected by the Eq. (33), we only have
$f_{l}$ as our unknown functions. Now we expand each of them in terms of the
chebyshev polynomials and simplify them using the identities (32) (33) (34).
We also re-write the boundary conditions in the same way.

After expanding all the perturbation equations in terms of Chebyshev
polynomials, we extract the co-efficients using the identities (32) (33) (34).
This leads to a system of $2(l_{max}-3)i_{max}$ linear equations for $\kappa$
with the unknown functions $f_{l}$. We represent the system of equations as

$A(\kappa)x=0,$ (35)

where x is the vector of ${f_{l}}$s containing Chebyshev coefficients of the
perturbation variables. To incorporate the boundary conditions, we replace the
equation that came from the highest order Chebyshev coefficient $f_{i_{max}}$
by a boundary condition (Idrisy et al., 2015). For example, to implement the
boundary condition $W_{l}(R)=0$ (26), we replace the highest order Chebyshev
coefficient $f_{i_{max}}$ as

$w_{l_{i_{max}}}=\frac{1}{2}w_{l_{0}}-\sum_{i=0}^{i_{max}-1}w_{l_{i}},$ (36)

where $w_{l}$’s are the Chebyshev coefficients of the variable $W_{l}$. We
have used at at $r=R$, $y=1$ and $T_{i}(1)=1$ $\forall i$.

###  3.2 Root finding method

To find the eigenfrequency $\kappa$, we set $det(A(\kappa))=0$. This leads to
a very high degree of polynomial in $\kappa$ for any reasonable value of
$l_{max}$ and $i_{max}$ which is not solvable by any standard root-finding
techniques. We use the second root-finding algorithm given in Idrisy et al.
(2015) which uses the Singular Value Decomposition (SVD) of the matrix A,
SVD(A) = $U\Sigma V$. We vary $\kappa$ in the physically possible range
($0.67$-$0.4$) and look for the value of $\kappa$ that results in the smallest
value for the last element on the diagonal of $\Sigma$. For a particular
$l_{max}$ and $i_{max}$, we find several roots for $\kappa$. The way to
determine the correct root is to start at small values of $l_{max}$ and
$i_{max}$ and increase them step by step. We will always converge to the
correct root for any $l_{max}$ and $i_{max}$ while the others will change
unpredictably.
In Fig. 1 we plot the last diagonal element in SVD of the matrix A vs $\kappa$
for different choice of $i_{max}$ with fixed $l_{max}$ and vice versa for an
$n=1$ polytrope with compactness of 0.153.

Figure 1: Convergence of the root with changing $l_{max}$ and $i_{max}$. The
upper panel shows the convergence for varying $i_{max}$ while $l_{max}=9$ and
the lower panel shows convergence for varying $l_{max}$ for $i_{max}=10$

For the polytropes we get convergence at minimum of $l_{max}=7$ and
$i_{max}=9$ and the convergence stops above $l_{max}=13$ and $i_{max}=13$. For
both polytropes and realistic equation of states, we get convergence upto 3
orders after decimal place. Sometimes, the converged root will vary in the 3rd
decimal place with different $l_{max}$ and $i_{max}$. We take the statistical
mode of the roots for different combinations of $l_{max}$ and $i_{max}$.

##  4 Results

###  4.1 Polytropic EOS

For the Newtonian stars, the $l=m=2$ r-mode is the one expected to dominate
the gravitational wave radiation from the hot and fast rotating Neutron stars
(Andersson et al., 1999; Lindblom et al., 1998). But, Lockitch et al. (2000)
showed that for the relativistic barotropes pure $l=m=2$ does not exist. The
corresponding modes are axial-led hybrid modes with $m=2$. To test the
accuracy of our numerical code, we first consider the case for polytropic EOS
and compare our results with previous studies (Idrisy et al., 2015; Lockitch
et al., 2003). In Fig 2, we plot the r-mode frequency as a function of
compactness for uniform density star model or $n=0$ polytropic EOS and compare
the same with Fig (1) in Lockitch et al. (2003). We get a significant match
and the relative error is around $0.1-0.7\%$. In Fig 3, we plot the r-mode
frequency as a function of compactness for $n=1$ polytropic EOS along with the
quadratic fit relation between r-mode frequency and compactness given by Eq.
(69) in Idrisy et al. (2015). In this case, we see a greater deviation (around
$2-3\%$) in the r-mode frequency from the results in Idrisy et al. (2015). We
discuss the implications and possible reasons for this difference in Sec. 5

Figure 2: Normalised eigenfrequency $\kappa/\kappa_{N}$ of $m=2$ mode with
their Newtonian counterparts for uniform density star as a function of
compactness. Relative differences are shown in the bottom panel. Figure 3:
Eigenvalue $\kappa$ for equilibrium stars of $n=1$ polytrope for different
compactness. Relative differences are shown in the bottom panel.

In Fig 4, we show the best linear and quadratic model fits using Least Squares
method to our results. For quadratic models, we also use a second model where
the zeroth order term of the polynomial is fixed to $2/3$ to constrain the
fact that as $M/R\longrightarrow 0$, we reach Newtonian limit where the r-mode
frequency for $l=m=2$ mode is given by

$\kappa_{N}=\frac{2}{m+1}=2/3~{}.$

Figure 4: Linear and Quadratic fits to $\kappa$ vs compactness results for
equilibrium stars of $n=1$ polytrope

From the $R^{2}$ value of the fits, we see that any quadratic model is a
better fit than the linear model. So, we use the quadratic model which
satisfies the Newtonian limit

$\kappa=0.667-0.461C-1.129C^{2}.$ (37)

We also see that along with the negative quadratic term, we also have the
linear term negative in our fit equation (37). This also satisfies the
physical constraint that the r-mode frequency should decrease as we increase
compactness (Lockitch et al., 2000) unlike the Eq. (69) in Idrisy et al.
(2015) which implies that $\kappa$ increases with increasing compactness for
$M/R<0.10$ due to the positive linear term.

###  4.2 Realistic equations of state

####  4.2.1 Tabulated equation of state

For realistic equation of states, (Idrisy et al., 2015) considered 14 EOSs
under the constraint that the EOS could support a minimum of $1.85M_{\odot}$.
From the recent observations, the heaviest known pulsar PSR J0740+6620 has a
maximum mass of $2.08^{+0.07}_{-0.07}$ $M_{\odot}$ (Fonseca et al., 2021).
Taking a upper limit of 1-$\sigma$ confidence interval of this maximum
observed pulsar mass, we only consider EOS that support a minimum of
$2.01M_{\odot}$ neutron star which rules out 3 EOSs - GNH3,BBB2 and ALF4 used
in Idrisy et al. (2015). Also, the recent analyses of the GW170817 event
(Abbott et al., 2019) apply a constraint on the upper bound of the effective
tidal deformability $\tilde{\Lambda}<$ 720 (Tong et al., 2020) using the
PhenomPNRTwaveform model and low-spin highest posterior density interval for
tidal deformability. Using these multi-messenger observations of neutron
stars, 2 very stiff EOSs (MS1,MS1b) considered in Idrisy et al. (2015) have
been ruled out with good confidence (Abbott et al., 2020a; Biswas, 2022).
Along with the 9 remaining EOSs from Idrisy et al. (2015), we consider 6
additional EOSs that satisfy the multi-messenger constraints. Out of these
additional 6 EOSs, one model QHC19 (Baym et al., 2019) incorporates a
transition between a hadronic phase in the crust and a quark matter phase in
the core and one other EOS model CMF5 (Dexheimer & Schramm, 2008) includes
nucleons and hyperons. The other equation of state models describe purely
nucleonic matter. All our EOS tables are obtained from either COMPOSE (Oertel
et al., 2017a, b) or an EOS catalogue from Özel & Freire (2016) used in
LALSuite (LIGO Scientific Collaboration, 2018). For the EOS tables obtained
from COMPOSE, we use their first-order interpolation option (Oertel et al.,
2017b) and for the tables from Özel & Freire (2016) we use standard cubic
spline interpolation (Abbott et al., 2020a). All the EOSs with their maximum
mass and radius at $1.4M_{\odot}$ are listed in Table 1.

We compute the r-mode frequency $\kappa$ for these EOS models for a range of
compactness varying from 0.1 to a maximum of 0.31. This compactness spans the
range of possible neutron stars. In Fig 5, we plot the r-mode frequency
$\kappa$ vs compactness for all the 15 EOSs and see that $\kappa$ does not
change much as a function of compactness for different EOSs.

Figure 5:  $\kappa$ vs compactness results for all 15 EOSs. The black line
represents the best fir to our quadratic model.

We use a quadratic fit with the zeroth term fixed to its Newtonian value $2/3$
and get the universal relation for $\kappa$ vs compactness C as

$\kappa=0.667-0.478C-1.11C^{2}.$ (38)

To compare the results from Gupta et al. (2022); Ma et al. (2021) where they
looked for possible resonant r-modes detection in the inspiral phase of the
binary mergers using the third generation detector Einstein telescope, in Fig
6, we plot $\kappa$ as a function of $log(\Lambda)$ for all the 15 EOSs and
fit it to a quadratic model.

Figure 6:  $\kappa$ vs $log(\Lambda)$ results for all 15 EOSs. The black line
represents the best fir to our quadratic model.

The corresponding universal relation is given by

$\kappa=0.3446+0.0471\>log(\Lambda)-0.002\>log^{2}(\Lambda).$ (39)

Comparing the result from Gupta et al. (2022), we get around $5$-$6\%$
difference in the r-mode frequency for a given tidal deformability. In Table
1, we also report the value of r-mode frequency $\kappa$ for a compactness of
0.15 and a tidal deformability of 400 for each of the 15 EOSs.

Table 1: List of all the tabulated EOSs with corresponding maximum mass, radius for $1.4M_{\odot}$(in km) star, $\kappa$ for a compactness of 0.15 (along with the same from Idrisy et al. (2015)) and tidal deformability of 400. EOS | $M_{max}$ | $R_{1.4M_{\odot}}$ |  $\kappa_{.15C}$( Idrisy et al. (2015)) | $\kappa_{400\Lambda}$
---|---|---|---|---
SLY | 2.05 | 11.77 | 0.573(0.587) | 0.558
AP3 | 2.39 | 12.06 | 0.572(0.588) | 0.555
AP4 | 2.21 | 11.40 | 0.571(0.587) | 0.557
ALF2 | 2.09 | 13.18 | 0.571(0.588) | 0.552
WFF1 | 2.13 | 10.40 | 0.572(0.587) | 0.551
WFF2 | 2.20 | 11.14 | 0.571(0.587) | 0.554
MPA1 | 2.46 | 12.48 | 0.571(0.588) | 0.553
ENG | 2.25 | 12.08 | 0.572(0.588) | 0.546
H4 | 2.03 | 12.95 | 0.573(0.591) | 0.552
DD2 | 2.42 | 13.04 | 0.572(-) | 0.554
BSK22 | 2.26 | 12.73 | 0.571(-) | 0.554
SLY9 | 2.16 | 12.20 | 0.574(-) | 0.555
QHC19 | 2.18 | 11.35 | 0.574(-) | 0.553
CMF5 | 2.07 | 12.87 | 0.573(-) | 0.554
SKI4 | 2.17 | 12.10 | 0.574(-) | 0.556

####  4.2.2 Nonparametric EOS model

Several different parametrizations of the neutron star EOS have been proposed
to constrain the EOS from multi-messenger observations of the neutron stars.
Generic parametrizations in terms of piecewise polytropes (Annala et al.,
2018; Hebeler et al., 2013; Read et al., 2009; Gamba et al., 2019), spectral
decomposition (Fasano et al., 2019; Lindblom, 2018) and speed-of-sound (Tews
et al., 2018; Greif et al., 2019) have been extensively used for such studies.
Here we consider the nonparametric representation constructed through Gaussian
process (Landry & Essick, 2019; Essick et al., 2020) rather than a
parametrization for the EOS. This EOS representation allows more model freedom
and can account for different degrees of freedom. In Legred et al. (2021),
using this model they studied the implications of the following multi-
messenger observations for the neutron star equation of state (EOS)-

  * •

the radio mass measurements for J0348+0432 (Antoniadis et al., 2013) and
J0740+6620 (Cromartie et al., 2019; Fonseca et al., 2021).

  * •

the GW mass and tidal deformability measurementsfrom GW170817 (Abbott et al.,
2017, 2018, 2019) and GW190425 (Abbott et al., 2020b)

  * •

mass and radius constraints from NICER observations of J0030+0451 (Miller et
al., 2019; Riley et al., 2019) and J0740+6620 (Miller et al., 2021; Riley et
al., 2021)

Using these combined radio, GW and x-ray data, they put constraints for the
microscopic EOS and the macroscopic NS properties, masses, radius and tidal
deformabilities. In Fig 7, we plot the $95\%$ confidence interval for the
pressure-energy density relation using the posterior obtained by combining all
the constraints (Legred et al., 2022).

Figure 7: 95$\%$ confidence interval for pressure-density relation with
constraints from GW, NICER and PSR mass measurement

From this posterior set of 10,000 EOSs, we choose 1000 EOS randomly and
calculate the r-mode frequency as a function of compactness and tidal
deformability. In Fig. 8, we plot the $95\%$ confidence interval for the
r-mode frequency $\kappa$ as a function of compactness and tidal
deformability. We find a substantial spread unlike the Fig. 5 for tabulated
EOSs owing to the fact that we now have considered 1000 EOSs which span a wide
region in the pressure-density relation 7.

Figure 8: The upper panel shows the $95\%$ level of the r-mode frequency for a
range of compactness and the best quadratic model fit to the posterior. The
lower panel shows the same as a function of tidal deformability

The best fits for the r-mode frquency $\kappa$ to our quadratic models are
given by

$\kappa=0.67-0.593C-0.772C^{2},$ (40)
$\kappa=0.3517+0.0463\>log(\Lambda)-0.0019\>log^{2}(\Lambda).$ (41)

##  5 Discussion and Conclusion

In this work, we solved perturbation equations for relativistic barotropic
stars consistently with the boundary conditions using a spectral method and
calculated the r-mode frequency as a function of compactness for various
equation of states. First we checked the convergence of our code and then
compared our results for $n=0$ polytrope with Fig. 1 from Lockitch et al.
(2003) and for $n=1$ polytropes with Fig. 2 from Idrisy et al. (2015). We see
that our results match up to $0.1\%$ for $n=0$ polytrope from Lockitch et al.
(2003). Since the precision of the converged $\kappa$ value in our case is
upto 3 decimal places, an error of the order $\leq 1\%$ is expected. But we
see a greater deviation $\sim 2-3\%$ for the results for $n=1$ polytrope from
Idrisy et al. (2015). Probable reason for this much larger deviation could be
an erroneous representation of the perturbation and boundary condition (Eq. 41
and Eq. 50) in Idrisy et al. (2015) in comparison with Eq. 20 and Eq. A4 in
Lockitch et al. (2003) respectively. From Table 1 we find that, while
considering the realistic EOS tables, there is also around $3{\%}$ deviation
in the r-mode frequency from Idrisy et al. (2015).

To get the frequency band to look for r-modes from astrophysical sources, we
should consider the spread of compactness in the neutron star population. The
compactness is estimated by the ratio of its stellar mass to the radius
(Idrisy et al., 2015)

$\frac{M}{R}\approx
0.207\left(\frac{M}{1.4M_{\odot}}\right)\left(\frac{10\text{ km}}{R}\right).$
(42)

The minimum value of neutron star mass is taken to be 1$M_{\odot}$ from Fig 1
of Lattimer & Prakash (2010). The limit is conservative in a sense that
stellar cores with lesser mass probably would not go through supernova
explosion to produce neutron stars. For our set of tabulated EOSs listed in
Table 1, we find the maximum radius for 1$M_{\odot}$ stars to be around
$14.5$km. This gives the lower limit of compactness to be $\approx 0.103$.
Assuming the casuality of EOSs, there is a upper limit to the compactness
$\frac{M}{R}\leq 0.35$ (Lattimer & Prakash, 2007) but none of the EOSs from
Table 1 reach this high value for any stable configuration. Idrisy et al.
(2015) uses the compactness range $0.11-0.31$, but we give our limits of the
r-mode frequency for both compactness ranges in Table 2.

Now, to calculate the narrow frequency band to search over in LIGO data, the
r-mode frequency($f$) is given in terms of the rotational frequency($\nu$)
(Caride et al., 2019; Abbott et al., 2021) as

$\frac{f}{\nu}=A-B\left(\frac{\nu}{\nu_{K}}\right)^{2},$ (43)

where $\nu_{K}$ is the Keplerian frequency = 506 Hz considered in Abbott et
al. (2021). The uncertainties in the range of $A$ and $B$ give the parameter
space for the GW signal model (Caride et al., 2019; Abbott et al., 2021).
General relativistic corrections for slowly rotating stars give the range of A
(Idrisy et al., 2015) and rapid rotation correction gives the range of B
(Yoshida et al., 2005). The ranges are given by (Caride et al., 2019)

$\displaystyle 1.39$ $\displaystyle\leq A\leq 1.57,$ $\displaystyle 0$
$\displaystyle\leq B\leq 0.195~{}.$ (44)

Since, the r-mode frequency in the inertial frame is given by
$f=|(\kappa-m)\Omega|$, we can use the limits on $\kappa$ to update the limits
on $A$ value which is also listed in Table 2. Considering all these multi-
messenger constraints and a conservative limit on the possible ranges of
compactness, we put a limit on the value of A to be $1.39\leq A\leq 1.64$
which should be used for narrow-band gravitational wave searches for known
pulsars.

Table 2: Proposed ranges of $\kappa$ and $A$ (Caride et al., 2019) for ranges of compactness comparing universal relations from Idrisy et al. (2015) with our results for both Tabulated EOS and Non-parametric EOSs Universal relation | $\kappa$ | $A$
---|---|---
(compactness range) |  |
(Idrisy et al., 2015)(0.11-0.31) | 0.614-0.433 | 1.39-1.57
Tabulated EOS(38)(0.11-0.31) | 0.601-0.412 | 1.40-1.59
Tabulated EOS(38)(0.10-0.35) | 0.608-0.364 | 1.39-1.64
Non-parametric EOS(40)(0.11-0.31) | 0.596-0.415 | 1.39-1.59
Non-parametric EOS(40)(0.10-0.35) | 0.604-0.371 | 1.40-1.63

PSR J0537-6910 is particularly interesting for r-mode searches because it is
the fastest-spinning known young pulsar with rotation frequency $\nu$ =62 Hz
which places gravitational-wave frequency in the LIGO sensitivity band and its
inter-glitch braking index $n\approx 7$ which is expected for GW emission via
r-mode (Andersson et al., 2018). In the latest LIGO r-mode search from PSR
J0537, values of $A$ and $B$ given in (5) gives the r-mode frequency band
86-98 Hz. If we use our universal relation (38) for a compactness range of
$0.11-0.31$ (same as in Idrisy et al. (2015)), we get the frequency band 87-99
Hz which is 1 Hz higher in both lower and upper bounds. This might not appear
very large but since we are looking at narrow frequency band searches for
these analyses, this new universal relation introduces a $8.3\%$ change (1 Hz
in a frequency band of 12 Hz) in the frequency band. For a compactness range
of $0.10-0.35$, this range becomes 86-101 Hz.

We also provided universal relations for the r-mode frequency $\kappa$ and the
dimsionless tidal deformability($\Lambda$) in Eq. (39) and Eq. (41). Recent
work by Gupta et al. (2022) used such a relation to reconstruct the EOS with
observations of the inspiral signal by Einstein Telescope (ET) with or without
r-modes. They used the universal relation for $\kappa$ vs compactness from
Idrisy et al. (2015) to derive the r-mode frequency as a function of
$log(\Lambda)$. If we compare with Eq. (7) from Gupta et al. (2022) against
our universal relation, we observe a $5-6\%$ difference in the r-mode
frequency which can be significant while constraining the nuclear EOS from
inspiral signal with r modes.

To conclude, we have derived the r-mode frequency as a function of compactness
for neutron stars and showed that, one can obtain a universal relation between
r-mode frequency and the compactness which is independent of the EOSs. With a
physically motivated range of compactness, we derived the frequency band to
search for r-modes in the LIGO data. For the particular interesting candidate
PSR J0537-6910, we showed that our narrow-band frequency range can vary from
the previous searches by $8-25\%$ depending on the compactness range chosen.
If a continuous wave from r-mode is detected, these universal relations can
also be used to constrain the nuclear EOS but distinguishing between different
EOSs will be difficult as the deviation in r-mode frequency is $<1\%$ for
different EOSs. Using ET, if r-mode excitation is observed in the inspiral
signal, we also can constrain the tidal deformability and the EOS using our
universal relations. Although in this study, we ignored rapid rotation and
other physical mechanisms inside the neutron star that can affect the r-mode
frequency, they might be important for some particular neutron stars and in
future, we would like to explore how these mechanisms affect r-mode frequency
for different EOSs.

##  6 Acknowledgements

The authors would like to thank David Ian Jones for his useful comments and
suggestions on the material of this paper. The authors acknowledge usage of
the IUCAA HPC computing facility for the numerical calculations.

## References

  * Aasi et al. (2015) Aasi, J., Abbott, B. P., Abbott, R., et al. 2015, Classical and Quantum Gravity, 32, 074001, doi: 10.1088/0264-9381/32/7/074001
  * Abbott et al. (2017) Abbott, B., Abbott, R., Abbott, T., et al. 2017, Physical Review Letters, 119, doi: 10.1103/physrevlett.119.161101
  * Abbott et al. (2018) —. 2018, Physical Review Letters, 121, doi: 10.1103/physrevlett.121.161101
  * Abbott et al. (2019) —. 2019, Phys. Rev. X, 9, 011001, doi: 10.1103/PhysRevX.9.011001
  * Abbott et al. (2020a) —. 2020a, Classical and Quantum Gravity, 37, 045006, doi: 10.1088/1361-6382/ab5f7c
  * Abbott et al. (2020b) —. 2020b, The Astrophysical Journal Letters, 892, L3, doi: 10.3847/2041-8213/ab75f5
  * Abbott et al. (2022) —. 2022, Phys. Rev. D, 105, 022002, doi: 10.1103/PhysRevD.105.022002
  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 131103, doi: 10.1103/PhysRevLett.116.131103
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, The Astrophysical Journal, 922, 71, doi: 10.3847/1538-4357/ac0d52
  * Acernese et al. (2014) Acernese, F., Agathos, M., Agatsuma, K., et al. 2014, Classical and Quantum Gravity, 32, 024001, doi: 10.1088/0264-9381/32/2/024001
  * Akutsu et al. (2021) Akutsu, T., Ando, M., Arai, K., et al. 2021, Progress of Theoretical and Experimental Physics, 2021, 05A101, doi: 10.1093/ptep/ptaa125
  * Alford & Schwenzer (2014) Alford, M. G., & Schwenzer, K. 2014, The Astrophysical Journal, 781, 26, doi: 10.1088/0004-637x/781/1/26
  * Andersson (1998) Andersson, N. 1998, The Astrophysical Journal, 502, 708, doi: 10.1086/305919
  * Andersson (2003) Andersson, N. 2003, Classical and Quantum Gravity, 20, R105, doi: 10.1088/0264-9381/20/7/201
  * Andersson et al. (2018) Andersson, N., Antonopoulou, D., Espinoza, C. M., Haskell, B., & Ho, W. C. G. 2018, The Astrophysical Journal, 864, 137, doi: 10.3847/1538-4357/aad6eb
  * Andersson & Comer (2001) Andersson, N., & Comer, G. 2001, Monthly Notices of the Royal Astronomical Society, 328, 1129, doi: 10.1046/j.1365-8711.2001.04923.x
  * Andersson et al. (1999) Andersson, N., Kokkotas, K., & Schutz, B. F. 1999, The Astrophysical Journal, 510, 846, doi: 10.1086/306625
  * Annala et al. (2018) Annala, E., Gorda, T., Kurkela, A., & Vuorinen, A. 2018, Phys. Rev. Lett., 120, 172703, doi: 10.1103/PhysRevLett.120.172703
  * Antoniadis et al. (2013) Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 1233232, doi: 10.1126/science.1233232
  * Baym et al. (2019) Baym, G., Furusawa, S., Hatsuda, T., Kojo, T., & Togashi, H. 2019, The Astrophysical Journal, 885, 42, doi: 10.3847/1538-4357/ab441e
  * Biswas (2022) Biswas, B. 2022, The Astrophysical Journal, 926, 75, doi: 10.3847/1538-4357/ac447b
  * Biswas et al. (2021) Biswas, B., Char, P., Nandi, R., & Bose, S. 2021, Physical Review D, 103, doi: 10.1103/physrevd.103.103015
  * Caride et al. (2019) Caride, S., Inta, R., Owen, B. J., & Rajbhandari, B. 2019, Phys. Rev. D, 100, 064013, doi: 10.1103/PhysRevD.100.064013
  * Chandrasekhar (1970) Chandrasekhar, S. 1970, Phys. Rev. Lett., 24, 611, doi: 10.1103/PhysRevLett.24.611
  * Chatterjee & Bandyopadhyay (2006) Chatterjee, D., & Bandyopadhyay, D. 2006, Phys. Rev. D, 74, 023003, doi: 10.1103/PhysRevD.74.023003
  * Chatterjee & Bandyopadhyay (2007) —. 2007, Phys. Rev. D, 75, 123006, doi: 10.1103/PhysRevD.75.123006
  * Cromartie et al. (2019) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2019, Nature Astronomy, 4, 72–76, doi: 10.1038/s41550-019-0880-2
  * Dexheimer & Schramm (2008) Dexheimer, V., & Schramm, S. 2008, The Astrophysical Journal, 683, 943, doi: 10.1086/589735
  * Dietrich et al. (2020) Dietrich, T., Coughlin, M. W., Pang, P. T. H., et al. 2020, Science, 370, 1450–1453, doi: 10.1126/science.abb4317
  * Essick et al. (2020) Essick, R., Landry, P., & Holz, D. E. 2020, Phys. Rev. D, 101, 063007, doi: 10.1103/PhysRevD.101.063007
  * Fasano et al. (2019) Fasano, M., Abdelsalhin, T., Maselli, A., & Ferrari, V. 2019, Phys. Rev. Lett., 123, 141101, doi: 10.1103/PhysRevLett.123.141101
  * Fesik & Papa (2020a) Fesik, L., & Papa, M. A. 2020a, The Astrophysical Journal, 895, 11, doi: 10.3847/1538-4357/ab8193
  * Fesik & Papa (2020b) —. 2020b, The Astrophysical Journal, 897, 185, doi: 10.3847/1538-4357/aba04e
  * Flanagan & Hinderer (2008) Flanagan, E. E., & Hinderer, T. 2008, Phys. Rev. D, 77, 021502, doi: 10.1103/PhysRevD.77.021502
  * Fonseca et al. (2021) Fonseca, E., et al. 2021, Astrophys. J. Lett., 915, L12, doi: 10.3847/2041-8213/ac03b8
  * Friedman & Schutz (1978) Friedman, J. L., & Schutz, B. F. 1978, The Astrophysical Journal, 222, 281, doi: 10.1086/156143
  * Gamba et al. (2019) Gamba, R., Read, J. S., & Wade, L. E. 2019, Classical and Quantum Gravity, 37, 025008, doi: 10.1088/1361-6382/ab5ba4
  * Ghosh et al. (2022) Ghosh, S., Chatterjee, D., & Schaffner-Bielich, J. 2022, The European Physical Journal A, 58, doi: 10.1140/epja/s10050-022-00679-w
  * Glendenning (2012) Glendenning, N. 2012, Compact Stars: Nuclear Physics, Particle Physics and General Relativity, Astronomy and Astrophysics Library (Springer New York).  https://books.google.co.in/books?id=cCDlBwAAQBAJ
  * Greif et al. (2019) Greif, S. K., Raaijmakers, G., Hebeler, K., Schwenk, A., & Watts, A. L. 2019, Monthly Notices of the Royal Astronomical Society, 485, 5363–5376, doi: 10.1093/mnras/stz654
  * Gupta et al. (2022) Gupta, P. K., Puecher, A., Pang, P. T. H., et al. 2022, Determining the equation of state of neutron stars with Einstein Telescope using tidal effects and r-mode excitations from a population of binary inspirals, arXiv, doi: 10.48550/ARXIV.2205.01182
  * Hartle (1967) Hartle, J. B. 1967, The Astrophysical Journal, 150, 1005, doi: 10.1086/149400
  * Hebeler et al. (2013) Hebeler, K., Lattimer, J. M., Pethick, C. J., & Schwenk, A. 2013, The Astrophysical Journal, 773, 11, doi: 10.1088/0004-637x/773/1/11
  * Hild et al. (2011) Hild, S., Abernathy, M., Acernese, F., et al. 2011, Classical and Quantum Gravity, 28, 094013, doi: 10.1088/0264-9381/28/9/094013
  * Hinderer (2008) Hinderer, T. 2008, The Astrophysical Journal, 677, 1216–1220, doi: 10.1086/533487
  * Ho et al. (2011) Ho, W. C. G., Andersson, N., & Haskell, B. 2011, Phys. Rev. Lett., 107, 101101, doi: 10.1103/PhysRevLett.107.101101
  * Ho & Lai (2000) Ho, W. C. G., & Lai, D. 2000, The Astrophysical Journal, 543, 386, doi: 10.1086/317085
  * Idrisy et al. (2015) Idrisy, A., Owen, B. J., & Jones, D. I. 2015, Physical Review D, 91, doi: 10.1103/physrevd.91.024001
  * Landry & Essick (2019) Landry, P., & Essick, R. 2019, Phys. Rev. D, 99, 084049, doi: 10.1103/PhysRevD.99.084049
  * Lattimer & Prakash (2007) Lattimer, J. M., & Prakash, M. 2007, Physics Reports, 442, 109, doi: https://doi.org/10.1016/j.physrep.2007.02.003
  * Lattimer & Prakash (2010) —. 2010, What a Two Solar Mass Neutron Star Really Means, arXiv, doi: 10.48550/ARXIV.1012.3208
  * Legred et al. (2021) Legred, I., Chatziioannou, K., Essick, R., Han, S., & Landry, P. 2021, Physical Review D, 104, doi: 10.1103/physrevd.104.063003
  * Legred et al. (2022) —. 2022, Impact of the PSR J0740+6620 radius constraint on the properties of high-density matter: Neutron star equation of state posterior samples, Zenodo, doi: 10.5281/zenodo.6502467
  * Levin & Ushomirsky (2001) Levin, Y., & Ushomirsky, G. 2001, Monthly Notices of the Royal Astronomical Society, 324, 917, doi: 10.1046/j.1365-8711.2001.04323.x
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration. 2018, LIGO Algorithm Library - LALSuite, free software (GPL), doi: 10.7935/GT1W-FZ16
  * Lindblom (2018) Lindblom, L. 2018, Phys. Rev. D, 97, 123019, doi: 10.1103/PhysRevD.97.123019
  * Lindblom & Mendell (2000) Lindblom, L., & Mendell, G. 2000, Phys. Rev. D, 61, 104003, doi: 10.1103/PhysRevD.61.104003
  * Lindblom et al. (1999) Lindblom, L., Mendell, G., & Owen, B. J. 1999, Phys. Rev. D, 60, 064006, doi: 10.1103/PhysRevD.60.064006
  * Lindblom & Owen (2002) Lindblom, L., & Owen, B. J. 2002, Phys. Rev. D, 65, 063006, doi: 10.1103/PhysRevD.65.063006
  * Lindblom et al. (1998) Lindblom, L., Owen, B. J., & Morsink, S. M. 1998, Phys. Rev. Lett., 80, 4843, doi: 10.1103/PhysRevLett.80.4843
  * Lockitch (1999) Lockitch, K. H. 1999, Stability and rotational mixing of modes in Newtonian and relativistic stars.  https://arxiv.org/abs/gr-qc/9909029
  * Lockitch et al. (2000) Lockitch, K. H., Andersson, N., & Friedman, J. L. 2000, Phys. Rev. D, 63, 024019, doi: 10.1103/PhysRevD.63.024019
  * Lockitch et al. (2003) Lockitch, K. H., Friedman, J. L., & Andersson, N. 2003, Phys. Rev. D, 68, 124010, doi: 10.1103/PhysRevD.68.124010
  * Ma et al. (2021) Ma, S., Yu, H., & Chen, Y. 2021, Physical Review D, 103, doi: 10.1103/physrevd.103.063020
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, The Astrophysical Journal, 887, L24, doi: 10.3847/2041-8213/ab50c5
  * Miller et al. (2021) —. 2021, The Astrophysical Journal Letters, 918, L28, doi: 10.3847/2041-8213/ac089b
  * Morsink & Rezania (2002) Morsink, S. M., & Rezania, V. 2002, The Astrophysical Journal, 574, 908, doi: 10.1086/341190
  * Oertel et al. (2017a) Oertel, M., Hempel, M., Klaehn, T., & Typel, S. 2017a, (https://compose.obspm.fr/)
  * Oertel et al. (2017b) —. 2017b, Rev. Mod. Phys., 89, 015007, doi: 10.1103/RevModPhys.89.015007
  * Özel & Freire (2016) Özel, F., & Freire, P. 2016, Masses, Radii, and the Equation of State of Neutron Stars, (xtreme.as.arizona.edu/NeutronStars)
  * Pang et al. (2021) Pang, P. T. H., Tews, I., Coughlin, M. W., et al. 2021, The Astrophysical Journal, 922, 14, doi: 10.3847/1538-4357/ac19ab
  * Papaloizou & Pringle (1978) Papaloizou, J., & Pringle, J. E. 1978, Monthly Notices of the Royal Astronomical Society, 182, 423, doi: 10.1093/mnras/182.3.423
  * Passamonti et al. (2009) Passamonti, A., Haskell, B., Andersson, N., Jones, D. I., & Hawke, I. 2009, Monthly Notices of the Royal Astronomical Society, 394, 730, doi: 10.1111/j.1365-2966.2009.14408.x
  * Provost et al. (1981) Provost, J., Berthomieu, G., & Rocca, A. 1981, Astronomy and Astrophysics, 94, 126
  * Punturo et al. (2010) Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Classical and Quantum Gravity, 27, 084007, doi: 10.1088/0264-9381/27/8/084007
  * Rajbhandari et al. (2021) Rajbhandari, B., Owen, B. J., Caride, S., & Inta, R. 2021, Physical Review D, 104, doi: 10.1103/physrevd.104.122008
  * Read et al. (2009) Read, J. S., Lackey, B. D., Owen, B. J., & Friedman, J. L. 2009, Physical Review D, 79, doi: 10.1103/physrevd.79.124032
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, The Astrophysical Journal, 887, L21, doi: 10.3847/2041-8213/ab481c
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, The Astrophysical Journal Letters, 918, L27, doi: 10.3847/2041-8213/ac0a81
  * Tews et al. (2018) Tews, I., Margueron, J., & Reddy, S. 2018, Physical Review C, 98, doi: 10.1103/physrevc.98.045804
  * Tong et al. (2020) Tong, H., Zhao, P., & Meng, J. 2020, Phys. Rev. C, 101, 035802, doi: 10.1103/PhysRevC.101.035802
  * Traversi et al. (2020) Traversi, S., Char, P., & Pagliara, G. 2020, The Astrophysical Journal, 897, 165, doi: 10.3847/1538-4357/ab99c1
  * Yoshida & Lee (2000) Yoshida, S., & Lee, U. 2000, The Astrophysical Journal Supplement Series, 129, 353, doi: 10.1086/313410
  * Yoshida et al. (2005) Yoshida, S., Yoshida, S., & Eriguchi, Y. 2005, Monthly Notices of the Royal Astronomical Society, 356, 217, doi: 10.1111/j.1365-2966.2004.08436.x