###### Abstract

The I-Love-Q relations are approximate equation-of-state independent relations
that connect the moment of inertia, the spin-induced quadrupole moment, and
the tidal deformability of neutron stars. In this paper, we study the I-Love-Q
relations for superfluid neutron stars for a general relativistic two-fluid
model: one fluid being the neutron superfluid and the other a conglomerate of
all charged components. We study to what extent the two-fluid dynamics might
affect the robustness of the I-Love-Q relations by using a simple two-
component polytropic model and a relativistic mean field model with
entrainment for the equation-of-state. Our results depend crucially on the
spin ratio $\Omega_{\rm n}/\Omega_{\rm p}$ between the angular velocities of
the neutron superfluid and the normal component. We find that the I-Love-Q
relations can still be satisfied to high accuracy for superfluid neutron stars
as long as the two fluids are nearly co-rotating $\Omega_{\rm n}/\Omega_{\rm
p}\approx 1$. However, the deviations from the I-Love-Q relations increase as
the spin ratio deviates from unity. In particular, the deviation of the Q-Love
relation can be as large as $O(10\%)$ if $\Omega_{\rm n}/\Omega_{\rm p}$
differ from unity by a few tens of percent. As $\Omega_{\rm n}/\Omega_{\rm
p}\approx 1$ is expected for realistic neutron stars, our results suggest that
the two-fluid dynamics should not affect the accuracy of any gravitational
waveform models for neutron star binaries that employ the relation to connect
the spin-induced quadrupole moment and the tidal deformability.

###### keywords:

neutron stars; superfluidity; general relativity; gravitational waves

1 1 0 https://doi.org/ The I-Love-Q Relations for Superfluid Neutron Stars The
I-Love-Q Relations for Superfluid Neutron Stars Cheung-Hei Yeung 1, Lap-Ming
Lin 1,*, Nils Andersson 2 and Greg Comer 3 Cheung-Hei Yeung, Lap-Ming Lin,
Nils Andersson and Greg Comer Yeung, C.-H.; Lin, L.-M.; Andersson, N.; Comer,
G.L<EMAIL_ADDRESS>
##  1 Introduction

The groundbreaking detection of the first gravitational wave signal from a
binary neutron star system GW170817 LIGO Scientific Collaboration et al.
(2017) has opened up a powerful channel to study the internal structures of
neutron stars and the poorly understood supranuclear equation-of-state (EOS).
In particular, the signal from the GW170817 event has been used to constrain
the tidal deformability and the implications for EOS models have also been
studied extensively (see, e.g., Annala et al. (2018); De et al. (2018);
Fattoyev et al. (2018); Most et al. (2018); Tews et al. (2018); Lim and Holt
(2018); Malik et al. (2018); Li et al. (2019); Carson et al. (2019)).
Observations of further gravitational wave events involving binary systems, or
even isolated neutron stars, with increasing precision will help us get closer
to unlocking the secrets of high density nuclear matter.

Observations of neutron stars can provide important information about the
poorly understood nuclear matter EOS due to the fact that the properties of
neutron stars, in general, depend sensitively on the matter model. It is thus
quite surprising that various approximately EOS-insensitive relations
connecting different neutron star properties have been discovered in the past
decade Tsui and Leung (2005); Lau et al. (2010); Yagi and Yunes (2013a, b);
Yagi (2014); Chan et al. (2014); Chakrabarti et al. (2014); Pappas and
Apostolatos (2014); Pappas (2015); Breu and Rezzolla (2016); Bozzola et al.
(2017); Luk and Lin (2018); Riahi et al. (2019); Sun et al. (2020). These
relations are “universal” in the sense that they are insensitive to the EOS
models to the $O(1\%)$ level (see Yagi and Yunes (2017); Doneva and Pappas
(2018) for reviews). These universal relations can be powerful tools to infer
the physical quantities of neutron stars. For instance, by making use of the
universal relation for the $f$-modes of neutron stars Lau et al. (2010), the
mass and radius of an isolated neutron star can be inferred accurately if
gravitational waves emitted from the $f$-mode oscillations of the star can be
detected.

Universal relations can also help reduce the number of parameters in
theoretical gravitational waveform models for binary neutron star inspirals
Lackey et al. (2019); Schmidt and Hinderer (2019); Andersson and Pnigouras
(2019); Barkett et al. (2020). It is noted that the I-Love-Q relations
discovered by Yagi and Yunes Yagi and Yunes (2013a, b) connect the moment of
inertia $I$, the $l=2$ tidal deformability, and the spin-induced quadrupole
moment $Q$ of slowly rotating neutron stars. The multipole Love relation
discovered by Yagi Yagi (2014) connects the $l=2$ and $l=3$ tidal
deformabilities. The f-mode-Love relations found by Chan et al. Chan et al.
(2014) relate the $f$-mode frequency to the tidal deformability. For waveform
models that include the quadrupolar ($l=2$) and octopolar ($l=3$) adiabatic
and dynamical tidal effects Andersson and Pnigouras (2019), these universal
relations Yagi and Yunes (2013a, b); Yagi (2014); Chan et al. (2014) can be
used to reduce the intrinsic matter parameters from ten to two, namely the
$l=2$ (scaled) tidal deformabilities of the two stars. The other eight matter
parameters include the spin-induced quadrupole moments, the $l=3$ tidal
deformabilities, and the $l=2$ and $l=3$ $f$-mode frequencies of the two
stars.

In view of their relevance for neutron-star astrophysics and gravitational-
wave physics, it is important to test the robustness of these universal
relations. Although they have been demonstrated to be insensitive to many EOS
models to the $O(1\%)$ level, it is not yet well established whether these
universal relations are also insensitive to more realistic neutron star
physics aspects, like the state of matter at high densities. Work in this
direction includes the extension of the I-Q relation for rapidly rotating
stars. While it was originally found that the relation becomes more EOS-
dependent when considering rapidly rotating stars with a fixed rotation
frequency Doneva et al. (2013), it was later found that the I-Q relation, in
fact, remains approximately EOS-insensitive if one uses dimensionless
parameters to characterize the rotation Pappas and Apostolatos (2014);
Chakrabarti et al. (2014). The I-Q relation is also found to fail for slowly
rotating neutron stars with very strong magnetic fields so that the stellar
deformation is dominated by the magnetic field instead of rotation Haskell et
al. (2013). Thermal effects relevant for protoneutron stars have also been
shown to break the I-Love-Q relations Martinon et al. (2014); Marques et al.
(2017). It has also been demonstrated that the I-Love-Q relations are still
satisfied to within $\sim 3$% for hybrid star models with a strong first-order
hadron-quark phase transition in the interior Paschalidis et al. (2018).
However, the I-Love relation can be broken if the quark matter core is in a
crystalline color-superconducting state Lau et al. (2017, 2019). In this
paper, we add a contribution to this line of research by studying the impact
of nucleon superfluid dynamics on the I-Love-Q relations.

Apart from newly born neutron stars in supernova explosions and possibly also
binary neutron stars during the merger phase, typical neutron stars are very
cold on the nuclear temperature scale ($\sim 10^{10}$ K). The internal
temperature of a newborn neutron star is expected to drop quickly below the
transition temperatures ($\sim 10^{9}$ K) for neutrons and protons to become
superfluid and superconducting (see, e.g., Lombardo and Schulze (2001);
Haskell and Sedrakian (2018); Sedrakian and Clark (2019)). It is thus expected
that nucleon superfluidity will exist in mature neutron stars. In fact, the
pulsar glitch phenomenon is generally thought to be a manifestation of the
neutron superfluid. While the exact mechanism is still not well understood,
the basic idea involves the transfer of angular momentum from the neutron
superfluid to the normal component, leading to a sudden spin up. In the
standard model for large glitches (see Haskell and Melatos (2015) for a
review) like those observed in the Vela pulsar, the neutron superfluid rotates
by forming a dense array of vortices inside the neutron star. As the star
spins down due to electromagnetic radiation, the vortices are pinned to the
crust Anderson and Itoh (1975) and the neutron superfluid essentially
decouples from the normal fluid and does not spin down. A lag between the
superfluid and the normal component develops, and a Magnus force is induced on
the vortices. When the lag is large enough, the vortices will suddenly unpin
and the superfluid will spin down. The crust then spins up, leading to a
glitch, due to the conservation of angular momentum. Furthermore, the long
relaxation time following pulsar glitches is also suggested to be a signature
of the superfluid component Baym et al. (1969)

Since there are strong theoretical and observational motivations to suggest
its existence in neutron stars, it is important to understand if and how
nucleon superfluidity may leave a signature on the gravitational waves emitted
from a binary neutron star system. After all, the stars are still expected to
be cold when the system sweeps through the sensitivity bands of ground-based
detectors. It has been estimated that the core temperature of superfluid
neutron stars is only heated up to $\sim 10^{7}$ K due to tidal heating in
inspiralling binaries Yu and Weinberg (2016), significantly below the
superfluid transition temperature ($\sim 10^{9}$ K). However, numerical
simulations Sekiguchi et al. (2011); Bernuzzi et al. (2016); Perego et al.
(2019) suggest that the temperature of a post-merger hypermassive neutron star
can be as high as several tens of $10^{10}$ K, at which superfluidity is
expected to be destroyed. As discussed above, the I-Love-Q relations can be
used to reduce the number of matter parameters in theoretical waveform models.
More specifically, it is the Q-Love relation that is used to express the spin-
induced quadrupole moment by the tidal deformability in waveform models.
However, this is possible only if the relations are insensitive to different
EOS models and physics input. In this work, we test whether the I-Love-Q
relations remain valid for superfluid neutron stars. If the relations are
broken by the superfluid dynamics, we may then (in principle) be able to probe
the existence and properties of nucleon superfluidity in neutron stars by
comparing observational data against waveform models with and without the
assumption of the I-Love-Q relations.

