Effect of superfluid matter of neutron star on the tidal deformability

TL;DR

This paper investigates how superfluidity in the neutron star core influences its tidal deformability, highlighting the importance of multiple tidal parameters for gravitational wave analysis.

Contribution

It introduces a dual-layer model with superfluid core and ordinary matter envelope to analyze tidal responses in neutron stars.

Findings

Superfluidity affects tidal Love numbers significantly.

Entrainment between fluids alters the tidal response.

Multiple tidal parameters are needed to probe superfluidity.

Abstract

We study the effect of superfluidity on the tidal response of a neutron star in a general relativistic framework. In this work, we take a dual-layer approach where the superfluid matter is confined in the core of the star. Then, the superfluid core is encapsulated with an envelope of ordinary matter fluid which acts effectively as the low-density crustal region of the star. In the core, the matter content is described by a two-fluid model where only the neutrons are taken as superfluid and the other fluid consists of protons and electrons making it charge neutral. We calculate the values of various tidal love numbers of a neutron star and discuss how they are affected due to the presence of entrainment between the two fluids in the core. We also emphasize that more than one tidal parameter is necessary to probe superfluidity with the gravitational wave from the binary inspiral.

Tables

Table 1. Table 1: Nucleon-meson coupling constants in the NL3 and GM1 sets are taken from Refs. Glendenning1991 ; Fattoyev2010 . The coupling constants are obtained by reproducing the saturation properties of symmetric nuclear matter as detailed in the text. All the parameters are in fm 2 , except b 𝑏 b and c 𝑐 c which are dimensionless.

	$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

Table 2. Table 2: Values of tidal deformabilities for 1.4 ​ M ⊙ 1.4 subscript 𝑀 direct-product 1.4M_{\odot}

	${\Lambda }_2^{el}$ 1-fl	${\Lambda }_2^{el}$ 2-fl	${\Lambda }_3^{el}$ 1-fl	${\Lambda }_3^{el}$ 2-fl	${\Lambda }_2^{mag}$ 1-fl	${\Lambda }_2^{mag}$ 2-fl
NL3	1268	1391	3455	4015.5	-7.9	-8.4
GM1	903	979	2241	2440.5	-6.2	-6.6

Table 3. Table 3: Estimated numerical coefficients for the fitting formula of the multipole Love relations

$y$	$x$	fluid type	$a$	$b$	$c$	$d$	$f$
$C$	${\Lambda }_2^{el}$	1-fl	0.364	-0.037	0.001	-	-
$C$	${\Lambda }_2^{el}$	2-fl	0.349	-0.031	0.0	-	-
$\log {\Lambda }_3^{el}$	${\Lambda }_2^{el}$	1-fl	-1.19050	1.14172	0.04698	-0.00447	0.00017
$\log {\Lambda }_3^{el}$	${\Lambda }_2^{el}$	2-fl	-1.09537	1.11772	0.04237	-0.00290	0.00007
$\log -{\Lambda }_2^{mag}$	${\Lambda }_2^{el}$	1-fl	-1.99175	0.44237	0.02082	-0.00039	-0.00000
$\log -{\Lambda }_2^{mag}$	${\Lambda }_2^{el}$	2-fl	-1.87476	0.35685	0.04403	-0.00317	0.00012

Equations

$$T_{\nu }^{\mu }=\Psi {\delta }_{\nu }^{\mu }+p^{\mu }{\chi }_{\nu }+n^{\mu }{\mu }_{\nu },$$

$$\Psi =\Lambda -n^{\rho }{\mu }_{\rho }-p^{\rho }{\chi }_{\rho },$$

$${\mu }_{\mu }=\mathcal{B}{n}_{\mu }+\mathcal{A}{p}_{\mu },{\chi }_{\mu }=\mathcal{C}{p}_{\mu }+\mathcal{A}{n}_{\mu },$$

$$\mathcal{A}=-\frac{\partial \Lambda }{\partial x^2},\mathcal{B}=-2\frac{\partial \Lambda }{\partial n^2},\mathcal{C}=-2\frac{\partial \Lambda }{\partial p^2}.$$

$${\nabla }_{\mu }{n}^{\mu }={\nabla }_{\mu }{p}^{\mu }=0.$$

$$n^{\mu }\nabla {}_{[\mu }{\mu }_{\nu ]}=p^{\mu }\nabla {}_{[\mu }{\chi }_{\nu ]}=0,$$

$${\mathcal{L}}_B$$

$$d{s}_0^2=g_{\alpha \beta }^{(0)}d{x}^{\alpha }d{x}^{\beta }=-e^{\nu (r)}d{t}^2+e^{\kappa (r)}d{r}^2+r^2(d{\theta }^2+{\sin }^2\theta d{\varphi }^2)$$

$${\kappa }^{\prime }$$

$${\nu }^{\prime }$$

$${\mathcal{A}}_0^0{|}_0{p}^{\prime }+{\mathcal{B}}_0^0{|}_0{n}^{\prime }+\frac{1}{2}\mu {|}_0{\nu }^{\prime }$$

$${\mathcal{C}}_0^0{|}_0{p}^{\prime }+{\mathcal{A}}_0^0{|}_0{n}^{\prime }+\frac{1}{2}\chi {|}_0{\nu }^{\prime }$$

