Temperature-dependent r-modes in superfluid neutron stars stratified by muons

TL;DR

This paper investigates how superfluidity and muons in neutron star cores influence the spectrum and stability of r-modes, revealing temperature-dependent mode interactions and potential suppression of instabilities.

Contribution

It introduces a detailed calculation of temperature-dependent superfluid r-modes in neutron stars with muons, highlighting the impact of muons on the mode spectrum and stability.

Findings

Existence of two main r-modes in minimal NS cores: normal and superfluid.

Muons cause an infinite set of superfluid r-modes with temperature-dependent frequencies.

Normal r-mode can strongly dissipate at resonance temperatures, suppressing instabilities.

Abstract

We calculate the finite-temperature r-mode spectrum of a slowly rotating superfluid Newtonian neutron star neglecting the entrainment between neutron and proton liquid components (i.e., neglecting the off-diagonal element of the entrainment matrix). We show that for `minimal' NS core composition (neutrons, protons, and electrons) only two $m=2$ r-modes exist --- normal mode, which is similar to ordinary r-mode in a nonsuperfluid star, and a superfluid temperature-dependent mode. Accounting for muons in the core dramatically modifies the oscillation spectrum, resulting in an infinite set of superfluid r-modes, whose frequencies vary with temperature. We demonstrate that the normal r-mode can exhibit avoided crossings with superfluid modes at certain `resonance' temperatures, where it dissipates strongly, which leads to substantial suppression of the r-mode instability near these temperatures. The corresponding instability windows are calculated and discussed.

Equations

$$T_{;\nu }^{\mu \nu }=0$$

$$T^{\mu \nu }=(P+ϵ)u^{\mu }{u}^{\nu }+P{g}^{\mu \nu }+Y_{ik}\left(w_{(i)}^{\mu }{w}_{(k)}^{\nu }+{\mu }_i{w}_{(k)}^{\mu }{u}^{\nu }+{\mu }_k{w}_{(i)}^{\nu }{u}^{\mu }\right)+{\tau }^{\mu \nu },$$

$${\partial }_{\mu }{j}_{(j)}^{\mu }=0,$$

$$j_{(k)}^{\mu }=n_k{u}^{\mu }+Y_{ik}{w}_{(i)}^{\mu },$$

$$j_{(\mathrm{l})}^{\mu }=n_{\mathrm{l}}{u}^{\mu }$$

$$u_{\nu }{\mathcal{V}}^{\mu \nu }={\mu }_{\mathrm{n}}{n}_{\mathrm{n}}{f}^{\mu }$$

$$f^{\mu }=\alpha {\perp }^{\mu \nu }{\mathcal{V}}_{\nu \lambda }\frac{1}{n_{\mathrm{n}}}{Y}_{\mathrm{n}k}{w}_{(k)\delta }{\perp }^{\lambda \delta }+\frac{{\beta }_{\mathrm{MF}}}{{\mathcal{V}}_{(\mathrm{M})}}{\perp }^{\mu \eta }{\perp }^{\nu \sigma }{\mathcal{V}}_{\eta \sigma }{\mathcal{V}}_{\lambda \nu }\frac{1}{n_{\mathrm{n}}}{Y}_{\mathrm{n}k}{w}_{(k)\delta }{\perp }^{\lambda \delta },$$

$$u_{\mu }{w}_{(i)}^{\mu }=0,$$

$$dϵ=TdS+\sum _{j=\mathrm{n},\mathrm{p},\mathrm{e},\mu }{\mu }_{j}d{n}_j+\frac{Y_{ik}}{2}d\left(w_{(i)}^{\alpha }{w}_{(k)\alpha }\right),$$

$$P=-ϵ+\sum _{j=\mathrm{n},\mathrm{p},\mathrm{e},\mu }{\mu }_j{n}_j+TS,$$

$$n_{\mathrm{p}}=n_{\mathrm{e}}+n_{\mu },$$

$$Y_{\mathrm{p}k}{w}_{(k)}^{\mu }=0,$$

$$-{\sigma }^2{\mathbf{\xi }}_{\mathrm{b}}+2ı\sigma \mathbf{\Omega }\times {\mathbf{\xi }}_{\mathrm{b}}=\frac{\delta w}{w_0^2}\nabla P_0-\frac{\nabla \delta P}{w_0},$$

$${\mathbf{\xi }}_{\mathrm{b}}\equiv \frac{{\mathbf{j}}_{\mathrm{b}}}{ı\sigma n_{\mathrm{b}}},$$

$$\delta n_{\mathrm{b}}+\mathrm{div}(n_{\mathrm{b}}{\mathbf{\xi }}_{\mathrm{b}})=0,$$

$$\delta n_{\mathrm{l}}+\mathrm{div}(n_{\mathrm{l}}\mathbf{\xi })=0.$$

$$h{\sigma }^2\mathbf{z}+2ı\sigma \mathbf{\Omega }\times \mathbf{z}{n}_{\mathrm{b}}{\mu }_{\mathrm{n}}\left(1+\alpha \frac{n_{\mathrm{b}}{\mu }_{\mathrm{n}}}{c^2}\right)=h{\sigma }^2\mathbf{z}-2ıh_1\sigma \mathbf{\Omega }\times \mathbf{z}=c^2{n}_{\mathrm{e}}\nabla \Delta {\mu }_{\mathrm{e}}+c^2{n}_{\mu }\nabla \Delta {\mu }_{\mu },$$

$$h=n_{\mathrm{b}}{\mu }_{\mathrm{n}}y,$$

$$h_1={\mu }_{\mathrm{n}}{n}_{\mathrm{b}}\left(\frac{n_{\mathrm{b}}}{Y_{\mathrm{n}k}{\mu }_k}-1\right),$$

$$y=\frac{n_{\mathrm{b}}{Y}_{\mathrm{pp}}}{{\mu }_{\mathrm{n}}(Y_{\mathrm{nn}}{Y}_{\mathrm{pp}}-Y_{\mathrm{np}}^2)}-1.$$

$$\delta n_i=\frac{\partial n_i}{\partial P}\delta P+\frac{\partial n_i}{\partial \Delta {\mu }_{\mathrm{e}}}\Delta {\mu }_{\mathrm{e}}+\frac{\partial n_i}{\partial \Delta {\mu }_{\mu }}\Delta {\mu }_{\mu }.$$

$$r=x\left[1-{\Omega }^2\alpha (x){\cos }^2\theta \right],$$

$$\sigma =\Omega {\sigma }_0\left(1+{\Omega }^2{\sigma }_1\right),$$

$$\delta f={\Omega }^2\delta f_1\exp (ı\sigma t+ım\varphi ).$$

$$d_r={\Omega }^2{d}_r^1\exp (ı\sigma t+ım\varphi ),$$

$$d_{\theta }=(d_{\theta }^0+{\Omega }^2{d}_{\theta }^1)\exp (ı\sigma t+ım\varphi ),$$

$$d_{\varphi }=(d_{\varphi }^0+{\Omega }^2{d}_{\varphi }^1)\exp (ı\sigma t+ım\varphi ),$$

$$\frac{\partial }{\partial \theta }\sin \theta {\xi }_{\mathrm{b}\theta }^0+ım{\xi }_{\mathrm{b}\varphi }^0=0,$$

$$\frac{\partial }{\partial \theta }\sin \theta {\xi }_{\theta }^0+ım{\xi }_{\varphi }^0=0.$$

$$\frac{\partial }{\partial \theta }\sin \theta z_{\theta }^0+ım{z}_{\varphi }^0=0.$$

$${\sigma }_0{\xi }_{\mathrm{b}\theta }^0+2ı\cos \theta {\xi }_{\mathrm{b}\varphi }^0=-\frac{ı}{m}\frac{\partial }{\partial \theta }\sin \theta \left({\sigma }_0{\xi }_{\mathrm{b}\varphi }^0-2ı\cos \theta {\xi }_{\mathrm{b}\theta }^0\right),$$

$$h(x){\sigma }_0{z}_{\theta }^0+2ıh_1(x)\cos \theta z_{\varphi }^0=-\frac{ı}{m}\frac{\partial }{\partial \theta }\sin \theta \left[h(x){\sigma }_0{z}_{\varphi }^0-2ıh_1(x)\cos \theta z_{\theta }^0\right].$$

$${\sigma }_0=\frac{2m}{l(l+1)},$$

$${\sigma }_0=\frac{2m}{l(l+1)}\frac{h_1(x)}{h(x)}$$

$${\xi }_{\mathrm{b}\theta }^0=ım{C}_{lm}(x)\frac{P_l^m(\cos \theta )}{\sin \theta },{\xi }_{\mathrm{b}\varphi }^0=-C_{lm}(x)\frac{d}{d\theta }{P}_l^m(\cos \theta ),$$

Full text

††pagerange: Temperature-dependent r-modes in superfluid neutron stars stratified by muons–References††pubyear: 2017

1 Introduction

Since it had been shown that r-modes in neutron stars (NSs) can be unstable (Andersson 1998; Friedman & Morsink 1998) and can effectively spin NSs down by transmitting their angular momentum to gravitational waves, a significant progress has been made in understanding the r-mode physics (Levin 1999; Lindblom & Mendell 2000; Andersson & Comer 2001; Andersson & Kokkotas 2001; Heyl 2002; Lee & Yoshida 2003; Ho, Andersson, & Haskell 2011; Haskell, Degenaar, & Ho 2012; Mahmoodifar & Strohmayer 2013; Alford & Schwenzer 2014; Gusakov, Chugunov, & Kantor 2014a, b; Haskell 2015). However, a number of important questions related to r-modes has remained unanswered. The main problem concerns the interpretation of observations of hot rapidly rotating NSs in low-mass X-ray binary (LMXB) systems. According to theoretical predictions the probability to observe such objects should be negligible since they should be heated and spun down by developing r-mode instability on a time-scale much smaller than the typical age of an LMXB system (Levin 1999). A number of possible solutions to this problem has been proposed (see, e.g., Haskell 2015 for a recent review). One of them (Gusakov et al. 2014a, b) appeals to baryon superfluidity in the NS core. The idea is that the instability may be damped out at some stellar temperatures by the resonance interaction of the normal r-mode with superfluid inertial modes 111Following Gusakov et al. (2014a), by the normal r-modes of a superfluid NS we mean the modes which, at a given temperature, have properties similar to ordinary r-modes in the non-superfluid NSs. They correspond to co-motion of superfluid and normal liquid components, and their eigenfrequencies do not depend on temperature. Normal r-modes are most unstable. By superfluid inertial modes we understand inertial modes which show superfluid-like behaviour at a given temperature. The oscillations of such modes are mostly counter-moving, baryon current is almost zero, their eigenfrequencies depend on temperature. These modes have no analogue in non-superfluid NSs. Note that, with temperature variation, a given mode can switch from a ‘superfluid’ regime to ‘normal’ regime and vice versa, experiencing avoided-crossings with neighbouring modes. in the vicinity of their avoided crossings. This mechanism was previously analysed qualitatively, but not quantitatively (Gusakov et al. 2014a, b; Chugunov et al. 2014; Kantor et al. 2016; Gusakov et al. 2016). To study the resonance interaction of normal r-mode with superfluid modes one needs to calculate the temperature-dependent oscillation spectrum of a rotating superfluid neutron star. This has not been done so far since most of the studies have been concentrated on the limiting case of vanishing stellar temperature (e.g., Lindblom & Mendell 2000; Andersson & Comer 2001; Andersson & Kokkotas 2001; Lee & Yoshida 2003; Andersson et al. 2009; Haskell et al. 2009).