Protons, electrons, and nuclei in neutron stars couple strongly (via the
electromagnetic interaction) on a very short timescale. When neutrons become
superfluid, they decouple from the charged components to a first
approximation. The interior of a superfluid neutron star can thus be
approximated by a two-fluid system: the neutron superfluid and the “proton”
fluid containing all charged particles. Besides being coupled via gravity, the
two fluids can also be coupled through the entrainment effect so that the flow
of one fluid induces a momentum in the other (as a result of the strong
interaction between neutrons and protons). The properties and dynamics of
superfluid neutron stars have been studied using this two-fluid model in a
general relativistic framework (see Andersson and Comer (2020) for a review).
In this paper, we use the two-fluid model to study the effects of superfluid
dynamics on the I-Love-Q relations. However, we neglect the effects of the
solid crust as it is expected to have only a tiny effect on global stellar
quantities like the tidal deformability Gittins et al. (2020).

The plan of this paper is as follows. In Section 2, we provide an outline of
the general relativistic two-fluid formalism that we employ to calculate the
moment of inertia, the spin-induced quadrupole moment, and the tidal
deformability of superfluid neutron stars. Section 3 describes the EOS models
used in this study. Section 4 presents our numerical results. Finally, we
summarize and discuss our results in Section 5. We use units where
$G=c=\hbar=1$ unless otherwise noted.

##  2 General Relativistic Two-Fluid Formalism

Our study is based on the general relativistic two-fluid formalism developed
by Carter and his collaborators (e.g., Carter (1989); Comer and Langlois
(1993, 1994); Carter and Langlois (1998); Langlois et al. (1998)). The
formalism is built around the master function $\Lambda(n^{2},p^{2},x^{2})$,
which is formed by the scalars constructed from the neutron $n^{\mu}$ and
proton $p^{\mu}$ number density currents: $n^{2}=-n_{\mu}n^{\mu}$,
$p^{2}=-p_{\mu}p^{\mu}$, and $x^{2}=-n_{\mu}p^{\mu}$. As already mentioned,
the “proton” fluid is used to refer to a conglomerate made of all charged
components. The master function plays the role of the EOS in the two-fluid
formulation. Given a master function $\Lambda$, the stress-energy tensor is
determined by

$T^{\mu}_{\nu}=\Psi\delta^{\mu}_{\nu}+n^{\mu}\mu_{\nu}+p^{\mu}\chi_{\nu},$ (1)

where the generalized pressure $\Psi$ is given by

$\Psi=\Lambda-n^{\alpha}\mu_{\alpha}-p^{\alpha}\chi_{\alpha}.$ (2)

The chemical potential covectors $\mu_{\alpha}$ and $\chi_{\alpha}$ (the fluid
four-momenta), respectively, for the neutrons and (conglomerate) protons are

$\mu_{\alpha}={\cal B}n_{\alpha}+{\cal A}p_{\alpha},\ \ \ \chi_{\alpha}={\cal
C}p_{\alpha}+{\cal A}n_{\alpha},$ (3)

where

${\cal A}=-{\partial\Lambda\over\partial x^{2}},\ \ \ {\cal
B}=-2{\partial\Lambda\over\partial n^{2}},\ \ \ {\cal
C}=-2{\partial\Lambda\over\partial p^{2}}.$ (4)

The coefficient $\cal A$ captures the entrainment effect Andersson and Comer
(2020); Andersson (2021) between the two fluids through which the current of
one fluid will induce a momentum in the other fluid. The equations of motion
for the two fluids consist of two conservation equations,

$\nabla_{\mu}n^{\mu}=0,\ \ \ \nabla_{\mu}p^{\mu}=0,$ (5)

and two Euler equations,

$n^{\mu}\nabla_{[\mu}\mu_{\nu]}=0,\ \ \ p^{\mu}\nabla_{[\mu}\chi_{\nu]}=0.$
(6)

The I-Love-Q relations connect the moment of inertia $I$ and quadrupole moment
$Q$ of slowly rotating stars to the tidal deformability $\lambda_{\rm tid}$ of
nonrotating stars. These quantities are determined by perturbative
calculations starting from a nonrotating equilibrium background solution. In
the following, we shall outline the main steps involved in deriving these
quantities and summarize the relevant equations for our discussion. We refer
the reader to Comer et al. (1999); Andersson and Comer (2001); Char and Datta
(2018) for more details.

###  2.1 Nonrotating Stars

The unperturbed background solution is assumed to be a spherically symmetric
and static spacetime described by the metric

$ds^{2}=-e^{\nu(r)}dt^{2}+e^{\lambda(r)}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2}).$ (7)

The equilibrium structure of a nonrotating superfluid neutron star, assuming
that the neutron and proton fluids coexist throughout the whole star, is
studied in Comer et al. (1999). We also make this assumption in our study.
This is obviously a simplified model as a realistic neutron star is expected
to be a multilayer system due to the physics of the superfluid phase
transition Andersson et al. (2002); Lin et al. (2008). For instance, due to
the density dependence of the superfluid gap function, the two-fluid region
may be sandwiched by single-fluid layers in a realistic neutron star (see
Figure 2 of Lin et al. (2008) for an illustration). However, we expect that
the multilayer aspects will reduce the two-fluid region in the star and make
the I-Love-Q relations less sensitive to the two-fluid dynamics. As we shall
see later, the I-Love-Q relations remain valid to high accuracy for our
simplified two-fluid stellar model unless the value of the spin ratio between
the two fluids becomes quite unrealistic. We expect that a more realistic
multilayer model would strengthen this conclusion.

The resulting metric and hydrostatic equilibrium equations are given by
Equations (25)–(27) of Comer et al. (1999):

$\displaystyle\lambda^{\prime}={1-e^{\lambda}\over r}-8\pi
re^{\lambda}\Lambda,$ (8) $\displaystyle\nu^{\prime}=-{1-e^{\lambda}\over
r}+8\pi re^{\lambda}\Psi,$ (9) $\displaystyle{\cal A}^{0}_{0}p^{\prime}+{\cal
B}^{0}_{0}n^{\prime}+{1\over 2}({\cal B}n+{\cal A}p)\nu^{\prime}=0,$ (10)
$\displaystyle{\cal C}^{0}_{0}p^{\prime}+{\cal A}^{0}_{0}n^{\prime}+{1\over
2}({\cal A}n+{\cal C}p)\nu^{\prime}=0,$ (11)

where primes denote derivatives with respect to $r$ and the coefficients
${\cal A}^{0}_{0}$, ${\cal B}^{0}_{0}$, and ${\cal C}^{0}_{0}$ are determined
from the master function $\Lambda$ (see Comer et al. (1999) for the explicit
expressions). The radius $R$ of the background star is defined by the
condition that the pressure vanishes at the surface (i.e., $\Psi(R)=0$) and
the total mass $M$ is obtained by

$M=-4\pi\int^{R}_{0}r^{2}\Lambda(r)dr.$ (12)

###  2.2 Slowly Rotating Stars

Rotating superfluid neutron stars have been studied using the general
relativistic two-fluid formalism Andersson and Comer (2001); Prix et al.
(2005); Sourie et al. (2016). In this work, we determine the moment of inertia
and quadrupole moment of slowly rotating superfluid neutron stars by following
the formalism developed in Andersson and Comer (2001). The spacetime for a
rotating star is given by an axisymmetric and stationary metric of the form

$ds^{2}=-\left[N^{2}-\sin^{2}\theta
K(N^{\phi})^{2}\right]dt^{2}+Vdr^{2}-2\sin^{2}\theta
KN^{\phi}dtd\phi+K\left(d\theta^{2}+\sin^{2}\theta d\phi^{2}\right).$ (13)

The neutron and proton fluids are assumed to be uniformly rotating with
angular velocities $\Omega_{\rm n}$ and $\Omega_{\rm p}$, respectively. In the
slow rotation approximation, the frequencies $\left(\Omega_{\rm n},\Omega_{\rm
p},\sqrt{\Omega_{\rm n}\Omega_{\rm p}}\right)$ are assumed to be much smaller
than the characteristic Kepler frequency, which is proportional to
$\sqrt{M/R^{3}}$ in Newtonian theory, and one can then expand the metric
functions up to second order in the rotation as

$\displaystyle N=e^{\nu(r)/2}(1+h_{0}(r)+h_{2}(r)P_{2}(\cos\theta)),$ (14)
$\displaystyle V=e^{\lambda(r)}(1+2v_{0}(r)+2v_{2}(r)P_{2}(\cos\theta)),$ (15)
$\displaystyle K=r^{2}(1+2k_{2}(r)P_{2}(\cos\theta)),$ (16) $\displaystyle
N^{\phi}=\omega(r),$ (17)

where $P_{2}(\cos\theta)=(3\cos^{2}\theta-1)/2$ is the second-order Legendre
polynomial. The perturbed metric functions $h_{0}$, $h_{2}$, $v_{0}$, $v_{2}$,
and $k_{2}$ are second order in angular velocities. The frame-dragging term
$\omega(r)$ is a first-order quantity determined by

${1\over r^{4}}\left(r^{4}e^{-(\lambda+\nu)/2}{\tilde{L}}^{\prime}_{\rm
n}\right)^{\prime}-16\pi e^{(\lambda-\nu)/2}(\Psi-\Lambda){\tilde{L}}_{\rm
n}=16\pi e^{(\lambda-\nu)/2}\chi p(\Omega_{\rm n}-\Omega_{\rm p}),$ (18)

where ${\tilde{L}}_{\rm n}=\omega-\Omega_{\rm n}$ and ${\tilde{L}}_{\rm
p}=\omega-\Omega_{\rm p}$. It should be noted that the equation reduces to the
corresponding frame-dragging equation for a single-fluid model Hartle (1967)
when the two fluids are co-rotating (i.e., $\Omega_{\rm n}=\Omega_{\rm p}$).
Equation (18) is solved by choosing a central value ${\tilde{L}}_{\rm n}(0)$
so that the interior and exterior solutions match at the stellar surface to
satisfy the condition

${\tilde{L}}_{\rm n}(R)=-\Omega_{\rm n}+{2J\over R^{3}},$ (19)

where $J=J_{\rm n}+J_{\rm p}$ is the total angular momentum of the star. The
neutron and proton angular momenta are given, respectively, by

$\displaystyle J_{\rm n}=-{8\pi\over
3}\int^{R}_{0}drr^{4}e^{(\lambda-\nu)/2}\left[\mu n{\tilde{L}}_{\rm n}+{\cal
A}np(\Omega_{\rm n}-\Omega_{\rm p})\right],$ (20) $\displaystyle J_{\rm
p}=-{8\pi\over 3}\int^{R}_{0}drr^{4}e^{(\lambda-\nu)/2}\left[\chi
p{\tilde{L}}_{\rm p}+{\cal A}np(\Omega_{\rm p}-\Omega_{\rm n})\right].$ (21)