$$\begin{array}{rl}{\mathcal{A}}_0^0& =\mathcal{A}+2\frac{\partial \mathcal{B}}{\partial p^2}np+2\frac{\partial \mathcal{A}}{\partial n^2}{n}^2+2\frac{\partial \mathcal{A}}{\partial p^2}{p}^2+\frac{\partial \mathcal{A}}{\partial x^2}np,\\ {\mathcal{B}}_0^0& =\mathcal{B}+2\frac{\partial \mathcal{B}}{\partial n^2}{n}^2+4\frac{\partial \mathcal{A}}{\partial n^2}np+\frac{\partial \mathcal{A}}{\partial x^2}{p}^2,\\ {\mathcal{C}}_0^0& =\mathcal{C}+2\frac{\partial \mathcal{C}}{\partial p^2}{p}^2+4\frac{\partial \mathcal{A}}{\partial p^2}np+\frac{\partial \mathcal{A}}{\partial x^2}{n}^2.\end{array}$$

$$m_{*}^{\prime}|_{0}=\frac{\partial m_{*}}{\partial k_{n}}\bigg{|}_{0}k_{n}^{\prime}+\frac{\partial m_{*}}{\partial k_{p}}\bigg{|}_{0}k_{p}^{\prime},$$

$$T_{\nu }^{\mu }=p{\delta }_{\nu }^{\mu }+(\rho +p)u^{\mu }{u}_{\nu },$$

$${\kappa }^{\prime }$$

$${\nu }^{\prime }$$

$$M=-4\pi {\int }_0^{R_c}dr{r}^2\Lambda {|}_0(r)+4\pi {\int }_{R_c}^{R}dr{r}^2\rho (r),$$

$$g_{\alpha \beta }=g_{\alpha \beta }^{(0)}+\delta g_{\alpha \beta },$$

$$\begin{array}{rl}& \sum _{l}diag[-e^{\nu (r)}{H}_0^{(l)}(r),e^{\kappa (r)}{H}_2^{(l)}(r),r^2{K}^{(l)}(r),r^2{\sin }^2\theta K^{(l)}(r)]\\ & \times P_l(\theta )\end{array}$$

$$\begin{array}{rl}\delta & G_{\theta }^{\theta }-\delta G_{\varphi }^{\varphi }=0\\ ⟹& H_0^{(l)}=-H_2^{(l)}\equiv H^{(l)}\end{array}$$

$$\begin{array}{rl}\delta & G_{\theta }^{\theta }+\delta G_{\varphi }^{\varphi }=16\pi \delta \Psi \\ ⟹& 2\delta \Psi =P_l(\theta )H^{(l)}(\Lambda -\Psi )\end{array}$$

$$\begin{array}{rl}\delta & G_{\theta }^r=0\\ ⟹& K^{(l)\prime }+H^{(l)\prime }+H^{(l)}{\nu }^{\prime }=0,\end{array}$$

$$\begin{array}{rl}{K}^{(l)}=& \frac{-r^2{\nu }^{\prime }{H}^{(l)\prime }}{(l^2+l-2)e^{\kappa }}\\ & +\frac{H^{(l)}\{2-r^2{\nu }^{\prime 2}+e^{\kappa }(8\pi r^2(\Psi -\Lambda )-l(l+1))\}}{(l^2+l-2)e^{\kappa }}\end{array}$$

$${\partial }_t\delta {\mu }_i={\partial }_i\delta {\mu }_t,{\partial }_t\delta {\chi }_i={\partial }_i\delta {\chi }_t.$$

$$\begin{array}{rl}& \delta {\mu }_0=(A_0^0\delta p+B_0^0\delta n)u^0\delta g_{00}+u^0\frac{\mu }{2}\delta g_{00}\\ & \delta {\chi }_0=(A_0^0\delta n+C_0^0\delta p)u^0\delta g_{00}+u^0\frac{\chi }{2}\delta g_{00}.\end{array}$$

$$\begin{array}{rl}& \delta n=\frac{(\chi A_0^0-\mu C_0^0)}{(B_0^0{C}_0^0-A_0^{02})}\frac{H^{(l)}{P}_l(\theta )}{2}\\ & \delta p=\frac{(\mu A_0^0-\chi B_0^0)}{(B_0^0{C}_0^0-A_0^{02})}\frac{H^{(l)}{P}_l(\theta )}{2}\end{array}$$

$$\begin{array}{rl}\delta \Lambda & =\frac{\partial \Lambda }{\partial x^2}\delta x^2+\frac{\partial \Lambda }{\partial p^2}\delta p^2+\frac{\partial \Lambda }{\partial n^2}\delta n^2\\ & =-[(An+Cp)\delta p+(Ap+Bn)\delta n]\\ & =-g\frac{H^{(l)}}{2}{P}_l(\theta ).\end{array}$$

$$\begin{array}{r}g=\frac{{\mu }^2{C}_0^0+{\chi }^2{B}_0^0-2\mu \chi A_0^0}{A_0^{02}-B_0^0{C}_0^0}\end{array}$$

$$\delta G_t^t-\delta G_r^r=-4\pi g{H}^{(l)}{P}_l(\theta )+4P_l(\theta )\pi H^{(l)}(\Psi -\Lambda )$$

$$\begin{split}&H^{(l)\prime\prime}+H^{(l)\prime}\bigg{[}4\pi re^{\kappa}(\Lambda+\Psi)+\frac{e^{\kappa}+1}{r}\bigg{]}\\ &+H^{(l)}\bigg{[}4\pi e^{\kappa}(-5\Lambda+9\Psi-g)-\nu^{\prime 2}-\frac{l(l+1)e^{\kappa}}{r^{2}}\bigg{]}=0\end{split}$$

$$\begin{split}&H^{(l)\prime\prime}+H^{(l)\prime}\bigg{[}4\pi re^{\kappa}(-\rho+p)+\frac{e^{\kappa}+1}{r}\bigg{]}\\ &+H^{(l)}\bigg{[}4\pi e^{\kappa}(5\rho+9p+\frac{\rho+p}{dp/d\rho})-\nu^{\prime 2}-\frac{l(l+1)e^{\kappa}}{r^{2}}\bigg{]}=0\end{split}$$