In this paper we present the first calculation of the temperature-dependent spectrum of a rotating superfluid NS. The calculations are performed for two NS models: The first one is standard (see, e.g., Lee & Yoshida 2003) and assumes minimal composition of the NS core (neutrons, protons, and electrons), while the second one accounts for admixture of muons. The presence of muons ‘stratifies’ the NS core (Kantor & Gusakov 2014; Dommes & Gusakov 2016; Passamonti, Andersson, & Ho 2016), which, as we show, has a dramatic effect on the low-frequency spectrum of rotating NSs. Using these results, we study the effect of resonance interaction of normal r-mode with superfluid modes on the stability of NSs and plot the ‘instability windows’ (the regions of stellar rotation frequency and temperature, in which NSs are unstable, see, e.g., Andersson & Kokkotas 2001) accounting for such interaction.

To simplify calculations we shall work in the Newtonian limit of the relativistic hydrodynamics, describing superfluid mixtures at a finite temperature (see, e.g., Mendell & Lindblom 1991; Gusakov & Andersson 2006; Gusakov, Kantor, Chugunov, & Gualtieri 2013; Gusakov 2016; Gusakov & Dommes 2016). We also use a number of approximations. First, we adopt the Cowling approximation (i.e., gravitational potential is not perturbed, Cowling 1941). Second, we consider only slowly rotating NSs, expanding the corresponding eigenmodes in the Taylor series in the rotation frequency, and allowing only for the first two terms in the expansion. Third, we assume that there is no entrainment, i.e., we neglect the off-diagonal element $Y_{\rm np}$ of the entrainment matrix (see Andreev & Bashkin 1976; Gusakov & Haensel 2005; Gusakov & Andersson 2006; Gusakov, Kantor, & Haensel 2009b, a, 2014c), $Y_{\rm np}=0$. This assumption is not unreasonable since numerical calculations show that entrainment is indeed small (Gusakov et al. 2009b, a, 2014c). Moreover, in the region where protons are non-superfluid, $Y_{\rm np}$ is strictly zero (see, e.g., Gusakov et al. 2014c). We plan to relax this simplifying assumption in the subsequent publication. And finally, we calculate only those modes which are purely toroidal to the leading order in rotation.

The paper is organized as follows. In section 2 we briefly discuss the hydrodynamics of superfluid mixtures (section 2.1) and apply it to derive the system of equations describing NS toroidal oscillations (section 2.2). Together with the boundary conditions (section 3), this system is solved numerically for a physics input from section 4. The numerical results are presented in section 5 for two models of the NS core with $\rm npe$ (section 5.1) and $\rm npe\mu$ compositions (section 5.2). The oscillation spectra for these two models are compared and interpreted in section 6. In section 7 we apply the calculated spectra to study damping of r-modes in the presence of avoided crossings of modes and plot the corresponding instability windows. Finally, in section 8 we discuss our results and try to understand how the oscillation spectrum will be modified in a more realistic situation of non-zero entrainment.

The paper also contains four appendices, which include some details on the derivation of the oscillation equations (Appendix A), an analytical solution for the toroidal modes of a superfluid NS with $\rm npe$ core composition (Appendix B), the expression for the mutual friction coefficient (Appendix C), and a list of typos in the illuminating paper by Provost, Berthomieu, & Rocca (1981), whose analysis of r-modes in normal stars is similar to that presented here (Appendix D).

2 Equations governing oscillations of a rotating NS

2.1 General equations of superfluid hydrodynamics

In our analysis we shall use the Newtonian limit of the relativistic hydrodynamics, formulated in Gusakov (2016); Gusakov & Dommes (2016) and describing rotating superfluid mixture at a finite temperature. Below we formulate the main equations of this hydrodynamics. We assume that the mixture is composed of nucleons (neutrons and protons), which can be in superfluid state, and of normal leptons (electrons and muons). We ignore all the dissipative effects except for mutual friction, and shear viscosity, which turn out to be important for calculating the instability windows in section 7. We also neglect the vortex-related contribution to the energy density, which is an accurate approximation for our problem (Sauls 1989), and do not consider the magnetic field effects. The equations governing superfluid dynamics consist of:

(i) Energy-momentum conservation law

[TABLE]

with the energy-momentum tensor

[TABLE]

where $P$ is the pressure, $\epsilon$ is the energy density; $g^{\mu\nu}$ is the metric tensor; small dissipative terms due to shear viscosity are contained in the tensor $\tau^{\mu\nu}$, see Gusakov et al. (2013) for details; $Y_{ik}$ is the relativistic symmetric entrainment matrix (Gusakov & Andersson 2006; Gusakov et al. 2009b, a, 2014c), analogue of the superfluid or mass-density matrix $\rho_{ik}$ of the non-relativistic theory (Andreev & Bashkin 1976; Borumand et al. 1996; Gusakov & Haensel 2005; Chamel & Haensel 2006; Gusakov 2010). Here and below indices $i,k$ run over neutrons and protons: $i$, $k={\rm n}$, ${\rm p}$. Note that, in Eq. (2) and further in the text we assume summation over the repeated indices $i$ and $k$. Finally, $u^{\mu}$ is the four-velocity of the non-superfluid component (leptons and baryon thermal excitations), normalized by the condition $u_{\mu}u^{\mu}=-1$, and $w^{\nu}_{(k)}$ ($k={\rm n},{\rm p}$) is the four-vector that describes superfluid degrees of freedom (in particular, it is related to the superfluid velocity $v_{({\rm s}k)}^{\nu}$ by $v_{({\rm s}k)}^{\nu}=(w^{\nu}_{(k)}+\mu_{k}u^{\nu})/m_{k}$, where $\mu_{k}$ and $m_{k}$ are, respectively, the relativistic chemical potential and the bare mass of particle species $k$).

(ii) Continuity equations for all particle species $j$

[TABLE]

where the particle current density for neutrons and protons is

[TABLE]

and for leptons (electrons and muons)

[TABLE]

(here and below ${\rm l}={\rm e}$, ${\rm\mu}$). In equations (4) and (5) $n_{j}$ is the number density for particle species $j$.

(iii) Neutron ‘superfluid’ equation

[TABLE]

where $\mathcal{V}^{\mu\nu}\equiv\partial^{\mu}(w_{({\rm n})}^{\nu}+\mu_{\rm n}u^{\nu})-\partial^{\nu}(w_{({\rm n})}^{\mu}+\mu_{\rm n}u^{\mu})$, and $f^{\mu}$ is given by the expression

[TABLE]

in which $\perp^{\mu\nu}\equiv g^{\mu\nu}+u^{\mu}u^{\nu}$, $\mathcal{V}_{({\rm M})}\equiv\sqrt{\mathcal{V}^{\mu}_{({\rm M})}\mathcal{V}_{({\rm M})\mu}}$, $\mathcal{V}^{\mu}_{({\rm M})}\equiv\frac{1}{2}\,\epsilon^{\mu\nu\gamma\delta}\,u_{\nu}\,\mathcal{V}_{\gamma\delta}$ (here $\epsilon^{\mu\nu\gamma\delta}$ is the Levi-Civita tensor, $\epsilon^{0123}=1$). Kinetic coefficients $\alpha$ and $\beta_{\rm MF}$ depend on the details of interaction between the non-superfluid component and vortices. In what follows we will be interested in the so-called weak-drag regime, when vortices weakly interact with the normal component (this regime takes place in the NS cores, see, e.g., Mendell 1991; Andersson et al. 2006). In this limit $\alpha=-1/(\mu_{\rm n}\mu_{k}Y_{{\rm n}k})$, while the coefficient $\beta_{\rm MF}$ is small and describes dissipation due to mutual friction (the dissipative mechanism related to scattering of electrons off the vortices).

These equations should be supplemented by the condition ensuring that all the thermodynamic quantities are defined in the comoving frame in which $u^{\mu}=(1,\,0,\,0,\,0)$,

[TABLE]

as well as by the second law of thermodynamics,

[TABLE]

the expression for pressure,

[TABLE]

and the conditions

[TABLE]

which are always satisfied (Mendell 1991) for the low-frequency hydrodynamic oscillations we are interested in; these conditions indicate that protons are effectively locked to electrons and muons by the electromagnetic forces. In equations (9) and (10) $T$ is the temperature and $S$ is the entropy density.

2.2 Oscillation equations

Using the Newtonian limit of the above equations, let us consider small oscillations of a slowly rotating (with the spin frequency $\Omega$) non-dissipative Newtonian NS in the Cowling approximation. In what follows we shall allow for muons in the inner layers of NSs, assuming $\rm npe\mu$-composition, and also take into account possible superfluidity of baryons (neutrons and protons) in the core. Let all the quantities depend on time $t$ as ${\rm e}^{\imath\sigma t}$ in the coordinate frame rotating with the star. Then the linearised equations governing small oscillations of superfluid (hereafter SFL) NSs in that frame consist of:

(i) Euler equation

[TABLE]

where $w=(P+\epsilon)/c^{2}$, $c$ is speed of light. Here and hereafter, the subscript [math] denotes the equilibrium value of some quantity (e.g., $P_{0}$) and $\delta$ stands for its Euler perturbation (e.g., $\delta P$). Finally, ${\boldsymbol{\xi}}_{\rm b}$ in equation (13) is the Lagrangian displacement of baryons, it is defined as

[TABLE]

where $n_{\rm b}\equiv n_{\rm n}+n_{\rm p}$ and ${\boldsymbol{j}}_{\rm b}\equiv{\boldsymbol{j}}_{\rm n}+{\boldsymbol{j}}_{\rm p}$ are the baryon number density and baryon current density, respectively.

(ii) Continuity equations for baryons and leptons (electrons and muons)

[TABLE]