The moment of inertia of a single-fluid star in general relativity is defined
by the ratio between the angular momentum and angular velocity of the star. We
generalize the definition to a two-fluid star so that the moments of inertia
of the neutron and proton fluids are defined, respectively, by

$I_{\rm n}={J_{\rm n}\over\Omega_{\rm n}},\ \ \ I_{\rm p}={J_{\rm
p}\over\Omega_{\rm p}}.$ (22)

The total moment of inertia is then given by $I=I_{\rm n}+I_{\rm p}$, which
reduces to the definition for a single-fluid star if the two fluids are co-
rotating (see also Sourie et al. (2016)). This completes the discussion of a
slowly rotating two-fluid star up to first order in the angular velocities.

In order to study the spin-induced quadrupole moment $Q$ of a slowly rotating
star, one needs to consider the metric and fluid perturbations up to second
order in the angular velocities. The second-order metric perturbation
functions defined in Equation (17) consist of the $l=0$ ($h_{0},v_{0}$) and
$l=2$ ($h_{2},v_{2},k_{2}$) contributions, where $l$ is the order of the
Legendre polynomials $P_{l}(\cos\theta)$. The equations and numerical schemes
for determining the interior solutions of these quantities can be found in
Andersson and Comer (2001). Outside the star, the problem is identical to the
single-fluid case and the metric functions can be obtained analytically Hartle
(1967). In particular, the exterior solution for $h_{2}$ is given by

$\displaystyle h_{2}(r)$ $\displaystyle=$ $\displaystyle-A\left[{3\over
2}\left({r\over M}\right)^{2}\left(1-{2M\over r}\right)\ln\left(1-{2M\over
r}\right)+{(r-M)(3-6M/r-2(M/r)^{2})\over M(1-2M/r)}\right]$ (24)
$\displaystyle+{J^{2}\over Mr^{3}}\left(1+{M\over r}\right),$

where the constant $A$ is determined by matching the interior and exterior
solutions at the stellar surface. Far from the star, the function
$h_{2}(r)P_{2}(\cos\theta)$ becomes the perturbed Newtonian potential and the
quadrupole moment can be read off from the coefficient of the
$P_{2}(\cos\theta)/r^{3}$ term:

$Q=-{8\over 5}AM^{3}-{J^{2}\over M}.$ (25)

The normalized (dimensionless) moment of inertia $\bar{I}$ and quadrupole
moment $\bar{Q}$ that appear in the I-Love-Q relations are defined by

${\bar{I}}={I\over M^{3}},\ \ \ {\bar{Q}}=-{Q\over a^{2}M^{3}},$ (26)

where $a=J/M^{2}$ is the dimensionless spin parameter.

###  2.3 Tidally Deformed Nonrotating Stars

The tidal deformability $\lambda_{\rm tid}$ measures the deformation of a
neutron star due to the tidal field produced by the companion in a binary
system. In the static-tide limit and when the separation between the two stars
is large compared to the radii of the stars, the computation of $\lambda_{\rm
tid}$ for nonrotating single-fluid neutron stars is well established Hinderer
(2008); Damour and Nagar (2009). The tidal deformations of slowly rotating
single-fluid neutron stars have also been studied in Pani et al. (2015);
Landry (2017). The formulation has recently been extended to nonrotating
superfluid neutron stars within the general relativistic two-fluid formalism
Char and Datta (2018); Datta and Char (2020). Similar to the study of slowly
rotating stars discussed in Section 2.2, the tidal deformability is determined
by perturbing a nonrotating background solution. We focus on the dominant
quadrupolar ($l=2$) static tidal field and consider the even-parity
perturbations in the Regge–Wheeler gauge so that the full metric is given by

$\displaystyle ds^{2}$ $\displaystyle=$
$\displaystyle-e^{\nu(r)}\left[1+H_{0}(r)P_{2}(\cos\theta)\right]dt^{2}+e^{\lambda(r)}\left[1+H_{2}(r)P_{2}(\cos\theta)\right]dr^{2}$
(28)
$\displaystyle+r^{2}\left[1+K(r)P_{2}(\cos\theta)\right](d\theta^{2}+\sin^{2}\theta
d\phi^{2}).$

The perturbed Einstein equations impose the condition $H_{0}=-H_{2}\equiv H$.
The linearized metric and fluid equations inside the star yield the following
equation for determining the tidal deformability Char and Datta (2018)

$H^{{}^{\prime\prime}}+H^{{}^{\prime}}\left\\{{2\over
r}+e^{\lambda}\left[{2m\over r^{2}}+4\pi
r(\Psi+\Lambda)\right]\right\\}+H\left[-{6e^{\lambda}\over r^{2}}+4\pi
e^{\lambda}(9\Psi-5\Lambda-g)-\nu^{{}^{\prime}2}\right]=0,$ (29)

where the function $m$ is defined by $e^{\lambda(r)}=(1-2m(r)/r)^{-1}$ and

$g={{\mu^{2}{\cal C}^{0}_{0}+\chi^{2}{\cal B}^{0}_{0}-2\mu\chi{\cal
A}^{0}_{0}}\over({\cal A}^{0}_{0})^{2}-{\cal B}^{0}_{0}{\cal C}^{0}_{0}}.$
(30)

Noticing that $\Psi$ is the pressure and $\Lambda$ is the negative energy
density of the unperturbed nonrotating background star, one can compare
Equation (29) with the corresponding equation in a single-fluid situation (see
Equation (15) of Hinderer (2008)) for a barotropic EOS $P(\rho)$ and observe
that the term $(\rho+P)/(dP/d\rho)$ in the single-fluid equation is now
replaced by the function $-g$ in the two-fluid equation. It can also be shown
that Equation (29) reduces to Equation (15) of Hinderer (2008) when the master
function depends only on one particle density (e.g., $\Lambda(p^{2})$) in the
single-fluid limit. Once the interior problem (Equation (29)) is solved for
the metric perturbation function $H$, the remaining step to determine the
tidal deformability is the same as that for the single-fluid counterpart.

The response of a nonrotating star to an external quadrupolar tidal field
${\cal E}_{ij}$ results in the following expansion of the metric function
$g_{tt}$ in the star’s local asymptotic rest frame Hinderer (2008):

$-{1+g_{tt}\over 2}=-{M\over r}-{3Q_{ij}\over 2r^{3}}\left({x^{i}x^{j}\over
r^{2}}-{1\over 3}\delta^{ij}\right)+{1\over 2}{\cal E}_{ij}x^{i}x^{j},$ (31)

where $Q_{ij}$ is the tidally induced traceless quadrupole moment tensor of
the star. The tidal deformability $\lambda_{\rm tid}$ of the star is defined
by $Q_{ij}=-\lambda_{\rm tid}{\cal E}_{ij}$. It is also convenient to define
the dimensionless ($l=2$) tidal Love number $k_{2}\equiv 3\lambda_{\rm
tid}/2R^{5}$. Outside the star, the solutions of Equation (29) are given by
the associated Legendre functions $Q_{2}^{2}(x)$ and $P_{2}^{2}(x)$:

$H(r)=c_{1}Q_{l}^{2}(x)+c_{2}P_{l}^{2}(x),$ (32)

where $x=r/M-1$. Using Equation (31), the coefficients $c_{1}$ and $c_{2}$ are
related to the tidal Love number by $k_{2}=4c_{1}M^{5}/15c_{2}R^{5}$. The
value of $c_{1}/c_{2}$ can be obtained by matching the interior and exterior
solutions of Equation (29) at the surface. The tidal Love number can then be
expressed as Hinderer (2008)

$\displaystyle k_{2}$ $\displaystyle=$ $\displaystyle{8\over
5}C^{5}(1-2C)^{2}\left[2+2C\left(y-1\right)-y\right]\Bigg{\\{}2C\left[4(1+y)C^{4}+(6y-4)C^{3}+(26-22y)C^{2}\right.$
(33)
$\displaystyle\left.+3C(5y-8)-3y+6\right]+3(1-2C)^{2}\left[2-y+2C(y-1)\right]\text{log}(1-2C)\Bigg{\\}}^{-1},$

where $C=M/R$ is the compactness and $y=RH^{{}^{\prime}}(R)/H(R)$. Finally,
the normalized (dimensionless) tidal deformability ${\bar{\lambda}}_{\rm tid}$
that appears in the I-Love-Q relations is defined by

${\bar{\lambda}}_{\rm tid}={\lambda_{\rm tid}\over M^{5}}={2\over
3}k_{2}C^{-5}.$ (34)

##  3 Equations Of State

###  3.1 Two-Fluid Polytropic Model

As discussed in the previous section, the master function
$\Lambda(n^{2},p^{2},x^{2})$ plays the role of EOS in the two-fluid
formulation. We will use two different EOS models in this work. The master
function of the first model is taken to be Comer et al. (1999)

$\Lambda=-m_{\rm n}n-\sigma_{\rm n}n^{\beta_{\rm n}}-m_{\rm p}p-\sigma_{\rm
p}p^{\beta_{\rm p}},$ (35)

where the neutron and proton masses are assumed to be equal ($m_{\rm n}=m_{\rm
p}$) for simplicity. We choose the parameters $\sigma_{\rm n}=0.2m_{\rm n}$,
$\sigma_{\rm p}=2m_{\rm n}$, $\beta_{\rm n}=2.3$, and $\beta_{\rm p}=1.95$ as
in Andersson and Comer (2001) to construct background equilibrium models in
this work. As the master function does not depend on $x^{2}$, there is no
entrainment effect and the fluids behave as independent polytropes and couple
only through gravity. The various coefficients ${\cal A},{\cal B}$, etc., (see
Equation (11)) can be easily computed for this model Comer et al. (1999). In
particular, we have ${\cal A}={\cal A}^{0}_{0}=0$.

###  3.2 Relativistic Mean Field Model