Full text

I Introduction

The observation of a gravitational wave (GW) from the binary neutron star (BNS) merger event GW170817 has allowed us to study the physics of the extreme environment of highly dense matter at strong gravity Abbott2017 ; Abbott2018 . During the orbital evolution, the tidal interaction between the stars of the binary deforms both of them. These deformations can be measured in terms of the relativistic tidal Love numbers of the stars Flanagan2008 ; Hinderer2008 ; Damour2009 ; Binnington2009 ; Hinderer et al. (2010). Precise measurements of these parameters from the GW signal during the inspiral phase can be extremely useful to study the nature and the equation of state (EOS) of the supranuclear matter inside a neutron star (NS) Agathos2015 ; Takami2014 ; Bose2018 . This is why a huge effort has been made to understand the modification of waveforms due to the tidal Love numbers and their measurability and distinguishability of different EOSs Vines2011 ; Damour2012 ; Read2013 ; DelPozzo2013 ; Wade2014 ; Favata2014 ; Hotokezaka2016 . Moreover, one can also infer on the fluid nature of those objects. As these stars are supposedly very old, their core temperature should be below the critical transition temperature for the BCS-like pair formation Sedrakian2018 . Therefore, one can expect superfluid (SF) neutrons and superconducting protons to form at the core of the star and superfluid neutrons in the inner crust Migdal1959 ; Clark1992 . Pulsar glitches and the rapid cooling of the NS in Cassiopeia A are examples which are explicable invoking superfluid matter inside NS Baym1975 ; Anderson1975 ; Page2011 ; Shternin2011 . These changes in the fluid nature of the star from a single-fluid to a multi-fluid object can influence its deformability in a non-trivial way Char2018 . Recently, we have investigated the role of superfluidity for the $\ell=2$ electric-type tidal Love number $k_{2}$ and the corresponding tidal deformability $\Lambda_{k_{2}}$ Char2018 , (hereafter, paper I). In this work we have modeled the star as a non-rotating sphere of superfluid nuclear matter. We had adopted the two-fluid model where one fluid is the neutron superfluid and the other is the normal charge-neutral fluid comprising protons and electrons Carter1989 ; Comer1994 ; Carter1995 ; carter1998_1 ; carter1998_2 ; Langois1998 ; Prix2000 . We found that the inclusion of superfluidity manifests significant change in $\Lambda_{k_{2}}$ compared to the non-superfluid case.

However, a neutron star is also a multi-layered object i.e. the phases of matter differ significantly from the crust to the core. As has been known that the property of low density nuclear matter is correlated directly with the radius, one has to take into account a proper crust model in the calculation. To do so, we follow the method described in Ref.Andersson_2002 , where the properties of the superfluid region inside the core is appropriately matched to the normal fluid envelope encapsulating the core. Therefore, the superfluid neutrons are confined in the core where as the envelope acts as the low density region of the star. Although, we do not consider the elasticity of the crustal region in our formalism, this dual-layer core-envelop approach can approximate the structure of the star with a crust. Since crustal elasticity does not bring considerable change in the Love numbers it is unnecessary to include it hereBiswas2019elasticity . We also study the junction conditions for the perturbed quantities of interest in detail.

At this point, it is important to note that when we speak of the deviation of $\Lambda_{k_{2}}$ due to the superfluid nature, we bring an ambiguity in our interpretation of the observed $\Lambda_{k_{2}}$. The value of $\Lambda_{k_{2}}$ in two-fluid calculation for a particular EOS model can be similar to the value in a single-fluid calculation for another EOS. So, we cannot distinguish between the EOS and also probe the fluid nature of matter at the same time with the measurement of $\Lambda_{k_{2}}$. One possible way to break the degeneracy is to have measurements of other Love numbers which have much smaller effects on the waveform. This gives us a primary motivation to study higher order electric-type Love numbers and magnetic-type Love numbers in the case of a superfluid star.

The paper is organized as follows. In Sec. II, we first discuss the two-fluid formalism followed by the calculation of the equilibrium structure along with a brief overview of the RMF model of dense matter to calculate the assorted matter coefficients of the model. Next, in Secs. III and IV, we derive the framework for even and odd parity tidal perturbations in the two-fluid model respectively. In Sec. V we discuss how the tidal Love numbers are calculated. Then, in Sec. VI we discuss our results. We assume $c=G=1$ and use the metric signature $(-,+,+,+)$ throughout the article.

II General relativistic superfluid neutron star

The main ingredients of the superfluid formalism have been developed and discussed in several works Carter1989 ; Comer1994 ; Carter1995 ; carter1998_1 ; carter1998_2 ; Langois1998 ; Comer1999 ; Prix2000 ; Andersson2001 . To incorporate SF matter inside NSs we follow a two-fluid model with entrainment. The central quantity of this formalism is the master function, $\Lambda$. It depends on three scalars, $n^{2}=-n^{\mu}n_{\mu}$, $p^{2}=-p^{\mu}p_{\mu}$, and $x^{2}=-n^{\mu}p_{\mu}$, where $n^{\mu}$ and $p^{\mu}$ are the number density currents of the neutron and proton, respectively. When the fluids are co-moving, $-\Lambda(n^{2},p^{2},x^{2})$ represents the total thermodynamic energy density. The energy-momentum tensor takes the following form,

[TABLE]

where, $\Psi$ is the generalized pressure, and it can be expressed as,

[TABLE]

where, $\chi_{\nu}$ and $\mu_{\nu}$ are, respectively, the chemical potential co-vectors of the proton and the neutron fluids.

[TABLE]

where the ${\cal A,B}$ and ${\cal C}$ coefficients are defined as follows,

[TABLE]