Here ${\boldsymbol{\xi}}\equiv{\boldsymbol{j}}_{\rm e}/(\imath\sigma n_{\rm e})$ is the Lagrangian displacement of the normal liquid component [we assume that all the normal-matter constituents (i.e., leptons and baryon thermal excitations) move with one and the same normal velocity due to efficient particle collisions]. If neutrons are non-superfluid, then ${\boldsymbol{\xi}}={\boldsymbol{\xi}}_{\rm b}$ and hydrodynamic equations become essentially the same as in the normal matter (even if protons are SFL, see equations 11 and 12), see e.g. Gusakov & Andersson 2006. Taking this into account, we, for brevity, shall call ‘normal’ (or ‘non-superfluid’) the liquid with non-superfluid neutrons, irrespective of the actual state of protons.

(iii) The ‘superfluid’ equation, analogue of the Euler equation for superfluid (neutron) liquid component

[TABLE]

where ${\boldsymbol{z}}\equiv{\boldsymbol{\xi}}_{\rm b}-{\boldsymbol{\xi}}$ is the superfluid Lagrangian displacement; and $\Delta\mu_{\rm l}\equiv\mu_{\rm n}-\mu_{\rm p}-\mu_{\rm l}$ is the chemical potential imbalance (note that in equilibrium $\Delta\mu_{\rm l}=0$, see Haensel et al. 2007, thus $\delta\Delta\mu_{\rm l}=\Delta\mu_{\rm l}$). Further,

[TABLE]

Note that in equation (19) summation over $k={\rm n}$, ${\rm p}$ is assumed. Superfluid equation in the form (17) is valid in the weak-drag regime only, when the interaction between the neutron vortices and normal component (e.g., electrons) is weak, which is a typical situation in NSs (see, e.g., Mendell 1991; Andersson et al. 2006). In addition, equation (17) assumes that Newtonian limit is justified. Since in that limit the redshift ${\rm e}^{\nu/2}\approx 1$ and $\nu^{\prime}/2=\nabla P_{0}/(w_{0}c^{2})\ll 1/R$ ($R$ is the stellar radius, which is a typical length-scale of the problem; we also used the fact that $P_{0}\ll\epsilon_{0}$ in the Newtonian limit), we skipped all redshifts in equation (17), replacing, e.g., $\Delta\mu_{\rm l}^{\infty}\equiv\Delta\mu_{\rm l}{\rm e}^{\nu/2}$ with $\Delta\mu_{\rm l}$ in the right-hand side of this equation. Moreover, the Newtonian limit implies that one can neglect (i) the frame drag effect, which is $\sim\Omega_{0}^{2}R^{3}\Omega/c^{3}$, where $\Omega_{0}\equiv\left(GM/R^{3}\right)^{1/2}$ is of the order of the Kepler frequency (and $M$ is the stellar mass); (ii) all the terms $\sim\Omega^{2}R^{2}/c^{2}$, since $c\rightarrow\infty$ in the Newtonian limit. Although for NSs it is not a well justified approximation (because they are essentially relativistic objects), however, it generally gives qualitatively correct results (see, e.g. Idrisy, Owen, & Jones 2015).222A detailed discussion of r-modes in relativistic neutron stars and the related problems (e.g., the problem of a continuous r-mode spectrum in non-barotropic relativistic stars) can be found in Kojima (1998); Kojima & Hosonuma (1999); Lockitch, Andersson, & Friedman (2001, 2003); Yoshida & Lee (2002, 2003); Lockitch, Andersson, & Watts (2004). The above equations should be supplemented with the ‘equation of state’ (EOS),

[TABLE]

In what follows we shall use $P$, $\Delta\mu_{\rm e}$ and $\Delta\mu_{\rm\mu}$ as independent thermodynamic variables. We shall consider a slowly rotating NS with $\Omega\ll\Omega_{0}$, and shall expand all the relevant quantities in a power series in small parameter $\Omega/\Omega_{0}$, keeping the first two terms in the expansion. This means, that the NS rotational oblateness, that has an order $(\Omega/\Omega_{0})^{2}$, should be accounted for. The isobaric surfaces of a rotating NS are not spherical as in non-rotating NSs. When the rotation is slow, the coordinates $r$ and $\theta$ of the isobaric surfaces are related by (Chandrasekhar 1933; Chandrasekhar & Roberts 1963; Hartle 1967; Provost, Berthomieu, & Rocca 1981; Saio 1982)

[TABLE]

where $r$ and $\theta$ are the radial coordinate and polar angle in the spherical coordinate system with the origin at the stellar centre and the axis z parallel to ${\boldsymbol{\Omega}}$. Every isobaric surface is characterized by the $\theta$-independent parameter $x$, $x={\rm const}$. The function $\alpha(x)$ (not to be confused with the mutual friction coefficient $\alpha$, see equation 7) is to be determined from the hydrostatic equilibrium equations. To calculate $\alpha(x)$ we used the Hartle scheme (Hartle 1967; Hartle & Thorne 1968). It allows us to find the structure of the relativistic slowly rotating star up to (and including) the terms $\sim(\Omega/\Omega_{0})^{2}$. In what follows it turns out to be more convenient to work in the coordinates $x$, $\theta$ and $\phi$ instead of $r$, $\theta$ and $\phi$ ($\phi$ is the azimuthal angle). Then the oblateness appears in the oscillation equations through the function $\alpha(x)$.

In the present paper we are interested in the small-amplitude (linear) oscillations, which depend on time $t$ and azimuthal angle $\phi$ as ${\rm e}^{\imath\sigma t+\imath m\phi}$ in the frame rotating with the star, and which have the eigenfrequencies $\sigma$ vanishing at $\Omega\rightarrow 0$. Thus, up to the terms $\sim(\Omega/\Omega_{0})^{2}$, $\sigma$ and the Euler perturbation of any (scalar) thermodynamic parameter $f$ (e.g., $P$, $\mu_{\rm l}$, $n_{\rm b}$, etc.) can be presented as (e.g., Provost, Berthomieu, & Rocca 1981; Lockitch & Friedman 1999; Lindblom & Mendell 2000)

[TABLE]

Concerning Lagrangian displacements ${\boldsymbol{\xi}}_{\rm b}$, $\boldsymbol{\xi}$, and $\boldsymbol{z}$, we look for a purely toroidal (to leading order in $\Omega/\Omega_{0}$) oscillation modes which assume the following ordering (Lockitch & Friedman 1999)

[TABLE]

where ${\boldsymbol{d}}$ stands for the displacements ${\boldsymbol{\xi}}_{\rm b}$, ${\boldsymbol{\xi}}$, or ${\boldsymbol{z}}$. Then the continuity equations (15) and (16) for baryons ($\rm b$) and for leptons ($\rm l=e$ or $\rm\mu$) can be presented, to the leading order in $\Omega/\Omega_{0}$, as

[TABLE]

Below we express all the perturbations through the displacements ${\boldsymbol{\xi}}_{\rm b}$ and ${\boldsymbol{z}}$, and the functions $\delta P$, $\Delta\mu_{\rm e}$, and $\Delta\mu_{\rm\mu}$. Thus, it is convenient to rewrite equation (29) with the help of equation (28) and the definition ${\boldsymbol{z}}={\boldsymbol{\xi}}_{\rm b}-{\boldsymbol{\xi}}$ as

[TABLE]

The Euler equation (13) and the superfluid equation (17) give us, to the leading order in $\Omega/\Omega_{0}$,

[TABLE]

These equations decouple into two systems, equations (28), (31) and equations (30), (32) (Andersson & Comer 2001; Lee & Yoshida 2003; Andersson et al. 2009). Equations (28) and (31) describe the normal r-modes, analogous to ordinary r-modes of non-superfluid NSs, while equations (30) and (32) describe superfluid modes driven by the relative motion (represented by the vector ${\boldsymbol{z}}$) of superfluid and normal (non-superfluid) liquid components. The solution to these two systems gives the following formulas for eigenfrequencies

[TABLE]

and eigenfunctions

[TABLE]

of normal and superfluid modes, respectively. In equations (35) and (36) $P_{l}^{m}({\rm cos}\theta)$ are the Legendre polynomials.

Since the function $h_{1}(x)/h(x)$ in equation (34) generally varies throughout the star, the frequency (34) cannot be a global frequency for the star as a whole — each layer has its own different eigenfrequency. This indicates that there are no purely toroidal superfluid modes in the system, an admixture of poloidal component is required. In other words, for superfluid modes $z_{r}$ should be non-zero at $\Omega\rightarrow 0$, making them (mixed) inertial modes rather than r-modes. The same conclusion has been arrived at in Andersson et al. (2009), although these authors obtained this result in the limit of vanishing stellar temperature ($T=0$).

However, if we assume that $Y_{np}=0$ then $h_{1}(x)=h(x)$ (see equations 18–20) and the superfluid equation (17) simplifies,

[TABLE]

Then the eigenfrequency (34) of superfluid modes,

[TABLE]

is independent of $x$ (Andersson & Comer 2001; Lee & Yoshida 2003; Andersson et al. 2009), becoming a global solution. Consequently, purely toroidal (to leading order in $\Omega/\Omega_{0}$) superfluid modes are possible in the limit $Y_{\rm np}=0$. Generally $|Y_{\rm np}|\ll Y_{\rm pp},Y_{\rm nn}$ (see section 4 for a numerical example); also, when protons are non-superconducting, one has $Y_{\rm np}=0$. Thus, in this paper we restrict ourselves to analysis of only the limiting case of vanishing entrainment, $Y_{\rm np}=0$, which may be not too bad approximation, see a corresponding discussion in section 8.

We see that the eigenfrequencies of all oscillation modes with $\sigma\propto\Omega$ and vanishing radial displacement coincide (to the lowest order in $\Omega/\Omega_{0}$) at $Y_{\rm np}=0$, while their eigenfunctions depend on the functions $C_{lm}(x)$ and $C_{zlm}(x)$ (see equations 35 and 36). To find these functions we need to proceed to the next order in $\Omega/\Omega_{0}$. Then the following relation is useful (it follows from equation 22),

[TABLE]

Using it along with expansions (24)–(27), we find from the continuity equations for, respectively, baryons and leptons:

(I) Baryons

[TABLE]

(II) Leptons (${\rm l=e}$, ${\rm\mu}$)

[TABLE]

where to calculate the thermodynamic derivatives in these equations it is assumed that the independent variables are $P$, $\Delta\mu_{\rm e}$ and $\Delta\mu_{\rm\mu}$; $g=-\nabla P_{0}/w_{0}$; and we used the fact that $\nabla n_{i0}=(\partial n_{i}/\partial P)\,\nabla P_{0}$ (we remind that the chemical potential imbalances vanish in equilibrium, $\Delta\mu_{\rm l}=0$, see, e.g., Haensel et al. 2007).

The Euler and superfluid equations written up to the second order give (see Appendix A for some details on the derivation of these equations):

(III) $r$-component of the Euler equation

[TABLE]