The second EOS that we will consider is a more realistic model based on
relativistic mean field (RMF) approximation. RMF based models have been used
previously in the study of general relativistic superfluid neutron stars.
Comer and Joynt Comer and Joynt (2003); Comer (2004) have calculated the
entrainment effect in a relativistic $\sigma-\omega$ model. Kheto and
Bandyopadhyay Kheto and Bandyopadhyay (2014) later extended the study to
include the isospin dependence by using a relativistic $\sigma-\omega-\rho$
model. More recently, Char and Datta Char and Datta (2018); Datta and Char
(2020) used the same RMF model to study the tidal deformability of superfluid
neutron stars. In this work, we also use the relativistic $\sigma-\omega-\rho$
model to construct a master function. The Lagrangian density for the system is
given by

$\displaystyle{\cal L}$ $\displaystyle=$ $\displaystyle\sum_{B={\rm
n,p}}{\bar{\Psi}}_{B}\left(i\gamma_{\mu}\partial^{\mu}-m_{B}+g_{\sigma{B}}\sigma-
g_{\omega{B}}\gamma_{\mu}\omega^{\mu}-g_{\rho{B}}\gamma_{\mu}\boldsymbol{\tau}_{B}\cdot\boldsymbol{\rho}^{\mu}\right){\Psi}_{B}$
(38) $\displaystyle-{1\over 2}\partial_{\mu}\sigma\partial^{\mu}\sigma-{1\over
2}m^{2}_{\sigma}\sigma^{2}-{1\over 3}bm(g_{\sigma}\sigma)^{3}-{1\over
4}c(g_{\sigma}\sigma)^{4}-{1\over 4}\omega_{\mu\nu}\omega^{\mu\nu}$
$\displaystyle-{1\over 2}m^{2}_{\omega}\omega_{\mu}\omega^{\mu}-{1\over
4}\boldsymbol{\rho}_{\mu\nu}\cdot\boldsymbol{\rho}^{\mu\nu}-{1\over
2}m^{2}_{\rho}\boldsymbol{\rho}_{\mu}\cdot\boldsymbol{\rho}^{\mu},$

where $\Psi_{B}$ is the Dirac spinor for baryons ${B}$ with baryon mass
$m_{B}$ and $\gamma_{\mu}$ denotes the Dirac matrices. The isospin operator is
denoted by $\boldsymbol{\tau}_{B}$. In this theory, we have the self-
interacting scalar $\sigma$ field, the vector omega field $\omega^{\mu}$, and
the isovector rho field $\boldsymbol{\rho}^{\mu}$. The latter two fields have
the corresponding field tensors $\omega_{\mu\nu}$ and
$\boldsymbol{\rho}_{\mu\nu}$. The nucleon mass $m$ is taken to be the average
of the bare neutron and proton masses. The resulting equations of motion from
the Lagrangian are solved in the RMF approximation. The entrainment effect is
incorporated by choosing a frame in which the neutrons have zero spatial
momentum and the proton momentum is given by $k_{\rm p}^{\mu}=(k_{0},0,0,K)$
Comer and Joynt (2003). In the limit $K\rightarrow 0$ that is relevant to the
slow rotation approximation considered in this work, the master function for
determining the equilibrium background neutron stars is given explicitly by
Char and Datta (2018)

$\displaystyle\Lambda$ $\displaystyle=$ $\displaystyle-{c_{\omega}^{2}\over
18\pi^{4}}(k_{\rm n}^{3}+k_{\rm p}^{3})^{2}-{c_{\rho}^{2}\over
72\pi^{4}}(k_{\rm p}^{3}-k_{\rm n}^{3})^{2}-{1\over 4\pi^{2}}\left(k_{\rm
n}^{3}\sqrt{k_{\rm n}^{2}+m_{*}^{2}|_{0}}+k_{\rm p}^{3}\sqrt{k_{\rm
p}^{2}+m_{*}^{2}|_{0}}\right)$ (42) $\displaystyle-{1\over
4c_{\sigma}^{2}}\Big{\\{}(2m-m_{*}|_{0})(m-m_{*}|_{0})+m_{*}|_{0}\left[bmc_{\sigma}^{2}(m-m_{*}|_{0})^{2}+cc_{\sigma}^{2}(m-m_{*}|_{0})^{3}\right]\Big{\\}}$
$\displaystyle-{1\over 3}bm(m-m_{*}|_{0})^{3}-{1\over
4}c(m-m_{*}|_{0})^{4}-{1\over 8\pi^{2}}\Bigg{\\{}k_{\rm p}(2k_{\rm
p}^{2}+m_{\rm e}^{2})\sqrt{k_{\rm p}^{2}+m_{\rm e}^{2}}$ $\displaystyle-m_{\rm
e}^{4}\ln\left({k_{\rm p}+\sqrt{k_{\rm p}^{2}+m_{\rm e}^{2}}\over m_{\rm
e}}\right)\Bigg{\\}},$

where $k_{\rm n}=(3\pi^{2}n)^{1/3}$ and $k_{\rm p}=(3\pi^{2}p)^{1/3}$ are the
neutron and proton Fermi momenta, respectively. We have also added the
contributions due to the electrons in the master function (the last term
containing the electron mass $m_{\rm e}$). The background value for the Dirac
effective mass $m_{*}|_{0}$ is determined by the following transcendental
equation

$\displaystyle m_{*}|_{0}$ $\displaystyle=$ $\displaystyle
m-m_{*}|_{0}{c_{\sigma}^{2}\over 2\pi^{2}}\Bigg{\\{}k_{\rm n}\sqrt{k_{\rm
n}^{2}+m^{2}_{*}|_{0}}+k_{\rm p}\sqrt{k_{\rm p}^{2}+m^{2}_{*}|_{0}}$ (45)
$\displaystyle+{1\over 2}m^{2}_{*}|_{0}\ln\left[{-k_{\rm n}+\sqrt{k_{\rm
n}^{2}+m^{2}_{*}|_{0}}\over k_{\rm n}+\sqrt{k_{\rm
n}^{2}+m^{2}_{*}|_{0}}}\right]+{1\over 2}m^{2}_{*}|_{0}\ln\left[{-k_{\rm
p}+\sqrt{k_{\rm p}^{2}+m^{2}_{*}|_{0}}\over k_{\rm p}+\sqrt{k_{\rm
p}^{2}+m^{2}_{*}|_{0}}}\right]\Bigg{\\}}$
$\displaystyle+bmc_{\sigma}^{2}(m-m_{*}|_{0})^{2}+cc_{\sigma}^{2}(m-m_{*}|_{0})^{3}.$

We refer the reader to Char and Datta (2018) for the expressions of the
various coefficients $\cal A$, $\cal B$, etc. (see Equation (11)). In contrast
to the two-fluid polytropic model, it should be noted that the coefficient
$\cal A$, which is responsible for the entrainment effect, is nonzero in this
RMF model. The coupling constants
$c_{\sigma}^{2}\equiv(g_{\sigma}/m_{\sigma})^{2}$,
$c_{\omega}^{2}\equiv(g_{\omega}/m_{\omega})^{2}$,
$c_{\rho}^{2}\equiv(g_{\rho}/m_{\rho})^{2}$, $b$, and $c$ in the master
function are determined by the nuclear matter saturation properties Kheto and
Bandyopadhyay (2014). In the following, we shall use the same NL3 Fattoyev et
al. (2010) and GM1 Glendenning and Moszkowski (1991) parameter sets (see Table
3.2) as employed in Char and Datta (2018).

Parameter sets for the NL3 and GM1 models. The nucleon-meson coupling
constants $c_{\sigma}^{2}$, $c_{\omega}^{2}$, and $c_{\rho}^{2}$ are expressed
in ${\rm fm}^{2}$, while $b$ and $c$ are dimensionless Char and Datta (2018).
Model $c_{\sigma}^{2}$ $c_{\omega}^{2}$ $c_{\rho}^{2}$ $b$ $c$ NL3 15.739
10.530 5.324 0.002055 $-0.002650$ GM1 11.785 7.148 4.410 0.002948 $-0.001071$

Before studying the I-Love-Q relations, it is worth pointing out that our two-
fluid stellar models are dominated by the neutron superfluid as the typical
central value of the proton fraction in our models is about 10%, which is
comparable to that of a realistic neutron star. For example, Figure 1 shows
the profiles of the proton fraction $p/(n+p)$ for $1.4M_{\odot}$ neutron stars
constructed from our EOS models.

Figure 1: Profiles of the proton fraction $p/(n+p)$ for $1.4M_{\odot}$ neutron
stars constructed from the polytropic (Poly), GM1, and NL3 EOS models.

##  4 Numerical Results

The I-Love-Q universal relations Yagi and Yunes (2013a, b) originally found
for single-fluid neutron stars connect the three dimensionless quantities
$\bar{I}$, $\bar{Q}$, and ${\bar{\lambda}}_{\rm tid}$. While one might naively
expect that $\bar{I}$ and $\bar{Q}$ are somehow related to each other as they
both characterize the effects of rotation, the fact that these quantities
together with ${\bar{\lambda}}_{\rm tid}$ are connected by the following
approximately EOS-independent relations are quite surprising Yagi and Yunes
(2013b):

$\ln y_{i}=a_{i}+b_{i}\ln x_{i}+c_{i}(\ln x_{i})^{2}+d_{i}(\ln
x_{i})^{3}+e_{i}(\ln x_{i})^{4},$ (46)

where $(x_{i},y_{i})$ are any two of $\bar{I}$, $\bar{Q}$, and
${\bar{\lambda}}_{\rm tid}$. The set of fitting coefficients
$(a_{i},b_{i},c_{i},d_{i},e_{i})$ are different for different pairs of $x_{i}$
and $y_{i}$ (see Table I in Yagi and Yunes (2013b)).