The expressions for $\mu_{\mu}$ and $\chi_{\mu}$ in Eq. 3 make the entrainment effect vivid. Momentum of the one component carries along some of the mass current of the other component when ${\cal A}\neq 0$. Thus, if ${\cal A}=0$ the master function becomes “entrainment-free” , implying that it is independent of $x^{2}$. The conservation equation for $n^{\mu}$ and $p^{\mu}$ implies,

[TABLE]

They also satisfy a set of Euler type equations Comer1999 ,

[TABLE]

where, the square brackets represent the antisymmetrization of the closed indices.

II.1 Equation of state of nuclear matter

We have calculated the master function $(\Lambda)$ using the $\sigma$-$\omega$-$\rho$ model with self-interaction in the RMF approximation Comer2003 ; Comer2004 ; Kheto2014 ; Kheto2015 . The Lagrangian of the theory is as follows,

[TABLE]

where, $m_{B}$ is the baryon mass. We use the nucleon mass $m$ as the average of the baryon masses. The Dirac effective mass $m_{*}$ has been defined as $m_{*}=m-g_{\sigma}\sigma$. The $\sigma$, $\omega$ and $\rho$ mesons represent the scalar, vector and vector-isovector interactions, respectively. {\mbox{\boldmath\tau}}_{B} is the isospin operator. $\Omega_{\mu\nu}$ and {\mbox{\boldmath\mathrm{P}}}_{\mu\nu} are the field tensors for $\omega$ and $\rho$ mesons respectively. For the two-fluid system, we choose a frame in such a way that the neutrons have zero spatial momentum and the proton momentum has a boost along the z-direction as $k_{p}^{\mu}=\left(k_{0},0,0,K\right)$. We follow the procedure as described in Refs. Kheto2014 ; Kheto2015 to solve the meson field equations and numerically evaluate the master function $\Lambda$, generalized pressure $\Psi$ etc. in the limit $K\rightarrow 0$.

We consider a normal fluid envelope around the superfluid core of the star to account for the behavior of the low density region of a NS. We assume this region to be free of superfluid neutrons. This assumption does not affect the macroscopic structure of the star. To describe the matter in this region, we employ the EOS for the inner crust calculated by Grill et al. grill2014 . We smoothly join the EOS by keeping the pressure continuous from the two-fluid region to the envelope. We also use the DH EOS haensel2007 for the outer part of the envelope.

II.2 Equilibrium configuration

We take the background metric of the star to be static and spherically symmetric. Under such assumptions, the metric can be written in the Schwarzschild form as follows,

[TABLE]

This metric structure is valid both in the core and the envelope. Only the energy-momentum tensor changes from one region to another.

II.2.1 Superfluid core

In the core the energy momentum tensor will take that of an SF matter, as has been described in Eq.(13). The two metric functions can then be evaluated from the Einstein’s equations as follows,

[TABLE]

By the following equations the radial profiles for $n(r)$ and $p(r)$ are determined,Comer1999 ,

[TABLE]

where,

[TABLE]

The two Fermi wave numbers $k_{n}$ and $k_{p}$ are the variables that are more appropriate for the RMF calculations. Thus, we substitute the number densities with the Fermi wave numbers using $n=\frac{k_{n}^{3}}{3\pi^{2}}$ and $p=\frac{k_{p}^{3}}{3\pi^{2}}$, and solve for $k_{n}$ and $k_{p}$ instead. We determine the Dirac effective mass $m_{*}|_{0}(k_{n},k_{p})$ using the method discussed in Comer2003 . The transcendental algebraic relation in Eq. 71 is turned into a differential equation using,

[TABLE]

where $k_{n}^{\prime}$ and $k_{p}^{\prime}$ are calculated from Eq. 10. The prime in the equation represents a radial derivative and a zero subscript represents that $K\rightarrow 0$ has been taken after the partial derivatives are calculated. We put the boundary condition at the center and the surface of the star. A non-singularity condition at the center imposes $\kappa(0)=0$ and $\kappa^{\prime}(0)$ and $\nu^{\prime}(0)$ vanishes. Together with Eq. 10 this condition imposes $k_{n}^{\prime}(0)=k_{p}^{\prime}(0)=0$. Necessary expressions for all the matter quantities used in our calculations \big{(}\Lambda|_{0},\Psi|_{0},\mu|_{0},\chi|_{0},m_{*}|_{0},{\cal A}|_{0},{\cal B}|_{0},{\cal C}|_{0},{\cal A}^{0}_{0}|_{0},{\cal B}^{0}_{0}|_{0},{\cal C}^{0}_{0}|_{0}, \left.\frac{\partial m_{*}}{\partial k_{n}}\right|_{0},\left.\frac{\partial m_{*}}{\partial k_{p}}\right|_{0}\big{)} can be found in Appendix B.

II.2.2 Normal fluid envelope

In the envelope the matter is modeled as one component normal fluid (NF). Therefore the energy momentum tensor can be written as,

[TABLE]

where $\rho$ and $p$ are the energy density and the pressure of the fluid in the envelope, respectively. And $u^{\mu}$ is the four velocity of the fluid.

Using this form of energy-momentum tensor equation for the two metric functions can be evaluated from the Einstein’s equations as follows,

[TABLE]

The continuity of the metric variables at the junction of the SF core and the normal fluid envelope has been discussed in appendix A.1. The surface of the star $r=R$ implies that the total mass of the star is,

[TABLE]

and $\Psi|_{0}(R_{c})=p(R_{c})$ and $p(R)=0$, where $R_{c}$ is the junction between the SF core and the NF envelope.

III Even Parity Perturbation Equations for zero frequency mode

To calculate the electric type tidal Love no., perturbation of the static and spherically symmetric background needs to be calculated. For this purpose we decompose the metric as Thorne1967 ,

[TABLE]

