The I-Love-Q Relations for Superfluid Neutron Stars

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.


Introduction
The groundbreaking detection of the first gravitational wave signal from a binary neutron star system GW170817 [1] 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., [2][3][4][5][6][7][8][9][10]). 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 [11][12][13][14][15][16][17][18][19][20][21][22][23][24]. These relations are "universal" in the sense that they are insensitive to the EOS models to the O(1%) level (see [25,26] 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 [12], 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 [27][28][29][30]. It is noted that the I-Love-Q arXiv:2105.00798v1 [astro-ph.HE] 23 Apr 2021 relations discovered by Yagi and Yunes [13,14] 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 [15] connects the l = 2 and l = 3 tidal deformabilities. The f-mode-Love relations found by Chan et al. [16] 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 [29], these universal relations [13][14][15][16] 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 [31], it was later found that the I-Q relation, in fact, remains approximately EOS-insensitive if one uses dimensionless parameters to characterize the rotation [17,18]. 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 [32]. Thermal effects relevant for protoneutron stars have also been shown to break the I-Love-Q relations [33,34]. It has also been demonstrated that the I-Love-Q relations are still satisfied to within ∼ 3% for hybrid star models with a strong first-order hadron-quark phase transition in the interior [35]. However, the I-Love relation can be broken if the quark matter core is in a crystalline color-superconducting state [36,37]. 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 (∼ 10 10 K). The internal temperature of a newborn neutron star is expected to drop quickly below the transition temperatures (∼ 10 9 K) for neutrons and protons to become superfluid and superconducting (see, e.g., [38][39][40]). 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 [41] 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 [42] 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 [43] 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 ∼ 10 7 K due to tidal heating in inspiralling binaries [44], significantly below the superfluid transition temperature (∼ 10 9 K). However, numerical simulations [45][46][47] 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 [48] 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 [49].
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 spininduced 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 =h = 1 unless otherwise noted.

General Relativistic Two-Fluid Formalism
Our study is based on the general relativistic two-fluid formalism developed by Carter and his collaborators (e.g., [50][51][52][53][54]). The formalism is built around the master function Λ(n 2 , p 2 , x 2 ), which is formed by the scalars constructed from the neutron n µ and proton p µ number density currents: n 2 = −n µ n µ , p 2 = −p µ p µ , and x 2 = −n µ p µ . 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 Λ, the stress-energy tensor is determined by where the generalized pressure Ψ is given by The chemical potential covectors µ α and χ α (the fluid four-momenta), respectively, for the neutrons and (conglomerate) protons are The coefficient A captures the entrainment effect [48,55] 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, and two Euler equations, The I-Love-Q relations connect the moment of inertia I and quadrupole moment Q of slowly rotating stars to the tidal deformability λ 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 [56][57][58] for more details.

Nonrotating Stars
The unperturbed background solution is assumed to be a spherically symmetric and static spacetime described by the metric The equilibrium structure of a nonrotating superfluid neutron star, assuming that the neutron and proton fluids coexist throughout the whole star, is studied in [56]. 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 [59,60]. 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 [60] 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 [56]: where primes denote derivatives with respect to r and the coefficients A 0 0 , B 0 0 , and C 0 0 are determined from the master function Λ (see [56] 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., Ψ(R) = 0) and the total mass M is obtained by

Slowly Rotating Stars
Rotating superfluid neutron stars have been studied using the general relativistic twofluid formalism [57,61,62]. In this work, we determine the moment of inertia and quadrupole moment of slowly rotating superfluid neutron stars by following the formalism developed in [57]. The spacetime for a rotating star is given by an axisymmetric and stationary metric of the form The neutron and proton fluids are assumed to be uniformly rotating with angular velocities Ω n and Ω p , respectively. In the slow rotation approximation, the frequencies Ω n , Ω p , Ω n Ω p are assumed to be much smaller than the characteristic Kepler frequency, which is proportional to √ M/R 3 in Newtonian theory, and one can then expand the metric functions up to second order in the rotation as where P 2 (cos θ) = (3 cos 2 θ − 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 framedragging term ω(r) is a first-order quantity determined by whereL n = ω − Ω n andL p = ω − Ω p . It should be noted that the equation reduces to the corresponding frame-dragging equation for a single-fluid model [63] when the two fluids are co-rotating (i.e., Ω n = Ω p ). Equation (12) is solved by choosing a central valueL n (0) so that the interior and exterior solutions match at the stellar surface to satisfy the conditioñ where J = J n + J p is the total angular momentum of the star. The neutron and proton angular momenta are given, respectively, by 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 The total moment of inertia is then given by I = I n + I p , which reduces to the definition for a single-fluid star if the two fluids are co-rotating (see also [62]). 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 (11) 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 θ). The equations and numerical schemes for determining the interior solutions of these quantities can be found in [57]. Outside the star, the problem is identical to the single-fluid case and the metric functions can be obtained analytically [63]. In particular, the exterior solution for h 2 is given by 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 θ) becomes the perturbed Newtonian potential and the quadrupole moment can be read off from the coefficient of the P 2 (cos θ)/r 3 term: The normalized (dimensionless) moment of inertiaĪ and quadrupole momentQ that appear in the I-Love-Q relations are defined bȳ where a = J/M 2 is the dimensionless spin parameter.