Our aim in this work is to study whether the I-Love-Q relations remain valid
for superfluid neutron stars. We have generalized the definitions of
$\bar{I}$, $\bar{Q}$, and ${\bar{\lambda}}_{\rm tid}$ to two-fluid stars in
Section 2. For a given master function (i.e., a EOS model) in the two-fluid
formulation, the procedure in our calculation is to first build a nonrotating
background star for given central number densities $n(0)$ and $p(0)$. Note
that the two densities are not independent as the background star is assumed
to be in chemical equilibrium Comer et al. (1999). A slowly rotating model is
then built by solving the perturbative equations with the nonrotating
background model as input. The neutron and proton fluids are assumed to be
rigidly rotating with angular velocities $\Omega_{\rm n}$ and $\Omega_{\rm
p}$. However, instead of providing $\Omega_{\rm n}$ and $\Omega_{\rm p}$ as
input parameters, the system of equations can be scaled in such a way that it
is sufficient to specify the spin ratio $\Omega_{\rm n}/\Omega_{\rm p}$ in the
calculation Andersson and Comer (2001). We will thus present our results for
stellar sequences specified by $\Omega_{\rm n}/\Omega_{\rm p}$. The spin ratio
is expected to be very close to unity as suggested by the pulsar glitch
phenomenon (see Section 5 for discussion). However, we will explore a much
wider parameter range to study the validity of the I-Love-Q relations. The
dimensionless moment of inertia $\bar{I}$ and quadrupole moment $\bar{Q}$ can
be determined for the slowly rotating star as outlined in Section 2.2. The
computation of the dimensionless tidal deformability ${\bar{\lambda}}_{\rm
tid}$ for the nonrotating star was also discussed in Section 2.3.

We first consider the results for the two-fluid polytropic EOS model. In
Figure 2a, we plot in the upper panel $\ln{\bar{I}}$ against
$\ln{\bar{\lambda}}_{\rm tid}$ (i.e., the I-Love relation) for sequences of
stars with spin ratios $\Omega_{\rm n}/\Omega_{\rm p}=0.4$, 0.7, 1, 1.3, and
1.6. The solid line is the fitting curve (Equation (46)) for the I-Love
relation for single-fluid ordinary neutron stars Yagi and Yunes (2013b). The
lower panel shows the relative error, $E=({\hat{y}}-y)/y$, between the
numerical data, $\hat{y}$, and the fitting curve, $y$. We see from the upper
panel that the numerical data can still match the fitting curve for single-
fluid stars very well. However, the effects of two-fluid dynamics become more
apparent in the lower panel. The co-rotating case $\Omega_{\rm n}/\Omega_{\rm
p}=1$ corresponds effectively to a single-fluid system and has the smallest
$|E|$ among the different sequences, as expected. It can also be seen that
$|E|$ becomes larger as the spin ratio deviates more away from unity. For the
cases where the neutron fluid rotates faster than the proton fluid (i.e.,
$\Omega_{\rm n}/\Omega_{\rm p}>1$), the error $|E|$ increases with
$\Omega_{\rm n}/\Omega_{\rm p}$. On the other hand, $|E|$ increases as
$\Omega_{\rm n}/\Omega_{\rm p}$ decreases from 0.7 to 0.4 for the opposite
situation where the proton fluid rotates faster. Since ${\bar{\lambda}}_{\rm
tid}$ decreases with increasing compactness, the numerical results thus
illustrate that the error $|E|$ increases with the compactness for a given
sequence.

((a))

((b))

((c))

Figure 2: The accuracy of the I-Love (a), Q-Love (b), and I-Q (c) relations
for our chosen two-fluid polytropic EOS using different spin ratios
$\Omega_{\rm n}/\Omega_{\rm p}$. The solid line in the upper panel of each
figure represents the corresponding fitting curve for ordinary single-fluid
neutron stars Yagi and Yunes (2013b). The lower panel of each figure shows the
relative error between the two-fluid results and the corresponding fitting
curve.

2

In Figure 2b,c, we plot $\ln{\bar{Q}}$ against $\ln{\bar{\lambda}}_{\rm tid}$
and $\ln{\bar{I}}$ against $\ln{\bar{Q}}$, respectively. Similar to Figure 2a,
the solid lines in these figures are the corresponding fitting curves for the
Q-Love and I-Q relations for single-fluid stars Yagi and Yunes (2013b). The
relative errors between the numerical data and the fitting curves are also
plotted in the lower panel of the figures. Similar to the I-Love relation in
Figure 2a, the errors $|E|$ for these two relations also increase as the spin
ratio deviates from unity. It should be pointed out that the I-Love-Q
relations (i.e., the solid lines) for single-fluid stars are approximately EOS
independent to within about the 1% level. This means that the I-Love-Q
relations may be considered to be broken by the effects of two-fluid dynamics
only for cases where $|E|>0.01$. For our chosen values of $\Omega_{\rm
n}/\Omega_{\rm p}$, we see that the I-Love-Q relations are broken only for the
case $\Omega_{\rm n}/\Omega_{\rm p}=0.4$. Taking our finding at face value, it
implies that accurate independent measurements of any pair of quantities in
the I-Love-Q relations can in principle be used to test the existence of two-
fluid dynamics in superfluid neutron stars, though quite extreme values of
$\Omega_{\rm n}/\Omega_{\rm p}$ are needed. The Q-Love relation plotted in
Figure 2b is the most promising in this aspect as the deviations between the
two-fluid data and the fitting curve are more significant than the other two
relations. It should be noted that the error for the Q-Love relation can reach
up to $|E|\sim 0.1$ level for the case $\Omega_{\rm n}/\Omega_{\rm p}=0.4$.

As discussed in Section 1, we can make use of various universal relations to
reduce the matter parameters in theoretical waveform models. In particular,
the Q-Love relation can be used to get rid of $\bar{Q}$ by expressing it in
terms of ${\bar{\lambda}}_{\rm tid}$, assuming that the relation is robust and
insensitive to EOS and various physics inputs to high accuracy. If the effects
of two-fluid dynamics can lead to $O(10\%)$ deviation of the Q-Love relation
as illustrated in our results, then this relation should be used with care in
waveform modeling. This also suggests that not imposing the Q-Love relation in
waveform modeling can, in principle, allow the tests of superfluid dynamics
using gravitational wave observations. These implications depend clearly on
whether the spin ratio $\Omega_{\rm n}/\Omega_{\rm p}$ can deviate
significantly from unity, as considered by us in Figure 2a–c. We will return
to this issue and assess the astrophysical relevance of our results in Section
5.

((a))

((b))

((c))

Figure 3: Similar to Figure 2, but for the GM1 and NL3 EOS models.

2

While we have only presented the results of one particular set of parameters
for the two-fluid polytropic EOS model in Figure 2a–c, we have checked that
the results obtained from different sets of EOS parameters are qualitatively
the same. In particular, the I-Love-Q relations remain valid to good accuracy
unless the spin ratio between the two fluids deviates significantly from
unity. However, the two-fluid polytropic model can at best only provide a
crude approximation to the properties of superfluid neutron stars since the
two fluids can couple only through gravity in this model. It is not obvious
that the conclusions derived from this simple model can be generalized
directly to more realistic EOS models. We fill this gap by considering the RMF
model. This model is more realistic in the sense that nucleon–nucleon
interactions are taken into account through the exchange of effective mesons
and the coupling parameters are determined by fitting to known nuclear matter
properties. The RMF model also has the advantage that the corresponding master
function is simple enough that the various thermodynamics coefficients that we
need in the two-fluid calculations can be obtained analytically. Most
importantly, in contrast to the polytropic EOS, the RMF model contains the
entrainment effect characterized by ${\cal A}\neq 0$, which is a unique
property of superfluid dynamics. In Figure 3a–c, we present the numerical
results for the NL3 and GM1 parametrizations of the RMF model (see Table 3.2).
In each figure, we consider three sequences of $\Omega_{\rm n}/\Omega_{\rm
p}=0.4$, 1, and 1.6 for each EOS. We see that the results of NL3 and GM1 EOSs
match the I-Love-Q relations to high accuracy for the co-rotating case
$\Omega_{\rm n}/\Omega_{\rm p}=1$ as expected. The data trends are
qualitatively the same as those of the polytropic EOS. We note that the
relative errors $|E|$ for all three relations increase as the spin ratio
deviates away from unity. For a given sequence of fixed EOS and spin ratio,
the errors also generally increase with the compactness. It can also be seen
from Figure 3b that the error for the Q-Love relation is the largest among the
three relations and can reach up to the $|E|\sim 0.1$ level.

Figure 4: The breaking of the Q-Love relation for extreme spin ratios as
illustrated by the polytropic (Poly), GM1, and NL3 EOS models. The numerical
values inside the parentheses in the legend box denote the values of
$\Omega_{\rm n}/\Omega_{\rm p}$. The solid line is the universal Q-Love
relation for ordinary single-fluid neutron stars Yagi and Yunes (2013b).

Let us now consider more extreme situations to explore the breaking of the
Q-Love relation when $\Omega_{\rm n}/\Omega_{\rm p}$ deviates significantly
from unity. We have studied the trend of the numerical results in the
$\ln{\bar{Q}}-\ln{\bar{\lambda}}_{\rm tid}$ plane. In Figure 4, we plot the
results for a few representative spin ratios as an illustration. The solid
line in the figure represents the original Q-Love relation for single-fluid
stars. We now consider how the data behave as the spin ratio deviates from
unity along the two opposite directions $\Omega_{\rm n}/\Omega_{\rm p}>1$ and
$<1$. In the former case, we find that the results converge to a large spin-
ratio limit, which can be approximated well by the case $\Omega_{\rm
n}/\Omega_{\rm p}=10$ (red data points). Increasing the spin ratio further
(e.g., $\Omega_{\rm n}/\Omega_{\rm p}=100$) has little effect on the results.
In this limit, the fluid motion is dominated by the neutron superfluid as the
proton fraction is small (see Figure 1). The possible range of deviation from
the Q-Love relation for $\Omega_{\rm n}/\Omega_{\rm p}>1$ is thus
approximately bounded by the solid line and the red data points. Note,
however, that the results for a given spin ratio are still EOS-insensitive in
this range.