(IV) $\phi$-component of the Euler equation

[TABLE]

(V) $r$-component of the superfluid equation

[TABLE]

(VI) $\phi$-component of the superfluid equation

[TABLE]

(VII) $\theta$-component of the Euler equation

[TABLE]

(VIII) $\theta$-component of the superfluid equation

[TABLE]

To solve this system, let us introduce the functions $\xi_{{\rm b}\theta}^{1}$, $\xi_{{\rm b}\phi}^{1}$, $z_{\theta}^{1}$, and $z_{\phi}^{1}$ as a sum of toroidal and poloidal components (e.g., Saio 1982)

[TABLE]

As in the normal (nonsuperfluid) stars (e.g., Lockitch & Friedman 1999), we can now expand the functions $z_{r}^{1}(x,\theta)$, $\xi_{{\rm b}r}^{1}(x,\theta)$, $Q(x,\theta)$, $T(x,\theta)$, $Q_{z}(x,\theta)$, and $T_{z}(x,\theta)$, as well as $\delta P_{1}(x,\theta)$, $\Delta\mu_{\rm e1}(x,\theta)$ and $\Delta\mu_{\rm\mu 1}(x,\theta)$, into Legendre polynomial series with fixed $m$,

[TABLE]

where the summation goes over $l_{2}$ and $l_{1}$ ($l_{2}=m+2k+1$ and $l_{1}=m+2k$, $k=0$, $1$, $2$, $\ldots$) and no summation over $m$ is assumed.

Substituting now these expansions, as well as expressions (35) and (36) into equations (I)–(VIII) and equating coefficients at the same Legendre polynomials, we first of all see that only the solutions with $l=m$ result in a non-contradictory system of equations (to see this, it is sufficient to consider equations (III)–(VI) only). Then it can be noted that only the terms in expansions (52)–(55) and (58)–(60) with $k=0$ contribute to the terms proportional to the ‘lowest’ Legendre polynomial in each equation, while the functions $T^{1}(x,\theta)$ and $T^{1}_{z}(x,\theta)$ completely drop out from those terms. The system of equations resulting from the lowest Legendre polynomials is thus closed, and since in what follows we will be interested only in that system, we can skip the summation in (52)–(55), accounting only for the first terms there. Other terms in the expansion [as well as the functions $T(x,\theta)$ and $T_{z}(x,\theta)$] cannot be constrained from equations (I)–(VIII); to constrain them one needs to work in the next order in $\Omega/\Omega_{0}$.

Setting $l=m$, we obtain a system of four first-order differential equations for the functions $\xi_{{\rm b}r}^{1}(x)$, $z_{r}^{1}(x)$, $C(x)$, and $C_{z}(x)$ (from here on we skip the indices $lm$ and $l_{2}m$) of the form

[TABLE]

where the vector $d(x)\equiv\left[\xi_{{\rm b}r}^{1}(x),z_{r}^{1}(x),C(x),C_{z}(x)\right]$; $A(\sigma_{1},x)$ is a certain matrix depending on various thermodynamic functions, entrainment matrix coefficients, NS oblateness, pressure profile etc. With the appropriate boundary conditions (see section 3), the system (61) constitute the eigenvalue problem for $\sigma_{1}$. It is remarkable, that for $\rm npe$-matter (when muons are absent) the solution to the system (61), as well as the value of $\sigma_{1}$ can be found analytically (see Appendix B).

3 Boundary conditions

At the stellar surface [$x=R+O(\Omega^{2}/\Omega_{0}^{2}$)] we require that the Lagrangian perturbation of the pressure to be zero,

[TABLE]

At the stellar centre ($x\rightarrow 0$) the system (61) gives the following asymptotes for $\xi_{{\rm b}r}^{1}(x)$ and $C(x)$,

[TABLE]

Additionally, if neutrons in the centre are superfluid, one also has

[TABLE]

Analysis of equations (III)–(VI) requires the functions $C(x)$ and $C_{z}(x)$ to be continuous throughout, respectively, the star and the superfluid region (in particular, there should be no discontinuities of $C(x)$ and $C_{z}(x)$ at the crust-core interface and at the threshold density where muons appear) 333We do not account for the density jumps in the crust and at the core-crust interface.. Next, the continuity equation for baryons, equation (I), requires the function $\xi_{{\rm b}r}^{1}(x)$ also to be continuous throughout the star. The continuity equation for leptons, equation (II), implies the continuity of $\xi_{r}^{1}(x)$. Thus, $z_{r}^{1}(x)$ should vanish at the superfluid/non-superfluid interface [since $z_{r}^{1}(x)=\xi_{{\rm b}r}^{1}(x)-\xi_{r}^{1}(x)$ and $z_{r}^{1}(x)=0$ in the non-superfluid region]. These boundary conditions together with the system (61) constitute an eigenvalue problem for $\sigma_{1}$, which should be solved numerically.

4 Physics input

In our numerical calculations we adopt the Heiselberg & Hjorth-Jensen (1999) parametrization of APR equation of state (Akmal, Pandharipande, & Ravenhall 1998) for the matter in the core (at densities $\rho>\rho_{\rm cc}$) and the equation of state BSk20 (Potekhin, Fantina, Chamel, Pearson, & Goriely 2013) to describe the crust ($\rho<\rho_{\rm cc}$). We choose $\rho_{\rm cc}\approx 1.30\times 10^{14}\,\rm g\,cm^{-3}$ to avoid discontinuity in the density profile (at this density both EOSs give the same pressure). Notice, that the parametrization of Heiselberg & Hjorth-Jensen (1999) allows for muons ($\rm npe\mu$-composition), that appear first at $n_{\rm b\mu}\approx 0.133\,\rm fm^{-3}$ ($\rho_{\rm\mu}\approx 2.26\times 10^{14}\,\rm g\,cm^{-3}$). To illustrate the effect of muons, in this paper we consider two NS models: model I is simplified, muons are artificially switched off at any density; model II is more realistic, it allows for the presence of muons at $\rho>\rho_{\rm\mu}$. The figures presented in this section are plotted for the model II.

All numerical examples considered in what follows (except for Fig. 6) are obtained for an NS with the mass $M=1.4M_{\odot}$. For the model II (accounting for muons) the circumferential radius for such a star is $R=12.1$ km, the central density is $\rho_{\rm c}=9.48\times 10^{14}$ g cm*-3*; for the model I (muons are switched off) $R=12.2$ km, $\rho_{\rm c}=9.26\times 10^{14}$ g cm*-3*.

To find the oblateness of a rotating NS we used the Hartle scheme (Hartle 1967; Hartle & Thorne 1968). The function $\alpha(x)$ (see equation 22) that parametrizes the oblateness in oscillation equations is shown in Fig. 1.

Note that, to solve the Newtonian oscillation equations we, somewhat inconsistently, used hydrostatic NS models, calculated in full General Relativity (for both rotating and non-rotating NS configurations); the EOS employed by us is also relativistic in a sense that $P\sim\epsilon$ (pressure is comparable to the energy density). One should bear in mind this inconsistency which can affect the results. To make the calculations more reliable it would be desirable to describe also NS oscillations within the fully relativistic framework.

When modelling the effects of superfluidity we allow for the neutron and proton superfluidity in the NS core, while neutron pairing in the crust is neglected (the crustal superfluidity should not affect strongly the global oscillations of NSs). The adopted model of nucleon superfluidity is presented in Fig. 2, where the critical temperature profiles $T_{{\rm c}i}(x)$ for neutrons and protons are shown as functions of $x$. This model qualitatively does not contradict the results of microscopic calculations (see, e.g., Lombardo & Schulze 2001; Yakovlev, Levenfish, & Shibanov 1999; Gezerlis, Pethick, & Schwenk 2014; Dong, Lombardo, & Zuo 2014) and is analogous to the nucleon pairing models used to explain observations of cooling isolated NSs (Page, Lattimer, Prakash, & Steiner 2004; Gusakov, Kaminker, Yakovlev, & Gnedin 2004, 2005; Shternin, Yakovlev, Heinke, Ho, & Patnaude 2011; Page, Lattimer, Prakash, & Steiner 2013).

With decrease of the stellar temperature $T$ the size of SFL-region [the region where neutrons are superfluid, it is given by the condition $T<T_{\rm cn}(x)$] either increases or, at sufficiently low $T$, remains unchanged. For instance, SFL-region corresponding to $T=4\times 10^{8}$ K, is shaded in Fig. 2. One can see that the three-layer configurations are possible in some temperature range, with no neutron superfluidity in the centre and in the outer region, but with the superfluid intermediate region. At lower temperatures the star becomes a two-layer one, with the inner superfluid and outer normal region.

The symmetric entrainment matrix $Y_{ik}$ that parametrizes the effects of superfluidity in the oscillation equations, is calculated in a way similar to how it was done in Kantor & Gusakov (2011). Its elements are presented in Fig. 3 as functions of the normalized radial coordinate $x/R$ at fixed temperature $T=2\times 10^{8}\,\rm K$ (panel a), and as functions of $T$ at fixed $x/R=0.6$ (panel b). One can see that the non-diagonal element $Y_{\rm np}$ is much smaller than the diagonal ones ($Y_{\rm np}\sim 0.1Y_{\rm pp}$ and $Y_{\rm np}\sim 0.02Y_{\rm nn}$) and, in addition, $Y_{\rm np}$ vanishes when protons are non-superfluid. More realistic mean-field models, (see, e.g., Gusakov et al. 2009b, a, 2014c), give similar or slightly higher ratios (by a factor of $2-3$ in the cited papers). This motivates us to consider the limit $Y_{\rm np}=0$ in this paper. In this limit $h(x)=h_{1}(x)$ and superfluid equation (17) reduces to (37), see section 2. To further justify the assumption $Y_{\rm np}=0$, we plot the ratio $h/h_{1}$ in Fig. 4, which illustrates that $h(x)$ indeed approximately coincides with $h_{1}(x)$. Panels (a) and (b) of Fig. 4 show $h/h_{1}$ as a function of $x/R$ at $T=2\times 10^{8}\,\rm K$ and as a function of $T$ at $x/R=0.6$, respectively. To plot the figure, we used equations (18), (19) and (20) with non-zero $Y_{\rm np}$. The parameter $h$, entering the superfluid equation (37), is plotted in Fig. 5 as a function of $x/R$ at $T=2\times 10^{8}\,\rm K$ (panel a) and as a function of $T$ at fixed $x/R=0.6$ (panel b). To calculate it we used equations (18) and (20) with $Y_{\rm np}=0$.

5 Results for the spectrum

5.1 Model I, no muons