where, $g^{(0)}_{\alpha\beta}$ and $\delta g_{\alpha\beta}$ are the background and the perturbed part of the metric respectively.

We decompose the metric and the fluid perturbation on the basis of spherical harmonics $Y^{m}_{l}(\theta,\phi)$. Because of the spherical symmetry of the background we take $m=0$ without breaking any generality Chandra . Therefore the basis is the Legendre polynomials $P_{l}(\theta)$.

It is well known that the perturbation can be decomposed into two kinds of classes according to their behavior under parity transformation. In this section, we will focus only on the even parity modes. For the even parity we focus on the static perturbations. Thus, the perturbations will have no explicit time dependence. After restricting ourselves in these conditions we choose the Regge-Wheeler gauge to fix the even parity perturbation $(\delta g^{(e)}_{\alpha\beta})$in the following form Regge1957 ,

[TABLE]

where, ${(e)}$ represents even parity sector.

III.1 Superfluid core

It is simple to calculate the perturbation in the energy momentum tensor. It can be expressed as, $\delta T^{0}_{0}=\delta\Lambda$ and $\delta T^{i}_{j}=\delta\Psi\delta^{i}_{j}$ . Using these in the Einstein equation and keeping only the first order of the perturbation, we can find the perturbed metric equations.

[TABLE]

[TABLE]

[TABLE]

where $\delta G^{r}_{r}=8\pi\delta T^{r}_{r}$ implies,

[TABLE]

From the linearized Euler equation we find,

[TABLE]

Staticity implies $\delta\mu_{0}=\delta\chi_{0}=0$. From [22] it is straightforward to show that,

[TABLE]

Using Eqs.(17), (22) and (23) we find,

[TABLE]

$\Lambda$ is a function of $n^{2},p^{2}$ and $x^{2}$. Therefore,

[TABLE]

where,

[TABLE]

We use the following Einstein equation along with the expression of $\delta\Lambda$ to calculate the final perturbation equation,

[TABLE]

After some calculation this reduces to,

[TABLE]

This is the central equation for the determination of the electric type tidal Love numbers. Note that Eq.(28) contains the coefficients $A_{\mu\nu}$, $B_{\mu\nu}$ and $C_{\mu\nu}$ which have been evaluated in the equilibrium configuration. The main difference between Eq. (28) and its non-superfluid single fluid counterpart Eq. (15) in Ref. Hinderer2008 is as follows. In the case of the normal fluid, it is assumed that the fluid is barotropic in nature. Therefore, it is possible to write $\delta\rho=\frac{d\rho}{dp}\delta p$ and substitute it in the perturbed Einstein equations. For any multi-fluid scenarios, this assumption is incorrect, in general. For this reason, we calculate $\delta\Lambda$ explicitly with respect to the fluid and the perturbed metric variables. Because of this, the final equation of even parity perturbation gets modified and so does the response to the perturbation subsequently.

III.2 Normal fluid envelope

We model the low density region as the one component normal fluid matter. Hence, it is simple to calculate the perturbation in the energy momentum tensor of the fluid. It can be expressed as, $\delta T^{0}_{0}=-\delta\rho=-\frac{d\rho}{dp}\delta p$ and $\delta T^{i}_{j}=\delta p\delta^{i}_{j}$ . Using these in the Einstein equation and keeping only the first order of the perturbation, perturbed metric equations have been found in several works Flanagan2008 ; Damour2009 . The equation is as follows,

[TABLE]

We take the initial condition for $H^{(l)}$ in the normal fluid region to be the value of the $H^{(l)}$ at the junction, found by solving Eq.(28). Then the solution of Eq.(29) gives the perturbation for the entire star.

IV Odd Parity Perturbation Equations for zero frequency mode

In this section we discuss the odd parity perturbation of the Einstein equation that will lead to the calculation of the magnetic type Love number. The zero frequency limit in the odd parity sector is discontinuous, as has been discussed in Ref. Pani_2018 . Keeping this in mind we take a time dependent perturbation of the metric and finally in the end we take the zero frequency limit carefully. After choosing the Regge-Wheeler gauge the metric perturbation $(\delta g^{(o)}_{\alpha\beta})$ can be written as follows,

[TABLE]

where $(o)$ represents odd parity.

IV.1 Superfluid core

For the odd parity modes $\delta n=0=\delta p$ where, $\delta p$ and $\delta n$ are the perturbed number density of the the proton and the neutron, respectively. If the perturbed velocity of the neutron and the proton are, respectively, $\delta u_{\mu}$ and $\delta v_{\mu}$ then only non-zero components can be written as Comer1999 ,

[TABLE]

where $U_{n}$ and $U_{p}$ are two arbitrary functions yet to be determined and $P_{l}$ is the Legendre polynomial.

Using the form of the velocity and metric perturbation in the Einstein equation, equation for the perturbations can be found. The equations relevant for our works are as follows:

[TABLE]

[TABLE]

A new master function is defined as, $\psi=e^{(\nu-\kappa)/2}\frac{h_{1}}{r}$ Pani_2018 . Equation (33) now can be written as,

[TABLE]

We take the time dependence of each mode as $e^{i\omega t}$. Putting everything together Eq.(32) can be written as,

[TABLE]

After taking the $\omega\rightarrow 0$ limit the zero frequency equation takes the following form,

[TABLE]

This is the central equation for the determination of the magnetic type tidal Love numbers. Note that Eq.(36) does not depend on the coefficients $A_{\mu\nu}$, $B_{\mu\nu}$ and $C_{\mu\nu}$ explicitly. But the effect of the SF nature enters through the dependence of $\Lambda$ on $x^{2}$. Because of this, the values of the magnetic Love numbers get modified even though the final equation of odd parity perturbation looks similar to the ones in Pani_2018 ; Damour2009 .