Tidally Deformed Nonrotating Stars
The tidal deformability λ 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 λ tid for nonrotating single-fluid neutron stars is well established [64,65]. The tidal deformations of slowly rotating single-fluid neutron stars have also been studied in [66,67]. The formulation has recently been extended to nonrotating superfluid neutron stars within the general relativistic two-fluid formalism [58,68]. 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 The perturbed Einstein equations impose the condition H 0 = −H 2 ≡ H. The linearized metric and fluid equations inside the star yield the following equation for determining the tidal deformability [58] where the function m is defined by e λ(r) = (1 − 2m(r)/r) −1 and Noticing that Ψ is the pressure and Λ is the negative energy density of the unperturbed nonrotating background star, one can compare Equation (20) with the corresponding equation in a single-fluid situation (see Equation (15) of [64]) for a barotropic EOS P(ρ) and observe that the term (ρ + P)/(dP/dρ) in the single-fluid equation is now replaced by the function −g in the two-fluid equation. It can also be shown that Equation (20) reduces to Equation (15) of [64] when the master function depends only on one particle density (e.g., Λ(p 2 )) in the singlefluid limit. Once the interior problem (Equation (20)) 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 E ij results in the following expansion of the metric function g tt in the star's local asymptotic rest frame [64]: where Q ij is the tidally induced traceless quadrupole moment tensor of the star. The tidal deformability λ tid of the star is defined by Q ij = −λ tid E ij . It is also convenient to define the dimensionless (l = 2) tidal Love number k 2 ≡ 3λ tid /2R 5 . Outside the star, the solutions of Equation (20) are given by the associated Legendre functions Q 2 2 (x) and P 2 2 (x): where x = r/M − 1. Using Equation (22), 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 (20) at the surface. The tidal Love number can then be expressed as [64] where C = M/R is the compactness and y = RH (R)/H(R). Finally, the normalized (dimensionless) tidal deformabilityλ tid that appears in the I-Love-Q relations is defined bȳ