Below we shall present the numerical results only for the modes with $m=2$. These modes are the most interesting for us, since they can interact resonantly (see section 8) with the normal r-mode with $m=2$, which is known to be the most unstable inertial mode in NSs. The results of our numerical solution to the system of differential equations (61) with the boundary conditions from section 3 are presented in Fig. 6, where the function $\sigma_{1}(T)$ is shown by solid lines. It is calculated by employing the physics input from section 4 (we remind that $\sigma_{0}$ is the same for all the toroidal modes with $\sigma\propto\Omega$). We find that a superfluid NS with the simplest $\rm npe$-composition of the core harbours only two toroidal modes with $\sigma\propto\Omega$ (marked I and II in Fig. 6). At any given temperature (except near avoided crossings) these two modes, are of clearly distinct nature, either normal or superfluid. Normal mode is practically independent of temperature and almost coincides with the ordinary r-mode in a normal (non-superfluid) star. In contrast, superfluid mode depends on temperature via (i) temperature-dependent function $h$ and (ii) shrinking of the superfluid region with increasing $T$. Figure 7 presents the eigenfunctions for these two modes versus radial coordinate $x$ in the low-temperature limit. One can clearly see that the normal mode corresponds to co-motion, while the superfluid mode — to counter-motion of superfluid and normal liquid components. However, the situation is not so distinct in the vicinity of the avoided-crossings ($T\approx 3.8\times 10^{7}\,\rm K$ and $T\approx 1.5\times 10^{8}\,\rm K$ in our numerical example), where the modes change their behaviour from normal-like to superfluid-like and vice-versa. Near these ‘resonance’ temperatures the modes are ‘hybrids’ with the mixed properties (see Chugunov & Gusakov 2011; Gusakov et al. 2013; Gualtieri et al. 2014; Gusakov et al. 2014a, b for a detailed discussion).

To illustrate the sensitivity of the oscillation spectrum to the model of superfluidity and the stellar mass, we carried out the same calculations for the density-independent critical temperatures $T_{\rm cn}=6\times 10^{8}\,\rm K$ and $T_{\rm cp}=5\times 10^{9}\,\rm K$, assuming two different stellar masses: (i) $M=1.4M_{\odot}$ (the results are shown by dots in Fig. 6), and (ii) $M=1.7M_{\odot}$ (dashes in Fig. 6). One can see that the superfluid mode is sensitive to the critical temperature profiles, while the normal mode is not. At the same time, the variation of the stellar mass affects both modes.

5.2 Model II, muons are present

The numerical results for the spectrum of a NS containing muons in the core are shown in two panels of Fig. 8. The right panel is a strongly zoomed in version of the left one. We find that an admixture of muons, populates the NS spectrum of toroidal modes (in addition to two nodeless r-modes, existing in $\rm npe$ NSs) with an infinite set of superfluid modes, which differ from one another by the number of nodes $n$ of the superfluid function $C_{z}(x)$. We show only the first four superfluid modes with $n=0$, $1$, $2$, $3$ in the left panel of the figure. The right panel shows the fragment of the left one which contains only two nodeless r-modes.

The eigenfunctions for the first four superfluid modes (one nodeless mode and three modes with nodes) are plotted in Fig. 9 assuming $T=0$. Each panel is marked with the number of nodes of the superfluid function $C_{z}(x)$.

Although the muon presence results in additional superfluid toroidal modes with $n\neq 0$ in the spectrum, they do not affect much the eigenfrequencies of the modes without nodes. To illustrate this point, we show by dots the spectrum of an NS with the npe-composition of the core in the right panel of Fig. 8. Notably, the difference is stronger at low temperatures, at which the Brunt-V$\ddot{\rm a}$is$\ddot{\rm a}$l$\ddot{\rm a}$ frequency is higher and stratification is more pronounced (Kantor & Gusakov 2014). Figure 10 illustrates how the admixture of muons affects the eigenfunctions of normal r-mode. Solid lines and dashed lines show, respectively, the functions $C(x)$ and $C_{z}(x)$ for $\rm npe$ and $\rm npe\mu$ core compositions. $C(x)$ appears to be almost insensitive to the composition (two lines are practically indistinguishable in the figure), while the superfluid displacement $C_{z}(x)$ is substantially larger in the $\rm npe\mu$ case. Thus, one can expect that the core composition may affect the r-mode dissipation due to mutual friction, which is sensitive to the difference between the normal and superfluid velocities (Mendell 1991; Lindblom & Mendell 2000; Andersson et al. 2006), that is to the function $C_{z}(x)$ (see section 7).

6 Why the results for $\rm npe$ and $\rm npe\mu$ core compositions are so different?

Our numerical analysis reveal that NSs with $\rm npe$ core composition host only two nodeless r-modes (i.e., normal and superfluid r-modes), while account for muons leads, in addition, to an infinite set of superfluid r-modes with the nodes. What is the reason for such a difference? To answer this question, consider the continuity equations for leptons. In the coordinates $r$, $\theta$, $\phi$ they read:

[TABLE]

For inertial modes $\delta n_{\rm l}$ vanishes to the leading order in $\Omega/\Omega_{0}$ (i.e., $\delta n_{\rm l}=O(\Omega^{2}/\Omega_{0}^{2})$, see equation 24), while $n_{\rm l}$ can be considered as $\theta$-independent and factored out of the corresponding derivative,

[TABLE]

Combining equations for electrons and muons, we find

[TABLE]

Since, generally, $\partial(n_{\rm e}/n_{\rm\mu})/\partial r\neq 0$, $\xi_{r}$ has to vanish in the leading order in $\Omega/\Omega_{0}$, $\xi_{r}=O(\Omega^{2}/\Omega_{0}^{2})$, for inertial modes in the superfluid $\rm npe\mu$ matter. This conclusion, of course, does not concern the functions $\xi_{{\rm b}r}$ and $z_{r}$. In the superfluid $\rm npe$ matter such constraint is absent, because superfluid neutrons do not move with the same velocity as electrons, while protons are non-stratified ($n_{\rm p}=n_{\rm e}$).

To put it differently, superfluidity makes $\rm npe$ matter non-stratified (the corresponding Brunt-V$\ddot{\rm a}$is$\ddot{\rm a}$l$\ddot{\rm a}$ frequency vanishes if we neglect small entropy contribution, see Gusakov & Kantor 2013), allowing for non-zero $\xi_{r}$, while admixture of muons again stratifies it (Kantor & Gusakov 2014; Dommes & Gusakov 2016; Passamonti et al. 2016), requiring $\xi_{r}=0$ at $\Omega\rightarrow 0$. The situation is very similar to what happens in non-superfluid NSs: In a stratified star $\xi_{r}$ must vanish at $\Omega\rightarrow 0$ (Yoshida & Lee 2000), but this is not necessary in non-stratified stars.

As in non-superfluid NSs, the condition $\xi_{r}=0$ transforms superfluid inertial modes (mixed modes with non-zero $\xi_{r}$) in $\rm npe$-matter to purely toroidal r-modes in $\rm npe\mu$-matter, producing thus an infinite set of superfluid r-modes which differ by the number of nodes $n$ of the eigenfunction $C_{z}(x)$ (see section 8 for more details).

7 Dissipative damping of r-modes

Let us now illustrate how the avoided-crossings of r-modes affect their dissipative damping. At resonance temperatures (near avoided-crossings) the most effective dissipation mechanism is the mutual friction (Alpar et al. 1984; Lindblom & Mendell 2000; Lee & Yoshida 2003; Gusakov et al. 2014a, b). According to Alpar et al. (1984), it arises from the electron scattering off the neutron vortices magnetized due to entrainment between the superfluid neutrons and superconducting protons. Note that, this mechanism is relevant only when $Y_{\rm np}\neq 0$. Thus, although we set $Y_{\rm np}=0$ when calculating the r-mode eigenfunctions and eigenfrequencies (which should be a reasonable approximation for sufficiently small $Y_{\rm np}$, see a corresponding discussion in section 8), we, at the same time, assume that $Y_{\rm np}$ is finite when studying the r-mode damping.

The general relativistic expression for the oscillation energy dissipation rate due to mutual friction (per unit volume), $\dot{\epsilon}_{\rm MF}$, can be deduced from Gusakov & Dommes (2016) (see also Kantor & Gusakov 2012; Dommes & Gusakov 2016; Gusakov 2016). In the Newtonian limit, we are interested in here, the corresponding expression can be obtained by making use of the equation (90) from Gusakov & Dommes (2016), as well as the formulas (6) and (7) and the definition of ${\boldsymbol{z}}$. As a result, we arrive at the following general formula (in dimensional units),

[TABLE]

which reduces to (up to the leading order in the rotation frequency)

[TABLE]

in the case of r-modes. Here $\beta_{\rm MF}$ is the mutual friction coefficient from the formula (7) (the same coefficient is introduced in Khalatnikov 1989; Sedrakian & Sedrakian 1995; Kantor & Gusakov 2012; Gusakov 2016; Dommes & Gusakov 2016; Gusakov & Dommes 2016). In the NS literature the dimensionless coefficient $B$ is often used instead of $\beta_{\rm MF}$ (see, e.g., Mendell 1991; Andersson et al. 2006; Haskell et al. 2009; note that the coefficient $B$ in Mendell 1991; Andersson et al. 2006; Haskell et al. 2009 differs from $B$ in Hall & Vinen 1956!). In the zero temperature limit they are related by $B=\beta_{\rm MF}\,n_{\rm n}\mu_{\rm n}/c^{2}$ 444To obtain this relation, compare the equation on the neutron superfluid velocity (I7) from Gusakov & Dommes (2016) with the mutual friction force given by (I8) and, e.g., corresponding equations (2) and (3) in Haskell et al. (2009). Note that the vector $\bf{w}_{\rm np}$ in Haskell et al. (2009) equals the vector $\bf{W}$ in Gusakov & Dommes (2016).. Using this relation one can check that equation (70) reduces to the equation (114) of Haskell et al. (2009).

7.1 Mutual friction coefficient

To calculate the mutual friction coefficient $\beta_{\rm MF}$ one needs to know the time-scale $\tau_{\rm v}$ on which electrons (and other species strongly coupled to electrons, i.e. muons, protons, and normal neutron excitations; hereafter the normal liquid component) relax to the motion of the neutron vortices. For cold $\rm npe$-matter it was derived in Alpar et al. (1984) (their equation 30b). At a finite temperature and in $\rm npe\mu$ matter $\tau_{\rm v}$ has not been estimated yet. In that case not only electrons, but also muons and proton thermal excitations can scatter off the magnetized neutron vortices. Moreover, at finite temperatures both neutron and proton thermal excitations can, in principle, scatter off the neutron excitations localized in the vortex cores. These issues clearly deserve a separate study. Here, for definiteness, we shall use the result of Alpar et al. (1984) for cold $\rm npe$-matter, only slightly modifying it. Namely, we shall respect that electrons are coupled not only to protons but also to muons and normal neutron excitations. Then the relaxation time can be presented as (Alpar et al. 1984)