IV.2 Normal fluid envelope

Details of the odd parity equations for normal fluid can be found in Ref.Pani_2018 . The final equation is as follows,

[TABLE]

We take the initial condition for $\psi^{(l)}$ in the normal fluid region to be the value of the $\psi^{(l)}$ at the junction, found by solving Eq.(36). Then the numerical solution of Eq.(37) gives the solution for odd mode perturbation for the entire star.

V Calculation of the tidal Love numbers

V.1 Electric type Love numbers

To calculate the tidal deformability, we solve Eq. 28 numerically inside the NS up to the junction between the SF core and the NF envelope. Using the junction conditions described in Appendix [A.1] we find the initial condition of $H^{(l)}$ in the envelope. This initial condition has been used to numerically evolve Eq.(29) up to the surface of the NS. After that the tidal Love numbers are calculated by matching the numerical value of $H^{(l)}$ found by integration with the external solution of the same equation on the surface of the star. Extensive discussion on this can be found in Refs.Hinderer2008 ; Binnington2009 ; Damour2009 . Here we focus only on the initial conditions. We integrate Eq. 28 for metric perturbation in core $H^{(l)}$ radially outward from the center using the profiles of the background quantities calculated from TOV equations. For numerical purposes, instead of starting from $r=0$, we use a very small cutoff radius $(r=r_{0}=10^{-6})$. The initial condition for Eq. 28 around the regular singular point $r=0$ can be taken to be $H^{(l)}(r)\sim\bar{h}r^{l}$, with $\bar{h}$ some arbitrary constant. Since this equation is homogeneous and the tidal deformability depends explicitly on the value of $y^{even(l)}~{}(=\frac{rH^{(l)\prime}}{H^{(l)}})$ at the surface, the scaling constant $\bar{h}$ does not hold any relevance. Therefore, we can choose the starting value for the metric variable as, $H^{(l)}(r_{0})=r_{0}^{l}$ and $H^{\prime}(r_{0})=lr_{0}^{l-1}$.

The deformability is expressed in terms of $y^{even(l)}$, found by solving Eq.(29) in the envelope, and the compactness $C=\frac{M}{R}$, by matching the internal and external value of $H^{(l)}$ at the surface. The tidal Love numbers $k_{2}$ and $k_{3}$ then take the following functional form Hinderer2008 ; Binnington2009 ; Damour2009 ,

[TABLE]

[TABLE]

The expression for dimensionless deformability can be found from Damour et al. to be Damour2009 ,

[TABLE]

Since the information of the fluid enters through $y^{(l)}|_{r=R}$ and $C$, these expressions of $k_{2}$ and $k_{3}$ are similar to the one fluid formalism. Two-fluid formalism does not change the external solution. It only changes the internal equation of $H^{(l)}$, resulting in a different value of $y^{(l)}|_{r=R}$, leading to the change in the value of $k_{l}$ but not their expressions.

V.2 Magnetic type Love numbers

To calculate the magnetic type tidal deformability, we solve Eq. (36) numerically inside the NS up to the junction between the SF core and the normal fluid envelope. Then using the junction conditions described in Appendix [A.1] we find the initial condition of $\psi^{(l)}$ in the envelope. Using this initial condition we numerically evolve Eq.(37) up to the surface of the NS. The tidal Love numbers are calculated by matching the numerical value of $\psi^{(l)}$ found by integration with the external solution of the same equation on the surface of the star. Details can be found in Ref. Damour2009 . We will integrate Eq. 36 for $\psi^{(l)}$ radially outward from the center using the profiles of the background quantities calculated from the TOV equations. Similar to the calculations of the electric-type Love number, we start from a very small cutoff radius $(r=r_{0}=10^{-6})$. The initial condition for Eq. 36 near the regular singular point $r=0$ can be taken to be $\psi^{(l)}(r)\sim\bar{\psi}r^{l+1}$, with $\bar{\psi}^{(l)}$ some constant. Since, this equation is homogeneous in $\Psi^{(l)}$ and the tidal deformability depends explicitly on the value of $y^{odd(l)}~{}(=\frac{r\psi^{(l)\prime}}{\psi^{(l)}})$ at the surface, the scaling constant $\bar{\psi}^{(l)}$ is not relevant. Therefore, the starting value for the metric variable can be chosen as, $\psi^{(l)}(r_{0})=r_{0}^{l+1}$ and $\psi^{(l)\prime}(r_{0})=(l+1)r_{0}^{l}$.

The deformability can be expressed in terms of $y^{odd(l)}$, found by solving Eq.(37) in the envelope, and the compactness $C=\frac{M}{R}$, by matching the internal and external values of $\psi^{(l)}$ at the surface. The tidal Love number $j_{2}$ takes the functional form Damour2009 ,

[TABLE]

The expression for dimensionless deformability can be found from Damour et al. to be Damour2009 ,

[TABLE]

This expression of $j_{2}$ is similar to the one fluid formalism because the information of the fluid enters through $y^{(l)}|_{r=R}$ and $C$. The two-fluid model does not change the external solution. It changes only the internal equation of $\psi^{(l)}$, that gives us a different value of $y^{(l)}|_{r=R}$, leading to the change in the value of $j_{l}$ but not its expression.

VI Results