Two-Fluid Polytropic Model
As discussed in the previous section, the master function Λ(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 [56] where the neutron and proton masses are assumed to be equal (m n = m p ) for simplicity. We choose the parameters σ n = 0.2m n , σ p = 2m n , β n = 2.3, and β p = 1.95 as in [57] 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 A, B, etc., (see Equation (8)) can be easily computed for this model [56]. In particular, we have A = A 0 0 = 0.

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 [69,70] have calculated the entrainment effect in a relativistic σ − ω model. Kheto and Bandyopadhyay [71] later extended the study to include the isospin dependence by using a relativistic σ − ω − ρ model. More recently, Char and Datta [58,68] used the same RMF model to study the tidal deformability of superfluid neutron stars. In this work, we also use the relativistic σ − ω − ρ model to construct a master function. The Lagrangian density for the system is given by where Ψ B is the Dirac spinor for baryons B with baryon mass m B and γ µ denotes the Dirac matrices. The isospin operator is denoted by τ B . In this theory, we have the self-interacting scalar σ field, the vector omega field ω µ , and the isovector rho field ρ µ . The latter two fields have the corresponding field tensors ω µν and ρ µν . 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 µ p = (k 0 , 0, 0, K) [69]. In the limit K → 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 [58] where k n = (3π 2 n) 1/3 and k p = (3π 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 e ). The background value for the Dirac effective mass m * | 0 is determined by the following transcendental equation We refer the reader to [58] for the expressions of the various coefficients A, B, etc. (see Equation (8)). In contrast to the two-fluid polytropic model, it should be noted that the coefficient A, which is responsible for the entrainment effect, is nonzero in this RMF model. The coupling constants c 2 σ ≡ (g σ /m σ ) 2 , c 2 ω ≡ (g ω /m ω ) 2 , c 2 ρ ≡ (g ρ /m ρ ) 2 , b, and c in the master function are determined by the nuclear matter saturation properties [71]. In the following, we shall use the same NL3 [72] and GM1 [73] parameter sets (see Table 1) as employed in [58].
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 neutron stars constructed from our EOS models.

Numerical Results
The I-Love-Q universal relations [13,14] originally found for single-fluid neutron stars connect the three dimensionless quantitiesĪ,Q, andλ tid . While one might naively expect that I andQ are somehow related to each other as they both characterize the effects of rotation, the fact that these quantities together withλ tid are connected by the following approximately EOS-independent relations are quite surprising [14]: where (x i , y i ) are any two ofĪ,Q, andλ 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 [14]). 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Ī,Q, andλ 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 [56]. 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 Ω n and Ω p . However, instead of providing Ω n and Ω p as input parameters, the system of equations can be scaled in such a way that it is sufficient to specify the spin ratio Ω n /Ω p in the calculation [57]. We will thus present our results for stellar sequences specified by Ω n /Ω 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 inertiā I and quadrupole momentQ can be determined for the slowly rotating star as outlined in Section 2.2. The computation of the dimensionless tidal deformabilityλ 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Ī against lnλ tid (i.e., the I-Love relation) for sequences of stars with spin ratios Ω n /Ω p = 0.4, 0.7, 1, 1.3, and 1.6. The solid line is the fitting curve (Equation (30)) for the I-Love relation for single-fluid ordinary neutron stars [14]. The lower panel shows the relative error, E = (ŷ − y)/y, between the numerical data,ŷ, and the fitting curve, y. We see from the upper panel that the numerical data can still match the fitting curve for singlefluid stars very well. However, the effects of two-fluid dynamics become more apparent in the lower panel. The co-rotating case Ω n /Ω 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., Ω n /Ω p > 1), the error |E| increases with Ω n /Ω p . On the other hand, |E| increases as Ω n /Ω p decreases from 0.7 to 0.4 for the opposite situation where the proton fluid rotates faster. Sinceλ tid decreases with increasing compactness, the numerical results thus illustrate that the error |E| increases with the compactness for a given sequence.
In Figure 2b,c, we plot lnQ against lnλ tid and lnĪ against lnQ, 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 [14]. 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 Ω n /Ω p , we see that the I-Love-Q relations are broken only for the case Ω n /Ω 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 Ω n /Ω 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| ∼ 0.1 level for the case Ω n /Ω 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 ofQ by expressing it in terms ofλ 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 Ω n /Ω 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. 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 A = 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 1). In each figure, we consider three sequences of Ω n /Ω 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 Ω n /Ω 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| ∼ 0.1 level.
Let us now consider more extreme situations to explore the breaking of the Q-Love relation when Ω n /Ω p deviates significantly from unity. We have studied the trend of the numerical results in the lnQ − lnλ 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 Ω n /Ω 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 Ω n /Ω p = 10 (red data points). Increasing the spin ratio further (e.g., Ω n /Ω 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 Ω n /Ω 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 Ω n /Ω 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 Ω n /Ω 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 Ω n /Ω 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 Ω n /Ω 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 Ω n /Ω p > 0. It can be seen from the results of Ω n /Ω p = −10 (green data points) that the data are insensitive to the EOS models and agree quite well with the results of Ω n /Ω p = 10. In fact, this should be expected as the system is dominated by the neutron superfluid if |Ω n /Ω p | 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 2, we compare the values of lnQ at lnλ 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 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 star (see Figure 1). In the table, the numerical values inside the parentheses are the percentage differences between the data and the value of lnQ 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. The results presented in Figure 4 may have little astrophysical relevance as |Ω n /Ω 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 Ω n /Ω 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 [74,75]. 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 [74], such as the case for a hot protoneutron star [33]. In Figure 5, we show the profiles of ellipticity e [76], normalized by the dimensionless spin parameter a, of constant energy-density surfaces for four stellar models constructed from the NL3 EOS with lnλ tid = 5 and different values of Ω n /Ω p . It is seen that, similar to the case of Ω n /Ω p = 1, the profiles for Ω n /Ω 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 Ω n /Ω p = −1 deviate largely from the single-fluid Q-Love relation. On the other hand, the case of Ω n /Ω 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 [74].

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 Ω n and Ω p , respectively. The spin ratio Ω n /Ω 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 [13,14] are still satisfied to high accuracy for superfluid neutron stars when the two fluids are nearly co-rotating so that Ω n /Ω p ≈ 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., Ω n /Ω p > 1 or < 1). If Ω n /Ω 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 Ω n /Ω 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., Ω n /Ω p > 1). When the lag δΩ = Ω n − Ω 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 δΩ is comparable to the change in the spin frequency during a glitch. The fractional change of spin frequency ∆ν/ν observed in pulsar glitches ranges from 10 −11 to 10 −5 [41]. This would imply that Ω n /Ω 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 Ω n /Ω p < 1. It is unclear whether this reverse condition, comparing to Ω n /Ω 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., ∆ν < 0) in the accreting pulsar NGC 300 ULX-1 have been observed by NICER [77] and their magnitudes |∆ν/ν| ∼ 10 −4 are significantly larger than typical glitches in isolated radio pulsars. Similarly, a glitch of magnitude |∆ν/ν| ∼ 10 −3 has also been seen in the accreting pulsar SXP 1062 [78]. The normal component is expected to spin faster than the neutron superfluid due to accretion in these systems and hence Ω n /Ω 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 δΩ between the two fluids is also characterized by the glitch magnitude, then these recent observations would suggest that the lag magnitude |δΩ| 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 Ω n /Ω 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.