[TABLE]

where $\mu_{\rm n}n_{\rm b}-\mu_{\rm n}\mu_{k}Y_{{\rm n}k}$ represents the density of the normal liquid component (Gusakov & Andersson 2006), and the last equality should be considered as the definition of the function $g(\beta)$. In equation (72) $J_{1}$ is the Bessel function of the first kind; $\alpha=2k_{\rm e}\xi$, $\beta=\xi/\Lambda_{\star}$, $\tau_{0}^{-1}=\pi N_{\tau}\Phi_{\star}^{2}$, where $\xi$ is the neutron coherence length, $\Lambda_{\star}=\left(4\pi e^{2}Y_{\rm pp}\right)^{-1/2}$ is the London penetration depth, $k_{\rm e}$ is the electron Fermi wave vector, $\Phi_{\star}=\frac{\pi\hbar c}{e}\frac{Y_{\rm np}}{Y_{\rm pp}}$ is the flux of the neutron vortex line, and, finally,

[TABLE]

where $n_{\rm v}$ is the vortex density. It is possible to express (see appendix C) $\beta_{\rm MF}$ through $\tau_{\rm v}$ by the formula

[TABLE]

Both the coefficient $\beta_{\rm MF}$ (or $B$) and the relaxation time $\tau_{\rm v}$ depend on the entrainment matrix elements $Y_{ik}$, which are, generally, functions of the Landau parameters $f_{1}^{ik}$ and temperature (Gusakov & Haensel 2005; Gusakov et al. 2009b, a). We would also like to note that an approximate formula for $B$, adopted in the literature (see, e.g., equation 66 in Andersson et al. 2006), can hardly be a good approximation even in the zero-temperature limit, since it assumes $g(\beta)=0$. We plot $1-g(\beta)$ as a function of $x/R$ in Fig. 11 to demonstrate that the coefficient $1-g(\beta)$ in equation (72) is not close to 1. The curves correspond to two stellar temperatures ($T=10^{7}\,\rm K$ and $T=5\times 10^{8}\,\rm K$); the figure is plotted for constant $T_{\rm cn}=6\times 10^{8}\,\rm K$ and $T_{\rm cp}=5\times 10^{9}\,\rm K$.

7.2 Dissipation rate and energy of oscillations

The mutual friction damping time for r-modes can be defined as (Lindblom & Mendell 2000; Lee & Yoshida 2003; Haskell et al. 2009; Andersson et al. 2009)

[TABLE]

where $E$ is the mechanical energy of oscillations coinciding, to the leading order in $\Omega/\Omega_{0}$, with their kinetic energy (e.g., Lindblom & Mendell 2000; Haskell et al. 2009). The latter is given by (Gusakov et al. 2013)

[TABLE]

Substituting the real parts of $\xi_{\rm b}$ and $z$ from (35) and (36), respectively, we get

[TABLE]

Integration of this formula over $\phi$ makes the energy independent of the phase of oscillations and introduces an additional factor $\pi$, while integration over $\theta$ results in a factor $2m(m+1)(2m)!/(2m+1)$, so that

[TABLE]

Next, $\dot{E}_{\rm MF}$ in equation (75) is the rate of oscillation energy damping due to mutual friction,

[TABLE]

(Actually, the integrand is non-zero only in the region where both neutrons and protons are superfluid.) Substituting (71) with $z_{\theta}$ and $z_{\phi}$ from (36), we obtain

[TABLE]

Integration over $\theta$ and $\phi$ gives

[TABLE]

which is independent of the phase of oscillations. Then

[TABLE]

where the last equality is the definition of $\tau_{0\,\rm MF}$ 555Note that both $C(x)$ and $C_{z}(x)$ excite in the leading order in $\Omega$, thus the ratio of integrals in equation (82) is independent of the rotation frequency..

7.3 Results for damping times and instability window

Figure 12 shows the inverse damping time, $\tau_{0\,\rm MF}$, (see equation 82) as a function of temperature $T$ for two nodeless r-modes (I and II, see Figs. 6 and 8). $C_{z}(x)$ is larger for the superfluid mode, hence it dissipates stronger, as it should be (see equation 82). Correspondingly, when a mode changes its behaviour from normal-like to superfluid-like, the efficiency of mutual friction dissipation for that mode rises sharply. The rapid decrease of mutual friction efficiency at high temperatures (at T\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;2\times 10^{8} K in the figure) is due to decrease of the entrainment matrix element $Y_{\rm np}$ at $T\rightarrow T_{\rm cn}$. Normal mode in $\rm npe\mu$ matter (solid lines) dissipates faster than in $\rm npe$ matter (dots) because $C_{z}(x)$ is higher in $\rm npe\mu$ matter (see Fig. 10).

Having determined the mutual friction damping time $\tau_{\rm MF}$, we can calculate the instability windows, that is the regions of stellar frequency and temperature at which a star becomes unstable with respect to excitation of r-modes. They are determined by the condition $1/\tau_{\rm GR}+1/\tau_{\rm Diss}<0$, where $\tau_{\rm GR}$ is the gravitational radiation timescale (see, e.g., Andersson & Kokkotas 2001; Gusakov et al. 2014a) and $\tau_{\rm Diss}$ is the dissipation timescale. As discussed, e.g., in Gusakov et al. (2014a), at the temperatures of interest the main dissipative mechanisms are the mutual friction dissipation and shear viscosity. Consequently, $1/\tau_{\rm Diss}=1/\tau_{\rm MF}+1/\tau_{\rm S}$, where $\tau_{\rm MF}$ is given by equation (82) and $\tau_{\rm S}$ (the dissipation timescale due to shear viscosity) is taken from Gusakov et al. (2014a), see their formula (13) 666Formula (13) of Gusakov et al. (2014a) assumes a polytropic NS with strong proton superconductivity. Here we apply the same formula to our NS models with critical temperature profiles from Fig. 2..

The instability windows, calculated as described above, are presented in Fig. 13 for two nodeless r-modes (I and II, see solid lines). In the region filled with grey the $\rm npe\mu$ NSs are stable (this is the ‘stability’ region), in the white region they are unstable. Circles with error bars show the observational data (taken from Gusakov et al. 2014a, 2016). Dots indicate the instability window for $\rm npe$ NSs.

Generally, the stability region, calculated not accounting for the resonant interaction of oscillation modes (but under realistic assumptions about kinetic properties of NS matter) is rather small and cannot contain the hottest and most rapidly rotating NSs (see, e.g., Ho et al. 2011; Haskell et al. 2012; Gusakov et al. 2014a; Chugunov et al. 2014; Gusakov et al. 2016 and the dashed curve in Fig. 13). Obviously, this fact poses a problem for the theory (Levin 1999; Ho et al. 2011; Haskell et al. 2012; Gusakov et al. 2014a, b). On the other hand, the realistic stability region easily accommodates all the observed sources, but it is obtained under a simplifying assumption $Y_{\rm np}=0$. What would happen if we relaxed this assumption is discussed in the next section.

8 Discussion

Prior to discussing the spectrum of SFL NSs let us first analyze the non-superfluid case. For that purpose the work of Yoshida & Lee (2000) can be very helpful (see also Lockitch & Friedman 1999). At high rotation frequencies the spectra of stratified and non-stratified NSs are similar — there is one nodeless r-mode and an infinite set of inertial modes. At low rotation frequencies stratification transforms inertial modes (the modes with mixed poloidal and toroidal components) into either g-modes or r-modes with nodes and $l=m$. In accordance with the connection rule (21)–(22) of Yoshida & Lee (2000), inertial modes $i_{k}$ with $k\geq 1$ transform into g-modes, and $i_{0}$ modes transform into r-modes 777Yoshida & Lee (2000) expand the equations describing inertial modes into series in terms of the Legendre polynomials, $P_{l}^{m}({\rm cos}\theta)$. For each $m$ there exists a system of independent equations that couples Legendre polynomials with different $l$. As a result, any inertial mode is a mixture of various Legendre polynomials $P_{l}^{m}({\rm cos}\theta)$ with fixed $m$. It turns out that for any such mode one can identify the leading toroidal harmonic [$P_{l_{0}-1}^{m}({\rm cos}\theta)$ in notations of Yoshida & Lee 2000], so that it can be characterized by two angular ‘quantum’ numbers, $l_{0}$ and $m$. However, after fixing $l_{0}$ and $m$ we still have some freedom resulting in a set of inertial modes with the frequencies $\sigma_{k}$. Yoshida & Lee (2000) enumerate these modes following the rule (for axial modes): $0\leq\sigma_{0}<\sigma_{1}<...<\sigma_{k}<...<\sigma_{(l_{0}-1-|m|)/2}$. Thus, generally, an inertial mode can be characterized by three numbers, $l_{0}$, $m$, and $k$. . This transformation occurs at rotation frequencies of the order of g-mode frequencies in a non-rotating NS. Thus, at slow rotation the toroidal modes with $\sigma\sim\Omega$ and $l=m$ of a stratified star are: (i) one nodeless r-mode and (ii) an infinite set of r-modes with nodes. In non-stratified stars there is only one nodeless r-mode.

The spectrum of SFL NSs is rather similar. In superfluid npe-matter the Brunt-V$\ddot{\rm a}$is$\ddot{\rm a}$l$\ddot{\rm a}$ frequency vanishes, making NSs with $\rm npe$ core composition (hereafter $\rm npe$ NSs) non-stratified (Gusakov & Kantor 2013), while even small admixture of muons stratifies them again (Kantor & Gusakov 2014; Dommes & Gusakov 2016; Passamonti et al. 2016). Our calculations show that, at low rotation rates (and at $Y_{\rm np}=0$) the only toroidal modes with $\sigma\sim\Omega$ in $\rm npe$ NSs are two nodeless r-modes (which behave as normal and superfluid modes at a given temperature), while $\rm npe\mu$ NSs harbours, in addition, an infinite set of $r_{n}$-modes with nodes (here $n\geq 1$ is the number of nodes of $r$-mode, see section 5.2). These modes are of purely superfluid nature (counter-motion dominates, see section 5.2 and Fig. 9), and hence we will denote them $r^{s}_{n}$-modes. It is quite natural to assume (although we have not checked this by direct calculation) that these $r^{s}_{n}$-modes correspond to the superfluid inertial $i^{s}_{0}$-modes of $\rm npe$ NSs (in analogy with non-superfluid NSs discussed above).