In this section, we discuss the numerical results for tidally deformed superfluid NS. At first, we calculate the static equilibrium configurations by solving the TOV equations using realistic EOS. Since only a few calculations are available for the two-fluid system in the literature, we choose a RMF type model with scalar self-interaction terms and use NL3 and GM1 parametrizations, as in paper I. We impose $\beta$-equilibrium at the center of the star by imposing $\mu|_{0}=\chi|_{0}$ to get a set of $k_{n}$, $k_{p}$ and $m_{*}$ for calculating the central number densities of the neutron and proton, energy density ($-\Lambda|_{0}$) and pressure ($\Psi|_{0}$). These quantities are used to solve Eqs. (9), (10), (14) and (15), to find the structure of the star and to generate profiles for various background quantities for several different sets of ($k_{n},k_{p},m_{*}$) that corresponds to the different central energy densities. The maximum mass, we have found to be 2.793 $\textnormal{M}_{\odot}$ for NL3 and the corresponding radius being 13.34 km. Similarly, for GM1, the maximum mass is calculated to be 2.384 $\textnormal{M}_{\odot}$ and the corresponding radius is 12.04 km. Details of those parameter sets can be found in Table 1. Moreover, for NL3 and GM1 sets, the crust-core transition pressures are 0.2698 and 0.2434 MeV$/$fm3 respectively. The two-fluid and the single fluid TOV integrations are smoothly joined at those pressures. Here, it is important to stress the fact that, these EOSs serve representative purposes only.

After getting the structure of the background, we find the numerical solution for $H^{(l)}$ for the entire star using Eqs. (28) and (29) and the junction conditions described in Appendix A.1. Using the background profiles mentioned earlier, find $y^{even(l)}$ at the surface of the stars and calculate the electric type Love numbers using Eqs. (38) and (39). Similarly we find the numerical solution for $\psi^{(l)}$ for the entire star using Eqs. (36) and (37) and the junction conditions described in Appendix A.2. Then we find $y^{odd(l)}$ at the surface of the stars and calculate the magnetic type Love number using Eq.(41). The behavior of $k_{2},k_{3}$ and $j_{2}$ w.r.t mass of the NS has been shown in the Figs.1,2 and 3 respectively, along with the case of normal fluid. We plot the dimensionless tidal deformabilities in Figs. 4, 5 and 6 along with the normal fluid case. The values of the tidal deformabilities for $1.4M_{\odot}$ is shown in Table 2. We show the percentage change in Fig. 7. For all the stellar configurations, we find the tidal deformabilities of the two fluid star are larger than the normal one fluid stars. To calculate the tidal deformabilities for the normal fluid case we used the unified EOS. As a result in both cases of NL3 and GM1, the crust is included in the calculation.

It is important to note that when we speak of the deviation of $\Lambda_{2}$ due to the superfluid nature, we bring an ambiguity in our interpretation of the observed $\Lambda_{2}$. The value of $\Lambda_{2}$ in the two-fluid calculation for a particular EOS model can be similar to the value in a single-fluid calculation for another EOS. So, we can not distinguish between the EOS and also probe the fluid nature of matter at the same time with the measurement of $\Lambda_{2}$. There are other possible degeneracies that can affect its value too Biswas2019anisotropy ; Raposo2019 . In Sec. VII we discuss how this degeneracy can be broken.

VII Universal relation

In this section we fit compactness $C$, $\Lambda^{el}_{3}$ and $\Lambda^{mag}_{2}$ calculated in the previous sections against $\Lambda^{el}_{2}$, to test the universal relation. In Fig. 8 we plot $C$ against $\Lambda^{el}_{2}$. The upper half of the left panel represents the case when the fluid has been taken to be one component normal fluid. The upper half of the right panel represents the case when the matter is modeled as a two component superfluid core and a normal fluid envelope. For all the cases we fit them with a fitting function. The lower halves of both panels show errors w.r.t. fitted curves. In Fig. 9 we plot $\Lambda^{el}_{3}$ against $\Lambda^{el}_{2}$. The upper half of the left panel represents the case when the fluid has been taken to be one component normal fluid. The upper half of the right panel represents the case when the matter is modeled as a two component superfluid core and a normal fluid envelope. For all cases we fit them with a fitting function. The lower halves of the both panels show errors w.r.t. fitted curves. For the $C-\Lambda^{el}_{2}$ relation, we fit the results for both the one fluid and the SF case with the following function Yagi_2014 ; Maselli_2013 ,

[TABLE]

For the other cases, we used the following fitting function Yagi_2014 ; Maselli_2013 ,

[TABLE]

The details of the fitted values of the parameters are described in Table 3. In Fig. 10 we show the differences between the $C-\Lambda^{el}_{2}$ universality curves for the different fluid scenario. For comparison, we have plotted the corresponding curve using the fitting parameters from Ref. Maselli_2013 , which has been named as $"1-fl$ Maselli$+"$. In Fig. 11 we show the differences between the $\Lambda^{el}_{3}-\Lambda^{el}_{2}$ universality curves for different fluid scenarios. For comparison, we have plotted the corresponding curve using the fitting parameters from Ref. Yagi_2014 , which has been named as $"1-fl$ Yagi$"$.

Interestingly one fluid formalism and SF formalism both show universal behavior, even though the values of tidal deformabilities change due to the inclusion of the SF. But a crucial feature that we have found is that the fitted curve for normal fluid formalism is different from the scenario when the matter is treated using two fluid formalism. For example in Table 3 it can be seen that the values of $a,b,c,d,f$ are different for two different formalisms. This is an important observation as it opens up the possibility to probe the SF nature of matter using this deviation. From the GW data, it is possible to estimate the values of $\Lambda^{el}_{2}$ and $\Lambda^{el}_{3}$. Having an estimation of such manner it is possible to check which universal relation is more suitable for the observed values of the deformabilities. As the two different fluid natures imply different universal relations, measured values of the deformabilities will be able to distinguish between the two different universality curves. This result will help us break the degeneracy between the fluid nature and the EOS, discussed in earlier sections. It is important to note that even though the universal curves are different, they are not “too different”. Therefore it remains to be seen whether this strategy will be able to break the degeneracy with the real data.