We now turn to the situation for $\Omega_{\rm n}/\Omega_{\rm p}<1$, which
turns out to be more interesting. As we decrease the spin ratio from unity,
the results initially deviate from the solid line gradually. As we have seen
in Figure 3b, the results for $\Omega_{\rm n}/\Omega_{\rm p}=0.4$ are still
quite close to the solid line. However, as the spin ratio is decreased
further, the results deviate from the solid line more significantly. We also
see that the results become more sensitive to the EOS models as the data shift
further away from the solid line. Consider for example the results for
$\Omega_{\rm n}/\Omega_{\rm p}=0.1$ (black data points) in the figure. If we
keep decreasing the spin ratio toward zero, the data move further upward in
the figure and depend more on the EOS models. However, if the spin ratio is
further decreased to become negative so that the two fluids are counter-
rotating, the data would then shift downward in the figure. The results for
$\Omega_{\rm n}/\Omega_{\rm p}=-0.4$ and $-1$ in the figure show that the
results move closer to the solid line and become less sensitive to the EOS
models as the spin ratio becomes more negative. However, the results would not
converge to the solid line as the spin ratio is decreased further. Instead,
the data converge to the large spin-ratio limit that we discussed above for
$\Omega_{\rm n}/\Omega_{\rm p}>0$. It can be seen from the results of
$\Omega_{\rm n}/\Omega_{\rm p}=-10$ (green data points) that the data are
insensitive to the EOS models and agree quite well with the results of
$\Omega_{\rm n}/\Omega_{\rm p}=10$. In fact, this should be expected as the
system is dominated by the neutron superfluid if $|\Omega_{\rm n}/\Omega_{\rm
p}|\gg 1$ and the spin-induced quadrupole moment $Q$ should be independent of
the sense of rotation of the neutron superfluid. There should thus be a single
large spin-ratio limit.

Figure 4 shows that the results become more sensitive to the EOS models as the
data are far away from the single-fluid universal relation. While it is beyond
the scope of this work to study the EOS dependence for these extreme cases in
general, we have used different parameters for the two-fluid polytropic model
to gain some understanding of the general trend. In Table 4, we compare the
values of $\ln{\bar{Q}}$ at $\ln{\bar{\lambda}}_{\rm tid}=5$ for three
different polytropic EOS models. The EOS parameters are chosen in such a way
that the central proton fractions for $1.4M_{\odot}$ neutron stars are 0.02
(Poly_2), 0.04 (Poly_4), and 0.08 (Poly_8). For comparison, our default
polytropic model used in Figure 4 has a central proton fraction about 0.09 for
a $1.4M_{\odot}$ star (see Figure 1). In the table, the numerical values
inside the parentheses are the percentage differences between the data and the
value of $\ln{\bar{Q}}$ predicted by the single-fluid Q-Love relation. For a
given EOS model, we see that the percentage differences increase significantly
as the spin ratio decreases from 0.2 to 0.1. This is consistent with the above
observation that the deviations of the superfluid data from the I-Love-Q
relations increase as the spin ratio deviates further away from unity. On the
other hand, the deviations decrease with the proton fraction inside the stars
for a fixed spin ratio. This agrees with the expectation that the effects of
two-fluid dynamics should become less important if the system is dominated by
one of the fluids.

[H] Values of $\ln{\bar{Q}}$ at $\ln{\bar{\lambda}}_{\rm tid}=5$ for three
different polytropic EOS models with spin ratios 0.1 and 0.2. The EOS
parameters are chosen in such a way that the central proton fractions for
$1.4M_{\odot}$ stellar models are 0.02 (Poly_2), 0.04 (Poly_4), and 0.08
(Poly_8). The numerical values inside the parentheses are the percentage
differences between the data and the predictions of the single-fluid Q-Love
relation.  $\Omega_{\rm n}/\Omega_{\rm p}$ Poly_2 Poly_4 Poly_8 0.1 2.22 (59%)
2.46 (76%) 2.66 (91%) 0.2 1.66 (18%) 1.79 (28%) 1.96 (40%)

The results presented in Figure 4 may have little astrophysical relevance as
$|\Omega_{\rm n}/\Omega_{\rm p}|$ is not expected to deviate from unity
significantly. It is also uncertain how the two fluids can sustain counter-
rotation globally so that the condition $\Omega_{\rm n}/\Omega_{\rm p}<0$ can
be achieved in neutron stars. Nevertheless, the breaking of the Q-Love
relation in these extreme situations is interesting theoretically as it is
still not properly understood why various universal relations can exist at
all, though some suggestions have been proposed Yagi et al. (2014); Sham et
al. (2015). Any examples of the breaking of universal relations might provide
hints on the origin of the universality. It has been found that the I-Love-Q
relations become less accurate if the ellipticity of isodensity surfaces of a
neutron star displays a large variation inside the star Yagi et al. (2014),
such as the case for a hot protoneutron star Martinon et al. (2014). In Figure
5, we show the profiles of ellipticity $e$ Hartle and Thorne (1968),
normalized by the dimensionless spin parameter $a$, of constant energy-density
surfaces for four stellar models constructed from the NL3 EOS with
$\ln{\bar{\lambda}}_{\rm tid}=5$ and different values of $\Omega_{\rm
n}/\Omega_{\rm p}$. It is seen that, similar to the case of $\Omega_{\rm
n}/\Omega_{\rm p}=1$, the profiles for $\Omega_{\rm n}/\Omega_{\rm p}=-1$ and
10 are nearly constant over a large part of the star. As we have seen in
Figure 4, the Q-Love data for these cases are still relatively insensitive to
the EOS models, though the data for $\Omega_{\rm n}/\Omega_{\rm p}=-1$ deviate
largely from the single-fluid Q-Love relation. On the other hand, the case of
$\Omega_{\rm n}/\Omega_{\rm p}=0.1$ has a more significant variation in the
ellipticity profile, which correlates to the observation that the Q-Love data
in this case are more sensitive to the EOS models as shown in Figure 4. Our
results are thus consistent with the suggestion that the breaking of the
I-Love-Q relations is correlated with a large variation of the ellipticity
Yagi et al. (2014).

Figure 5:  Profiles of the ellipticity $e$ (normalized by the spin parameter
$a$) of constant energy-density surfaces for stellar models constructed from
the NL3 EOS with $\ln{\bar{\lambda}}_{\rm tid}=5$ and different values of
$\Omega_{\rm n}/\Omega_{\rm p}$.

##  5 Discussion

In summary, we have studied the I-Love-Q relations for superfluid neutron
stars based on a general relativistic two-fluid formulation. When neutrons
become superfluid, they decouple from the charged components to a first
approximation. We model the interior of a superfluid neutron star by a two-
fluid system containing a neutron superfluid and a ”proton” fluid containing
all charged particles. The neutron and proton fluids are also assumed to be
uniformly rotating with angular velocities $\Omega_{\rm n}$ and $\Omega_{\rm
p}$, respectively. The spin ratio $\Omega_{\rm n}/\Omega_{\rm p}$ is the most
important parameter in this work. We see from the results of our chosen EOS
models that the I-Love-Q relations originally discovered for single-fluid
ordinary neutron stars Yagi and Yunes (2013a, b) are still satisfied to high
accuracy for superfluid neutron stars when the two fluids are nearly co-
rotating so that $\Omega_{\rm n}/\Omega_{\rm p}\approx 1$. However, it is also
seen that the errors between our superfluid data and the universal relations
increase as the spin ratio deviates from unity in both directions (i.e.,
$\Omega_{\rm n}/\Omega_{\rm p}>1$ or $<1$). If $\Omega_{\rm n}/\Omega_{\rm p}$
could be different from unity by a few tens of percent, the deviation of the
Q-Love relation can be as large as $O(10\%)$ and this may have implications
for gravitational waveform models that make use of this relation (as discussed
in Section 1).

How much can the spin ratio $\Omega_{\rm n}/\Omega_{\rm p}$ actually deviate
from unity? We can try to address this question by making a connection to the
pulsar glitch phenomenon. The neutron superfluid inside a spinning isolated
neutron star is nearly decoupled from the proton fluid, which spins down
continuously due to electromagnetic radiation. The neutron superfluid is thus
expected to spin faster (i.e., $\Omega_{\rm n}/\Omega_{\rm p}>1$). When the
lag $\delta\Omega=\Omega_{\rm n}-\Omega_{\rm p}$ increases to some critical
value, the Magnus force induced on the superfluid vortices becomes strong
enough to lead to an unpinning of the vortices. The neutron superfluid will
then spin down and the proton fluid will spin up as a result of conservation
of angular momentum, leading to the glitch phenomenon. As a first
approximation, it seems reasonable to expect that the lag $\delta\Omega$ is
comparable to the change in the spin frequency during a glitch. The fractional
change of spin frequency $\Delta\nu/\nu$ observed in pulsar glitches ranges
from $10^{-11}$ to $10^{-5}$ Haskell and Melatos (2015). This would imply that
$\Omega_{\rm n}/\Omega_{\rm p}=1$ is satisfied to high accuracy for isolated
pulsars and the deviations of the universal relations due to the two-fluid
dynamics will not be observed astrophysically.

However, as one of our goals is to study whether the Q-Love relation remain
robust enough to be used in waveform modeling for binary neutron star systems,
it is instructive to also discuss the situation in binary systems. A new
physical process that may come into play in binary systems comparing to
isolated stars is the possibility of mass transfer between the two stars.
While it is not expected that mass transfer can occur when a binary neutron
star system enters the LIGO sensitivity band during the early inspiral phase,
it is possible that the first-born neutron star of the system might have spun
up due to accretion of matter by Roche-lobe overflow when its companion (i.e.,
the second-born neutron star) entered the post main-sequence evolution phase.
In this case, the proton fluid can spin up and rotate faster than the internal
neutron superfluid, leading to the situation $\Omega_{\rm n}/\Omega_{\rm
p}<1$. It is unclear whether this reverse condition, comparing to $\Omega_{\rm
n}/\Omega_{\rm p}>1$, could potentially increase the deviation of the spin
ratio away from unity. However, there may already be some hints from
observations.

Recently, three anti-glitches (i.e., $\Delta\nu<0$) in the accreting pulsar
NGC 300 ULX-1 have been observed by NICER Ray et al. (2019) and their
magnitudes $|\Delta\nu/\nu|\sim 10^{-4}$ are significantly larger than typical
glitches in isolated radio pulsars. Similarly, a glitch of magnitude
$|\Delta\nu/\nu|\sim 10^{-3}$ has also been seen in the accreting pulsar SXP
1062 Serim et al. (2017). The normal component is expected to spin faster than
the neutron superfluid due to accretion in these systems and hence
$\Omega_{\rm n}/\Omega_{\rm p}<1$. If the glitches in these accreting pulsars
are caused by the transfer of angular momentum between the superfluid and
normal components of the star, with the assumption that the lag $\delta\Omega$
between the two fluids is also characterized by the glitch magnitude, then
these recent observations would suggest that the lag magnitude
$|\delta\Omega|$ that can be sustained between the two fluids could be much
larger than that suggested by the glitches in isolated pulsars.