So, why do the superfluid (mixed) inertial $i^{s}_{0}$-modes of $\rm npe$ NSs mutate into the superfluid $r^{s}_{n}$-modes under stratification by muons? Let us remind that for the modes with $\sigma\sim\Omega$ the electron and muon continuity equations imply that $\xi_{r}\rightarrow 0$ at $\Omega\rightarrow 0$ (see section 6), which means that $\xi_{{\rm b}r}=z_{r}$ at $\Omega\rightarrow 0$. Note also that there is a clear distinction between the normal and superfluid modes due to a very good decoupling of superfluid and normal degrees of freedom (see Haskell et al. 2009; Gusakov & Kantor 2011; Kantor & Gusakov 2012; Dommes & Gusakov 2016). In the limiting case of perfect decoupling the superfluid modes cannot excite baryon displacements $\xi_{\rm b}$, and thus can meet the condition $\xi_{r}=0$ only by vanishing $z_{r}$ (i.e., transforming into purely toroidal modes). Our numerical results indicate that in the situation of approximate decoupling this conclusion remains the same.

Moreover, our results imply that stratification by muons does not transform normal (mixed) inertial modes $i_{0}$ into normal r-modes with the nodes (we do not find such modes in the calculated spectrum). We expect that they survive in $\rm npe\mu$ NSs. This assumption seems to be reasonable, since stratification by muons put only one constraint, $\xi_{r}=0$ at $\Omega\rightarrow 0$, limiting the variety of (mixed) inertial modes, but not excluding them completely. Why do the normal inertial modes can retain their nature (in contrast to superfluid modes) under stratification by muons? Let us again illustrate our point by assuming perfect decoupling of superfluid and normal degrees of freedom. Perfect decoupling does not mean that normal modes cannot excite the superfluid displacement ${\boldsymbol{z}}$ (see Gusakov & Kantor 2011; Kantor & Gusakov 2012; Gusakov et al. 2013; Dommes & Gusakov 2016). Thus, the condition $\xi_{r}=0$ can be easily satisfied for normal inertial modes by choosing $z_{r}$ such that $\xi_{r}=\xi_{{\rm b}r}-z_{r}=0$ at $\Omega\rightarrow 0$.

By analogy with the non-superfluid star, it seems reasonable that, at low rotation rates (when the spin frequency is smaller than the corresponding g-mode frequency of a non-rotating NS, that is hundreds of Hz for superfluid $\rm npe\mu$ NSs, see Kantor & Gusakov 2014; Passamonti et al. 2016), stratification by muons transforms $i^{s}_{k}$-modes of $\rm npe$ NSs (with $k\geq 1$, see footnote 7 for the definition of $k$) into g-modes (they have finite $\sigma$ at $\Omega\rightarrow 0$ and hence are not required to have $\xi_{r}=0$ at $\Omega\rightarrow 0$). At the same time, normal $i_{k}$ ($k\geq 1$) modes do not transform into g-modes, retaining their nature. The scheme, illustrating this complex behavior, is presented in Fig. 14.

In accordance with the numerical results of Lee & Yoshida (2003), where SFL $\rm npe$ NSs were considered, it is $i^{s}_{0}$-modes (see footnote 7) that experience avoided-crossings with normal r-mode in $\rm npe$ NSs. As has already been discussed above, stratification by muons transforms $i^{s}_{0}$-modes into $r_{n}^{s}$-modes, studied here. In contrast, superfluid $i^{s}_{k}$-modes ($k\geq 1$) have rather high eigenfrequencies and do not meet the normal r-mode (Lee & Yoshida 2003) in $\rm npe$ NSs. Stratification by muons mixes $i^{s}_{k}$ modes ($k\geq 1$) with g-modes and thus (additionally) increases their eigenfrequencies to the corresponding g-mode eigenfrequencies in the limit $\Omega\rightarrow 0$. Consequently, $i^{s}_{k}$ modes ($k\geq 1$) should not interact (form avoided crossings) with normal r-mode in $\rm npe\mu$ NSs as well.

To sum up, the $r^{s}_{n}$-modes ($n\geq 1$) of $\rm npe\mu$ NSs calculated in this paper, being the counterparts of $i^{s}_{0}$-modes of $\rm npe$ NSs calculated in Lee & Yoshida (2003), were exactly the modes interesting for us since we were looking for the resonance interaction (avoided crossings) of normal r-mode with superfluid modes. Note that Lee & Yoshida (2003) does not discuss explicitly the interaction of superfluid nodeless r-mode with the normal r-mode. However, their results imply that it takes place in NSs with weak entrainment (when $Y_{\rm np}$ is small; see also section 5 and the work by Gusakov et al. 2014b).

In this paper we assumed that $Y_{\rm np}=0$. This limit implies $h(x)=h_{1}(x)$ (see equations 18–20) and it makes possible purely toroidal superfluid modes. In this simplified model we found that all purely toroidal modes have the same eigenfrequency at the leading order in rotation, $\sigma_{0}=\frac{2m}{l(l+1)}$, and only modes with $l=m$ exist. Next order approximation removes this degeneracy (see Figs. 6 and 8). As a result, normal r-mode exhibits avoided-crossings with the nodeless superfluid mode, while other $r^{s}_{n}$ modes (with $n\geq 1$) have substantially lower $\sigma_{1}$ and do not interact with the normal r-mode.

If one relaxes the assumption $Y_{\rm np}=0$ 888More rigorously, if one relaxes the assumption that $h_{1}(x)/h(x)$ is independent of $x$., one does not find global eigenfrequencies for purely toroidal superfluid modes in the leading order in $\Omega/\Omega_{0}$ (see the corresponding discussion in section 2 and in Andersson et al. 2009). What one finds, instead, is $x$-dependent eigenfrequency which should be regularized to some global value by admixing poloidal component to the solution, that is, by making superfluid radial displacement $z_{r}$ non-zero at $\Omega\rightarrow 0$. Still, the NS eigenmodes should vary continuously with variation of $Y_{\rm np}$. Thus, the limit $Y_{\rm np}=0$ should be a good approximation for small but non-zero $Y_{\rm np}$. How small should it be? The answer is not so straightforward.

At first glance, the real eigenfrequency should be of the order of the typical value of $\frac{2m}{l(l+1)}\frac{h_{1}(x)}{h(x)}$. While the variation of $\frac{2m}{l(l+1)}\frac{h_{1}(x)}{h(x)}$ (let us denote it $\delta\sigma_{0}$) throughout the star is small, the required admixture of the poloidal component should also be small, not affecting oscillations substantially. It seems that if $\sigma_{1}\Omega^{2}>\delta\sigma_{0}/\sigma_{0}$, the effect on eigenmodes due to second order rotational corrections should be more pronounced than the effect of non-zero $Y_{\rm np}$.

However, the situation may be not so simple. Let us remind that, due to the condition $\xi_{r}\rightarrow 0$ at $\Omega\rightarrow 0$ (valid for $\rm npe\mu$-matter, see section 6), $\xi_{{\rm b}r}=z_{r}$ at $\Omega\rightarrow 0$. Thus, the non-zero value of $Y_{\rm np}$ leads to the excitation of $\xi_{{\rm b}r}$ in superfluid modes. On the other hand, superfluid and normal modes are decoupled to a good accuracy (Gusakov & Kantor 2011; Kantor & Gusakov 2012; Dommes & Gusakov 2016). In the limiting case of perfect decoupling this means that superfluid modes cannot excite baryon displacements ${\boldsymbol{\xi}}_{\rm b}$. Thus, in this limiting case superfluid modes with $\sigma\sim\Omega$ should be completely pushed out of that part of $\rm npe\mu$ region where $Y_{\rm np}\neq 0$ (where protons are superconducting). In reality, the decoupling is not perfect and ${\boldsymbol{\xi}}_{\rm b}$ can be excited, to some extent, by superfluid modes in $\rm npe\mu$ region. However, the above consideration clearly indicates that decoupling leads to at least partial expulsion of the eigenfunctions from the $\rm npe\mu$ region at non-zero $Y_{\rm np}$. Consequently, even small value of $Y_{\rm np}$ can, in principle, substantially modify the superfluid eigenmodes if decoupling is good. All this complicated behaviour is driven by the interplay of three small parameters: rotation frequency, $Y_{\rm np}$, and decoupling parameter (see Gusakov & Kantor 2011 for the definition and more details on this parameter).

The above discussion clearly indicates that stratification by muons plays a key role in defining the low-frequency spectrum of superfluid NSs, whose rotation frequencies are smaller or comparable to the g-mode frequencies. Keeping in mind that typical g-mode frequencies in superfluid $\rm npe\mu$ NSs can reach hundreds of Hz (see Kantor & Gusakov 2014; Dommes & Gusakov 2016; Passamonti et al. 2016), account for muons appears to be crucial even for the most rapidly rotating NSs.

Summarizing, using the approximation of vanishing entrainment between the superfluid neutrons and superconducting protons, we calculated the spectrum of toroidal modes of slowly rotating superfluid neutron stars. We analysed NSs with both $\rm npe$ and $\rm npe\mu$ compositions of their core. We found that $\rm npe$ NSs harbour only two types of toroidal modes — nodeless normal and nodeless superfluid r-modes. In contrast, $\rm npe\mu$ NSs additionally possess an infinite set of superfluid r-modes with nodes. For both NS compositions, normal r-mode exhibits avoided-crossings with the nodeless superfluid r-mode. Other superfluid r-modes in $\rm npe\mu$ NSs have lower frequencies and do not approach (and form avoided crossings with) the normal r-mode. In addition, we calculated the ‘instability windows’ for r-modes, accounting for the mode dissipation due to mutual friction and shear viscosity. We found that the instability window is split up by two stability peaks so that all the observed sources can be interpreted as either lying inside the stability region or climbing the left edge of one of those peaks. This result confirms the phenomenological scenario of Gusakov et al. (2014a, b); Chugunov et al. (2014); Gusakov et al. (2016) which was suggested to explain stability of most rapidly rotating NSs. In the subsequent work we plan to relax the simplifying assumption $Y_{\rm np}=0$, adopted in the present paper, and to accurately calculate the temperature-dependent spectrum of inertial modes in rotating superfluid NSs.

Acknowledgments

We are grateful to Andrey Chugunov and Vasiliy Dommes for valuable comments and discussions.

Appendix A Derivation of equations of section 2.2

We start with $r$ component of Euler equation (13)

[TABLE]

Making use of

[TABLE]

that follows from the formula $\delta\epsilon=\mu_{\rm n}\delta n_{\rm b}$, and the fact that in the equilibrium $\Delta\mu_{\rm l}=0$, (83) can be rewritten as

[TABLE]

$\phi$ component of equation (13) gives us $\delta P$,

[TABLE]

where we introduced the vector ${\boldsymbol{a}}\equiv\sigma^{2}{\boldsymbol{\xi}}_{\rm b}-2\imath\sigma{\boldsymbol{\Omega}}\times{\boldsymbol{\xi}}_{\rm b}$.

$\theta$-component of the Euler equation is

[TABLE]