VIII Conclusion

Results found in the current work are very important in the context of constraining the dense matter EOS using the GW data. Values of the deformabilities for superfluid NS are higher than the normal fluid star, for a given RMF model. At present tight constraint has been put on the EOS from the BNS observation Abbott2017 ; Abbott2018 ; Nandi_2018 ; Zack_2019 ; Zack_2019_2 ; Bharat_2019 ; Luca_2019 . The results found here indicate that, more EOSs will be ruled out which are otherwise allowed if we do not consider superfluidity inside the NS. This provides us the opportunity to improve our understanding of the SF nature of the dense matter with better observational data in the future.

We find that the Love numbers are usually larger for a two-fluid system. Comer et al. Comer1999 found the existence of several superfluid oscillation modes that cannot be found otherwise in a single fluid star. This nature is very specific to the two-fluid formalism where different fluid modes can appear due to the existence of the two different types of fluid displacements. Flanagan and Hinderer discussed the fact that the tidal deformation of a star can be thought of as the sum of the deformations arising from different fluid modes that have been excited inside the star, due to the tidal perturbation Flanagan2008 . Therefore, we can say that, due to the appearance of extra fluid modes in the superfluid stars, we will get slightly larger deformations under tidal perturbation.

We argued that there is a degeneracy between the fluid nature and the EOS. Interestingly, we found that in the SF case the tidal deformabilities show a universal relation but the universal curve is different from the one fluid case. We discussed how measuring different tidal deformabilities and using universal relations can break the degeneracy between the fluid nature and the EOS.

Acknowledgements

We thank Paolo Pani, Bhaskar Biswas, Niels Andersson, Rana Nandi and Sukanta Bose for helpful discussions. We are also grateful to H. Pais for kindly providing us with the GM1 inner crust table. This work is supported in part by the Navajbai Ratan Tata Trust. S.D. would like to thank the University Grants Commission (UGC), India, for financial support.

Appendix A Junction condition

In our current work we have modeled the NS as a superfluid core with a normal fluid envelope. Crustal physics is encoded in the current model via the normal fluid envelope. As there are two layers of fluid in our model, it is necessary to find the junction condition across the boundary. For the purpose of simplicity in this section we will use $\Psi$ as the symbol for pressure both in the SF core and in the NF envelope, while we derive the junction conditions. To calculate the junction conditions we take the level surfaces of $\Psi$. As there are no “delta-function like” discontinuities in $\Psi$, the first and second fundamental forms are continuous everywhere inside the star MTW . Therefore, by imposing continuity in the first and second forms we can find the junction conditions.

The normal to the level surface of $\Psi$ is

[TABLE]

The induced three metric (first fundamental form) $\gamma_{\mu\nu}$ is,

[TABLE]

The extrinsic curvature (second fundamental form) $\mathcal{K}_{\mu\nu}$ is defined as follows,

[TABLE]

where parentheses imply symmetrization of the indices. Junction conditions will be found from the continuity of $\gamma_{\mu\nu}$ and $\mathcal{K}_{\mu\nu}$.

A.1 Equilibrium configuration and even parity sector

As we are mainly interested in the perturbation on the background, we write $\Psi$ as follows,

[TABLE]

As a smooth background is constructible even in the presence of perturbation, we assume that the background and the perturbed part of $\gamma_{\mu\nu}$ and $K_{\mu\nu}$ are separately continuous at the junction. We will discuss only those components of $\gamma_{\mu\nu}$ and $K_{\mu\nu}$ that are relevant for our purpose, for more details see Ref.Andersson_2002 . First we consider the components that are useful for the even mode perturbation added to the background quantities. In the zero frequency limit the relevant quantities can be expressed as,

[TABLE]

With $\delta\Psi(r,\theta)=\delta\Psi(r)P_{l}(\theta)$ these sets of equations imply,

[TABLE]

where, $R_{c}$ represents the radius of the boundary. Physical quantities with no tilde represent their value in the SF region just below the junction. A tilde represents the value of the physical quantity in the normal fluid region just above the junction.

A.2 Odd parity sector

For the continuity of the quantities of the odd mode perturbation we follow similar procedure. But as has been discussed earlier we consider the time dependent perturbation for that purpose. we find

[TABLE]

Taking $h_{i}(t,r)=\int d\omega\hat{h}_{i}(\omega,r)e^{-i\omega t}$ implies $\hat{h}_{i}$ is continuous implying $\psi$ is continuous (for the definition check IV). Continuity of $K_{03}$ implies $\omega e^{(\nu-\kappa)/2}r\psi+\frac{\{e^{(\nu-\kappa)/2}(\psi r)^{\prime}\}^{\prime}}{\omega}$ is continuous. Using Eq.(35) in SF region we find that the following expression is continuous:

[TABLE]

A similar expression can be found in the normal fluid region with $\Psi\rightarrow p$ and $\Lambda\rightarrow-\rho$. Since $\psi$ is continuous, this implies $\psi^{\prime}$ is continuous across the junction.

Appendix B Expressions for matter variables

In the limit $K\rightarrow 0$ the master function and the chemical potentials of the neutron and proton fluids can be expressed as,

[TABLE]

The generalized pressure $\Psi$ and the master function are realized by the following relationship,

[TABLE]

In the above expressions, $c_{\sigma}^{2}=(g_{\sigma}/m_{\sigma})^{2}$, $c_{\omega}^{2}=(g_{\omega}/m_{\omega})^{2}$, $c_{\rho}^{2}=(g_{\rho}/m_{\rho})^{2}$ and

[TABLE]

The expressions for the other matter coefficients (see Kheto2014 ; Kheto2015 ) that are used as the inputs in field equations are as follows:

[TABLE]

[TABLE]

where,

[TABLE]

respectively.