While it may be too optimistic (if not unrealistic) to expect that the spin
ratio $\Omega_{\rm n}/\Omega_{\rm p}$ for neutron stars in binary systems
could differ greatly from unity, to the level that the deviation of the Q-Love
relation can be as large as $O(10\%)$, a more detailed analysis of this issue
would require one to study not only the coupling between the normal and
superfluid components at the mesoscopic level inside a neutron star, but also
the formation and evolution of binary neutron star systems—issues that are
currently far from completely understood.

Methodology and formal analysis, C.H.Y.; supervision, L.M.L.; draft
preparation, C.H.Y., L.M.L.; review and editing, C.H.Y., L.M.L., N.A., G.C.
All authors have read and agreed to the published version of the
manuscript.This work is partially supported by a grant from the Research Grant
Council of the Hong Kong Special Administrative Region, China (Project No.
14300320). N.A. acknowledges funding by STFC through grant number
ST/R00045X/1. The authors declare no conflict of interest.

References

## References

  * LIGO Scientific Collaboration et al. (2017) LIGO Scientific Collaboration; Virgo Collaboration; Abbott, B.P.; Abbott, R.; Abbott, T.D.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; et al.  GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral.  Phys. Rev. Lett. 2017, 119, 161101. [CrossRef]
  * Annala et al. (2018) Annala, E.; Gorda, T.; Kurkela, A.; Vuorinen, A.  Gravitational-Wave Constraints on the Neutron-Star-Matter Equation of State.  Phys. Rev. Lett. 2018, 120, 172703. [CrossRef]
  * De et al. (2018) De, S.; Finstad, D.; Lattimer, J.M.; Brown, D.A.; Berger, E.; Biwer, C.M.  Tidal Deformabilities and Radii of Neutron Stars from the Observation of GW170817.  Phys. Rev. Lett. 2018, 121, 091102. [CrossRef]
  * Fattoyev et al. (2018) Fattoyev, F.J.; Piekarewicz, J.; Horowitz, C.J.  Neutron Skins and Neutron Stars in the Multimessenger Era.  Phys. Rev. Lett. 2018, 120, 172702. [CrossRef]
  * Most et al. (2018) Most, E.R.; Weih, L.R.; Rezzolla, L.; Schaffner-Bielich, J.  New Constraints on Radii and Tidal Deformabilities of Neutron Stars from GW170817.  Phys. Rev. Lett. 2018, 120, 261103. [CrossRef]
  * Tews et al. (2018) Tews, I.; Margueron, J.; Reddy, S.  Critical examination of constraints on the equation of state of dense matter obtained from GW170817.  Phys. Rev. C 2018, 98, 045804. [CrossRef]
  * Lim and Holt (2018) Lim, Y.; Holt, J.W.  Neutron Star Tidal Deformabilities Constrained by Nuclear Theory and Experiment.  Phys. Rev. Lett. 2018, 121, 062701. [CrossRef]
  * Malik et al. (2018) Malik, T.; Alam, N.; Fortin, M.; Providência, C.; Agrawal, B.K.; Jha, T.K.; Kumar, B.; Patra, S.K.  GW170817: Constraining the nuclear matter equation of state from the neutron star tidal deformability.  Phys. Rev. C 2018, 98, 035804. [CrossRef]
  * Li et al. (2019) Li, B.A.; Krastev, P.G.; Wen, D.H.; Zhang, N.B.  Towards understanding astrophysical effects of nuclear symmetry energy.  Eur. Phys. J. A 2019, 55, 117. [CrossRef]
  * Carson et al. (2019) Carson, Z.; Steiner, A.W.; Yagi, K.  Constraining nuclear matter parameters with GW170817.  Phys. Rev. D 2019, 99, 043010. [CrossRef]
  * Tsui and Leung (2005) Tsui, L.K.; Leung, P.T.  Probing the Interior of Neutron Stars with Gravitational Waves.  Phys. Rev. Lett. 2005, 95, 151101. [CrossRef]
  * Lau et al. (2010) Lau, H.K.; Leung, P.T.; Lin, L.M.  Inferring physical parameters of compact stars from their f-mode gravitational wave signals.  Astrophys. J. 2010, 714, 1234. [CrossRef]
  * Yagi and Yunes (2013a) Yagi, K.; Yunes, N.  I-Love-Q: Unexpected Universal Relations for Neutron Stars and Quark Stars.  Science 2013, 341, 365–368. [CrossRef]
  * Yagi and Yunes (2013b) Yagi, K.; Yunes, N.  I-Love-Q relations in neutron stars and their applications to astrophysics, gravitational waves, and fundamental physics.  Phys. Rev. D 2013, 88, 023009. [CrossRef]
  * Yagi (2014) Yagi, K.  Multipole Love relations.  Phys. Rev. D 2014, 89, 043011. [CrossRef]
  * Chan et al. (2014) Chan, T.K.; Sham, Y.H.; Leung, P.T.; Lin, L.M.  Multipolar universal relations between $f$-mode frequency and tidal deformability of compact stars.  Phys. Rev. D 2014, 90, 124023. [CrossRef]
  * Chakrabarti et al. (2014) Chakrabarti, S.; Delsate, T.; G‘̀urlebeck, N.; Steinhoff, J.  $I\text{$-$}Q$ Relation for Rapidly Rotating Neutron Stars.  Phys. Rev. Lett. 2014, 112, 201102. [CrossRef]
  * Pappas and Apostolatos (2014) Pappas, G.; Apostolatos, T.A.  Effectively Universal Behavior of Rotating Neutron Stars in General Relativity Makes Them Even Simpler than Their Newtonian Counterparts.  Phys. Rev. Lett. 2014, 112, 121101. [CrossRef]
  * Pappas (2015) Pappas, G.  Unified description of astrophysical properties of neutron stars independent of the equation of state.  Mon. Not. R. Astron. Soc. 2015, 454, 4066–4084. [CrossRef]
  * Breu and Rezzolla (2016) Breu, C.; Rezzolla, L.  Maximum mass, moment of inertia and compactness of relativistic stars.  Mon. Not. R. Astron. Soc. 2016, 459, 646–656. [CrossRef]
  * Bozzola et al. (2017) Bozzola, G.; Stergioulas, N.; Bauswein, A.  Universal relations for differentially rotating relativistic stars at the threshold to collapse.  Mon. Not. R. Astron. Soc. 2017, 474, 3557–3564. [CrossRef]
  * Luk and Lin (2018) Luk, S.S.; Lin, L.M.  Universal Relations for Innermost Stable Circular Orbits around Rapidly Rotating Neutron Stars.  Astrophys. J. 2018, 861, 141. [CrossRef]
  * Riahi et al. (2019) Riahi, R.; Kalantari, S.Z.; Rueda, J.A.  Universal relations for the Keplerian sequence of rotating neutron stars.  Phys. Rev. D 2019, 99, 043004. [CrossRef]
  * Sun et al. (2020) Sun, W.; Wen, D.; Wang, J.  New quasiuniversal relations for static and rapid rotating neutron stars.  Phys. Rev. D 2020, 102, 023039. [CrossRef]
  * Yagi and Yunes (2017) Yagi, K.; Yunes, N.  Approximate universal relations for neutron stars and quark stars.  Phys. Rep. 2017, 681, 1. [CrossRef]
  * Doneva and Pappas (2018) Doneva, D.D.; Pappas, G. Universal Relations and Alternative Gravity Theories.  In The Physics and Astrophysics of Neutron Stars; Rezzolla, L., Pizzochero, P., Jones, D.I., Rea, N., Vidaña, I., Eds.; Springer International Publishing: Cham, Switzerland, 2018; pp. 737–806. [CrossRef]
  * Lackey et al. (2019) Lackey, B.D.; Pürrer, M.; Taracchini, A.; Marsat, S.  Surrogate model for an aligned-spin effective-one-body waveform model of binary neutron star inspirals using Gaussian process regression.  Phys. Rev. D 2019, 100, 024002. [CrossRef]
  * Schmidt and Hinderer (2019) Schmidt, P.; Hinderer, T.  Frequency domain model of $f$-mode dynamic tides in gravitational waveforms from compact binary inspirals.  Phys. Rev. D 2019, 100, 021501. [CrossRef]
  * Andersson and Pnigouras (2019) Andersson, N.; Pnigouras, P.  The seismology of Love: An effective model for the neutron star tidal deformability.  arXiv 2019, arXiv:1905.00012.
  * Barkett et al. (2020) Barkett, K.; Chen, Y.; Scheel, M.A.; Varma, V.  Gravitational waveforms of binary neutron star inspirals using post-Newtonian tidal splicing.  Phys. Rev. D 2020, 102, 024031. [CrossRef]
  * Doneva et al. (2013) Doneva, D.D.; Yazadjiev, S.S.; Stergioulas, N.; Kokkotas, K.D.  Breakdown of I-Love-Q Universality in Rapidly Rotating Relativistic Stars.  Astrophys. J. 2013, 781, L6. [CrossRef]
  * Haskell et al. (2013) Haskell, B.; Ciolfi, R.; Pannarale, F.; Rezzolla, L.  On the universality of I-Love-Q relations in magnetized neutron stars.  Mon. Not. R. Astron. Soc. 2013, 438, L71–L75. [CrossRef]
  * Martinon et al. (2014) Martinon, G.; Maselli, A.; Gualtieri, L.; Ferrari, V.  Rotating protoneutron stars: Spin evolution, maximum mass, and I-Love-Q relations.  Phys. Rev. D 2014, 90, 064026. [CrossRef]
  * Marques et al. (2017) Marques, M.; Oertel, M.; Hempel, M.; Novak, J.  New temperature dependent hyperonic equation of state: Application to rotating neutron star models and $I\text{$-$}Q$ relations.  Phys. Rev. C 2017, 96, 045806. [CrossRef]
  * Paschalidis et al. (2018) Paschalidis, V.; Yagi, K.; Alvarez-Castillo, D.; Blaschke, D.B.; Sedrakian, A.  Implications from GW170817 and I-Love-Q relations for relativistic hybrid stars.  Phys. Rev. D 2018, 97, 084038. [CrossRef]
  * Lau et al. (2017) Lau, S.Y.; Leung, P.T.; Lin, L.M.  Tidal deformations of compact stars with crystalline quark matter.  Phys. Rev. D 2017, 95, 101302. [CrossRef]
  * Lau et al. (2019) Lau, S.Y.; Leung, P.T.; Lin, L.M.  Two-layer compact stars with crystalline quark matter: Screening effect on the tidal deformability.  Phys. Rev. D 2019, 99, 023018. [CrossRef]
  * Lombardo and Schulze (2001) Lombardo, U.; Schulze, H.J. Superfluidity in Neutron Star Matter.  In Physics of Neutron Star Interiors; Blaschke, D., Glendenning, N., Sedrakian, A., Eds.; Springer: Berlin, Germany, 2001; pp. 30–53. [CrossRef]
  * Haskell and Sedrakian (2018) Haskell, B.; Sedrakian, A. Superfluidity and Superconductivity in Neutron Stars.  In The Physics and Astrophysics of Neutron Stars; Rezzolla, L., Pizzochero, P., Jones, D.I., Rea, N., Vidaña, I., Eds.; Springer International Publishing: Cham, Switzerland, 2018; pp. 401–454. [CrossRef]
  * Sedrakian and Clark (2019) Sedrakian, A.; Clark, J.W.  Superfluidity in nuclear systems and neutron stars.  Eur. Phys. J. A 2019, 55, 1–56. [CrossRef]
  * Haskell and Melatos (2015) Haskell, B.; Melatos, A.  Models of pulsar glitches.  Int. J. Mod. Phys. D 2015, 24, 1530008. [CrossRef]
  * Anderson and Itoh (1975) Anderson, P.W.; Itoh, N.  Pulsar glitches and restlessness as a hard superfluidity phenomenon. Nature 1975, 256, 25–27. [CrossRef]
  * Baym et al. (1969) Baym, G.; Pethick, C.; Pines, D.  Superfluidity in Neutron Stars.  Nature 1969, 224, 673–674. [CrossRef]
  * Yu and Weinberg (2016) Yu, H.; Weinberg, N.N.  Resonant tidal excitation of superfluid neutron stars in coalescing binaries.  Mon. Not. R. Astron. Soc. 2016, 464, 2622–2637. [CrossRef]
  * Sekiguchi et al. (2011) Sekiguchi, Y.; Kiuchi, K.; Kyutoku, K.; Shibata, M.  Effects of Hyperons in Binary Neutron Star Mergers.  Phys. Rev. Lett. 2011, 107, 211101. [CrossRef]
  * Bernuzzi et al. (2016) Bernuzzi, S.; Radice, D.; Ott, C.D.; Roberts, L.F.; M‘̀osta, P.; Galeazzi, F.  How loud are neutron star mergers?  Phys. Rev. D 2016, 94, 024023. [CrossRef]
  * Perego et al. (2019) Perego, A.; Bernuzzi, S.; Radice, D.  Thermodynamics conditions of matter in neutron star mergers.  Eur. Phys. J. A 2019, 55. [CrossRef]
  * Andersson and Comer (2020) Andersson, N.; Comer, G.L.  Relativistic fluid dynamics: Physics for many different scales. _arXiv_ 2020, arXiv:gr-qc/2008.12069.
  * Gittins et al. (2020) Gittins, F.; Andersson, N.; Pereira, J.P.  Tidal deformations of neutron stars with elastic crusts.  Phys. Rev. D 2020, 101, 103025. [CrossRef]
  * Carter (1989) Carter, B.  Covariant theory of conductivity in ideal fluid or solid media.  In _Relativistic Fluid Dynamics_ ; Anile, A.M., Choquet-Bruhat, Y., Eds.; Springer: Berlin/Heidelberg, Germany, 1989; pp. 1–64. [CrossRef]
  * Comer and Langlois (1993) Comer, G.L.; Langlois, D.  Hamiltonian formulation for multi-constituent relativistic perfect fluids.  Class. Quantum Gravity 1993, 10, 2317–2327. [CrossRef]
  * Comer and Langlois (1994) Comer, G.L.; Langlois, D.  Hamiltonian formulation for relativistic superfluids.  Class. Quantum Gravity 1994, 11, 709–721. [CrossRef]
  * Carter and Langlois (1998) Carter, B.; Langlois, D.  Relativistic models for superconducting-superfluid mixtures.  Nucl. Phys. B 1998, 531, 478–504. [CrossRef]
  * Langlois et al. (1998) Langlois, D.; Sedrakian, D.M.; Carter, B.  Differential rotation of relativistic superfluid in neutron stars.  Mon. Not. R. Astron. Soc. 1998, 297, 1189–1201. [CrossRef]
  * Andersson (2021) Andersson, N.  A Superfluid Perspective on Neutron Star Dynamics.  Universe 2021, 7, 17. [CrossRef]
  * Comer et al. (1999) Comer, G.L.; Langlois, D.; Lin, L.M.  Quasinormal modes of general relativistic superfluid neutron stars.  Phys. Rev. D 1999, 60, 104025. [CrossRef]
  * Andersson and Comer (2001) Andersson, N.; Comer, G.L.  Slowly rotating general relativistic superfluid neutron stars.  Class. Quantum Gravity 2001, 18, 969–1002. [CrossRef]
  * Char and Datta (2018) Char, P.; Datta, S.  Relativistic tidal properties of superfluid neutron stars. Phys. Rev. D 2018, 98, 084010. 084010. [CrossRef]
  * Andersson et al. (2002) Andersson, N.; Comer, G.L.; Langlois, D.  Oscillations of general relativistic superfluid neutron stars.  Phys. Rev. D 2002, 66, 104002. [CrossRef]
  * Lin et al. (2008) Lin, L.M.; Andersson, N.; Comer, G.L.  Oscillations of general relativistic multifluid/multilayer compact stars.  Phys. Rev. D 2008, 78, 083008. [CrossRef]
  * Prix et al. (2005) Prix, R.; Novak, J.; Comer, G.L.  Relativistic numerical models for stationary superfluid neutron stars.  Phys. Rev. D 2005, 71, 043005. [CrossRef]
  * Sourie et al. (2016) Sourie, A.; Oertel, M.; Novak, J.  Numerical models for stationary superfluid neutron stars in general relativity with realistic equations of state.  Phys. Rev. D 2016, 93, 083004. [CrossRef]
  * Hartle (1967) Hartle, J.B.  Slowly Rotating Relativistic Stars. I. Equations of Structure.  Astrophys. J. 1967, 150, 1005. [CrossRef]
  * Hinderer (2008) Hinderer, T.  Tidal Love Numbers of Neutron Stars.  Astrophys. J. 2008, 677, 1216–1220. [CrossRef]
  * Damour and Nagar (2009) Damour, T.; Nagar, A.  Relativistic tidal properties of neutron stars.  Phys. Rev. D 2009, 80, 084035. [CrossRef]
  * Pani et al. (2015) Pani, P.; Gualtieri, L.; Ferrari, V.  Tidal Love numbers of a slowly spinning neutron star.  Phys. Rev. D 2015, 92, 124003. [CrossRef]
  * Landry (2017) Landry, P.  Tidal deformation of a slowly rotating material body: Interior metric and Love numbers.  Phys. Rev. D 2017, 95, 124058. [CrossRef]
  * Datta and Char (2020) Datta, S.; Char, P.  Effect of superfluid matter of a neutron star on the tidal deformability.  Phys. Rev. D 2020, 101, 064016. [CrossRef]
  * Comer and Joynt (2003) Comer, G.L.; Joynt, R.  Relativistic mean field model for entrainment in general relativistic superfluid neutron stars.  Phys. Rev. D 2003, 68, 023002. [CrossRef]
  * Comer (2004) Comer, G.L.  Slowly rotating general relativistic superfluid neutron stars with relativistic entrainment.  Phys. Rev. D 2004, 69, 123009. [CrossRef]
  * Kheto and Bandyopadhyay (2014) Kheto, A.; Bandyopadhyay, D.  Isospin dependence of entrainment in superfluid neutron stars in a relativistic model.  Phys. Rev. D 2014, 89, 023007. [CrossRef]
  * Fattoyev et al. (2010) Fattoyev, F.J.; Horowitz, C.J.; Piekarewicz, J.; Shen, G.  Relativistic effective interaction for nuclei, giant resonances, and neutron stars.  Phys. Rev. C 2010, 82, 055803. [CrossRef]
  * Glendenning and Moszkowski (1991) Glendenning, N.K.; Moszkowski, S.A.  Reconciliation of neutron-star masses and binding of the $\Lambda$ in hypernuclei.  Phys. Rev. Lett. 1991, 67, 2414–2417. [CrossRef]
  * Yagi et al. (2014) Yagi, K.; Stein, L.C.; Pappas, G.; Yunes, N.; Apostolatos, T.A.  Why I-Love-Q: Explaining why universality emerges in compact objects.  Phys. Rev. D 2014, 90, 063010. [CrossRef]
  * Sham et al. (2015) Sham, Y.H.; Chan, T.K.; Lin, L.M.; Leung, P.T.  Unveiling the university of I-Love-Q relations.  Astrophys. J. 2015, 798, 121. [CrossRef]
  * Hartle and Thorne (1968) Hartle, J.B.; Thorne, K.S.  Slowly Rotating Relativistic Stars. II. Models for Neutron Stars and Supermassive Stars.  Astrophys. J. 1968, 153, 807. [CrossRef]
  * Ray et al. (2019) Ray, P.S.; Guillot, S.; Ho, W.C.G.; Kerr, M.; Enoto, T.; Gendreau, K.C.; Arzoumanian, Z.; Altamirano, D.; Bogdanov, S.; Campion, R.; et al.  Anti-glitches in the Ultraluminous Accreting Pulsar NGC 300 ULX-1 Observed with NICER.  Astrophys. J. 2019, 879, 130. [CrossRef]
  * Serim et al. (2017) Serim, M.M.; Şahiner, Ş.; Çerri-Serim, D.; İnam, S.Ç.; Baykal, A.L.  Discovery of a glitch in the accretion-powered pulsar SXP 1062.  Mon. Not. R. Astron. Soc. 2017, 471, 4982–4989. [CrossRef]