We will be interested in the first two terms in the expansion of the right-hand side of this equation in small parameter $\Omega/\Omega_{0}$. Thus, the independent variables $r$, $\theta$, $\phi$ should be changed to $x$, $\theta$, $\phi$ more cautiously. Applying equation (39) to equation (87), we obtain

[TABLE]

Making use of equation (83) and substituting the expression (86) into (88), we get (note that here $x$, $\theta$, and $\phi$ are independent variables)

[TABLE]

Substituting now the expression (22) for $r$, we arrive at

[TABLE]

This equation, written to the leading order in $\Omega/\Omega_{0}$, gives us equation (31); the next-to-leading order in $\Omega/\Omega_{0}$ gives equation (VII). Proceeding in a similar way with equation (37), one finds equations (32) (with $h=h_{1}$) and (VIII).

Appendix B npe-matter: Analytic solution

Here we present the oscillation equations and discuss their solution for NSs, whose cores are composed of npe-matter. To the leading order in $\Omega/\Omega_{0}$ the oscillation equations (equations 28, 30, 31, and 32) are the same as for $\rm npe\mu$-matter. However, in the next order the system of equations for npe-matter will differ from $\rm npe\mu$-equations (I)–(VIII) and consist of:

(I) the continuity equation for baryons:

[TABLE]

where the independent variables in thermodynamic derivatives are $P$ and $\Delta\mu_{\rm e}=\mu_{\rm n}-\mu_{\rm p}-\mu_{\rm e}$.

(II) the continuity equation for electrons:

[TABLE]

(III) $r$-component of the Euler equation

[TABLE]

(IV) $\phi$-component of the Euler equation

[TABLE]

(V) $r$-component of the superfluid equation

[TABLE]

(VI) $\phi$-component of the superfluid equation

[TABLE]

These equations should be supplemented with $\theta$-components of the Euler and superfluid equations, which coincide with the corresponding equations in $\rm npe\mu$-matter [see equations (VII) and (VIII) in the main text].

The analysis of this system of equations ensures that, as in the $\rm npe\mu$-case, only modes with $l=m$ are possible in $\rm npe$-matter (we remind that this system is valid for toroidal modes only). To solve this system, let us express $\Delta\mu_{{\rm e}1}$ (defined by equation 24) from equation (96) and substitute the answer into equation (95). The solution to the resulting equation is remarkably simple,

[TABLE]

where we set to unity the dimensional integration constant. Expressing now $\delta P_{1}$ from equation (94) and substituting $\delta P_{1}$ and $\Delta\mu_{{\rm e}1}$ into equation (93), we obtain the first-order inhomogeneous differential equation for $C(x)$,

[TABLE]

whose solution is

[TABLE]

where $C_{0}$ is an integration constant to be determined below.

Substituting now the expansions (52)–(57) into (46) divided by ${\rm sin}\theta$, and equating coefficients at the terms proportional to $P_{m}^{m}$, we obtain an algebraic equation relating $Q^{1}(x)$ to $\xi_{{\rm b}r}^{1}(x)$ and $C(x)$ (note that the function $T^{1}(x)$ cancels out in the terms proportional to $P_{m}^{m}$). Similarly, using equation (47), we relate $Q_{z}^{1}(x)$ to $z_{r}^{1}(x)$ and $C_{z}(x)$.

Now we substitute the derived expression for $Q^{1}(x)$, as well as (98), into the continuity equation for baryons (91). Equating then the coefficients at $P_{m+1}^{m}$ and using the expressions (99) and (97), we obtain the following first-order inhomogeneous differential equation for $\xi_{{\rm b}r}^{1}(x)$,

[TABLE]

where $A(x)$, $B_{00}(x)$, $B_{10}(x)$, $B_{01}(x)$, and $B_{11}(x)$ are certain functions of $x$ defined by equation (91). The solution to this equation is

[TABLE]

where $\xi_{0}$ is the integration constant and

[TABLE]

Following the same strategy with the electron continuity equation (92), we find for $z_{r}^{1}(x)$,

[TABLE]

where $x_{\rm sfl1}$ corresponds to the superfluid/nonsuperfluid inner boundary, $z_{0}$ is the integration constant and

[TABLE]

Analysis of the hydrodynamic equations in the stellar centre implies that $\xi_{{\rm b}r}^{1}(x)\propto x^{m}$, $z_{r}^{1}(x)\propto x^{m}$ at $x\rightarrow 0$, and thus $\xi_{0}=0$ and $z_{0}=0$. If the stellar centre is non-superfluid then $z_{0}=0$ anyway, because at the superfluid/non-superfluid interfaces $z_{r}^{1}=0$ (see section 3), and application of this condition to the inner boundary of the superfluid region ($x=x_{\rm sfl1}$) gives us $z_{0}=0$. The finiteness of the function $\xi_{{\rm b}r}^{1}(x)$ at the surface ($x=R$) and vanishing of $z_{r}^{1}(x)$ at the outer superfluid/non-superfluid boundary ($x={x_{\rm sfl2}}$) imply

[TABLE]

The system of equations (105)–(106) for $\sigma_{1}$ and $C_{0}$ has two solutions corresponding to two oscillation modes discussed in the text.

Appendix C Relating mutual friction coefficient $\beta_{\rm MF}$ to the relaxation time $\tau_{\rm v}$ at finite temperature

To relate the mutual friction coefficient $\beta_{\rm MF}$ (or $B$) to the microphysics parameters, it is customary to use (e.g., Mendell 1991; Andersson et al. 2006) the notion of the relaxation time $\tau_{\rm v}$ (Feibelman 1971; Sauls et al. 1982; Alpar et al. 1984). This time characterizes relaxation of the electron velocity, $\boldsymbol{u}$ (which is also the velocity of other particle species, coupled to electrons on time-scales much shorter than $\tau_{\rm v}$, i.e., muons, protons, and neutron thermal excitations; hereafter normal liquid component), to the velocity of the static (fixed) neutron vortices. Note that only the component of $\boldsymbol{u}$ orthogonal to the direction of vortices experiences relaxation; let us denote it $\boldsymbol{u}_{\perp}$. In the reference frame of fixed vortices (where their velocity ${\boldsymbol{V}}_{\rm L}=0$) the dissipative force acting on the flow of normal component is

[TABLE]

Here $(\mu_{\rm n}n_{\rm b}-\mu_{\rm n}\mu_{k}Y_{{\rm n}k})/c^{2}$ represents the density of normal component (Gusakov & Andersson 2006). In an arbitrary frame, if we relax the assumption about static vortices, the friction force on the normal component is

[TABLE]

The same dissipative force, taken with the opposite sign, acts on neutron vortices. Since the total force on vortices should vanish (they are practically massless), this dissipative force must be compensated by the Magnus force acting on vortices from the neutron superfluid. Consequently, we can conclude that vortices act on the neutron superfluid with the (minus) Magnus force, or with

[TABLE]

Consider now the weak-drag limit relevant for NS matter (Mendell 1991; Andersson et al. 2006). In this limit the vortex velocity is (see equation 101 in Gusakov & Dommes 2016, where we skipped the term $\propto\beta$, since it is small)

[TABLE]

and $\alpha=-1/(\mu_{\rm n}\mu_{k}Y_{{\rm n}k})$. In the locally flat comoving reference frame [where $u^{\mu}=(1,0,0,0)$] the spatial components of $V_{\rm(L)}^{\mu}$ are

[TABLE]

The dissipative force on the neutron superfluid is then

[TABLE]

On the other hand, this force can be deduced from the equation (7). Using equation (I7) of Gusakov & Dommes (2016), one can show that the force on the neutron superfluid equals $-\boldsymbol{F}=-n_{\rm n}\boldsymbol{f}\mu_{\rm n}\mu_{k}Y_{{\rm n}k}/c^{2}$, where ${\boldsymbol{f}}$ is the spatial part of the four-vector $f^{\mu}$ (see equation 100 of Gusakov & Dommes 2016), taken in the comoving frame. Using equation (100) from Gusakov & Dommes (2016), the latter expression can be represented, in the comoving frame, as

[TABLE]

where $\mathcal{V}_{({\rm M})}=2\Omega\mu_{\rm n}/c^{2}$ in the case of solid-body rotation. Comparing now (112) and (113), we obtain (in dimensional form)

[TABLE]

Appendix D Typos in the paper by Provost et al. (1981)

When considering the r-mode problem in superfluid NSs we followed the approach of Provost et al. (1981). Although generally this paper is extremely illuminating, it contains a number of typos and some incorrect statements in its Appendix. Here we would like to list them. First of all, the definition of $g$ [between the formulas (7d) and (8)] reads $g=\rho^{-1}\left|\frac{dp}{dx}\right|$. Next, the first equation in (10) should read

[TABLE]

Let us now turn to the Appendix of Provost et al. (1981), discussing barotropic stars. We note that:

(i) the corrections $\xi^{1}_{\theta}$ and $\xi^{1}_{\phi}$ do not have to be purely poloidal; however, the toroidal component does not appear in the lowest harmonics in equations (7c) and (7d).

(ii) In the second term in equation (A4) $m+3$ should be replaced with $(m+3)/x$.

(iii) In the solution (A5) $1/\gamma$ in the exponents should be replaced with $-1/\gamma$.

(iv) Finally, we believe that the discussion of the behaviour of the solution near the surface is not correct. The integral entering (A5) with $x^{-m-3}p^{-1/\gamma}(x)$ factored out (let us denote it by $\int^{x}_{0}\ldots$) is regular near the stellar surface (since $1/\gamma>1$), while the factor $1/p^{1/\gamma}(x)$ diverges. The only way to have a regular solution at the surface ($x=R$) is to require $\int_{0}^{R}\ldots=0$. This condition determines the eigenvalues $\sigma_{1}$.999Similar approach is used in appendix B to calculate the spectrum of superfluid $\rm npe$ NSs.

We checked this interpretation and found that it is correct. Namely, we calculated the values of $\sigma_{1}$ in the Cowling approximation in the case of a polytropic star (polytropic index $n=1$). We verified that our results are in agreement with those listed in the last column of Table I in Andersson & Kokkotas (2001). This table presents numerical values for the parameter $\omega_{2}$, defined as

[TABLE]

and related to $\sigma_{1}$ by the expression: $\sigma_{1}=\omega_{2}/\sigma_{0}$. Note that, actually, Andersson & Kokkotas (2001) give different definition for $\omega_{2}$ (see their equation 32). However, comparing the first and third columns of Table I in Andersson & Kokkotas (2001) with Table I in Lindblom et al. (1999) [in which the parameter $\kappa_{2}$ is presented, defined by equations (3.7) and (3.13)], one can check that, to generate Table I, Andersson & Kokkotas (2001) employed the definition (116) for $\omega_{2}$.