Resonance suppression of the r-mode instability in superfluid neutron stars: Accounting for muons and entrainment
aa r X i v : . [ a s t r o - ph . H E ] F e b Resonance suppression of the r-mode instability in superfluid neutron stars:Accounting for muons and entrainment
Elena M. Kantor, Mikhail E. Gusakov, Vasiliy A. Dommes
Ioffe Institute, Politekhnicheskaya 26, 194021 Saint-Petersburg, Russia (Dated:)We calculate the finite-temperature r-mode spectrum of a superfluid neutron star accounting forboth muons in the core and the entrainment between neutrons and protons. We show that thestandard perturbation scheme, considering the rotation rate as an expansion parameter, breaksdown in this case. We develop an original perturbation scheme which circumvents this problem bytreating both the perturbations due to rotation and (weak) entrainment simultaneously. Applyingthis scheme, we propose a simple method for calculating the superfluid r-mode eigenfrequency in thelimit of vanishing rotation rate. We also calculate the r-mode spectrum at finite rotation rate forrealistic microphysics input (adopting, however, the Newtonian framework and Cowling approxima-tion when considering perturbed oscillation equations) and show that the normal r-mode exhibitsresonances with superfluid r-modes at certain values of temperatures and rotation frequencies inthe parameter range relevant to neutron stars in low-mass X-ray binaries (LMXBs). This turns therecently suggested phenomenological model of resonance r-mode stabilization into a quantitativetheory, capable of explaining observations. A strong dependence of resonance rotation rates andtemperatures on the neutron superfluidity model allows us to constrain the latter by confrontingour calculations with the observations of neutron stars in LMXBs.
PACS numbers:
I. INTRODUCTION
As it was shown in 1998 [1, 2], in the absence of dissipation r-modes (predominantly toroidal oscillations of rotatingstars restored by the Coriolis force [3]) are unstable with respect to radiation of gravitational waves at any rotationfrequency of a star. In practice, r-mode instability is mostly interesting for neutron stars (NSs), since only for NSs dor-modes have a reasonably fast growth rate. Being excited, r-modes emit gravitational waves, which carry off angularmomentum from the star. Gravitational radiation back-reaction excites the r-mode by increasing its amplitude.Dissipation opposes this process. Calculations show that in cold NSs r-mode instability is effectively damped bythe shear viscosity, while in hot NSs it is damped by the bulk viscosity. In slowly rotating NSs r-mode instabilityis also effectively suppressed, since gravitational radiation back-reaction is weak at slow rotation. As a result, onlywarm rapidly rotating NSs may fall into the so called “instability window” (i.e., the region of stellar temperaturesand rotation rates, where r-modes are unstable). Such NSs are observed in the low-mass X-ray binaries (LMXBs).Numerous observations of these NSs pose a challenge, because modeling shows that NSs should quickly leave theinstability window in the course of their evolution in LMXBs [4], since excited r-modes heat up and spin down thestars rapidly. Various proposals for reconciling theory with observations have been discussed in the literature (see, e.g.the reviews [5, 6]). Many of them involve some exotic physics (e.g. the presence of hyperons/quarks in the NS cores),or make some model-dependent assumptions about the mechanism of nonlinear saturation of r-modes [5, 7–10].Here we shall focus on the r-mode stabilization mechanism proposed in Refs. [11, 12], which appeals to resonancestabilization of r-modes by superfluid (hereafter SF) modes. Interestingly, this mechanism involves minimal assump-tions about the properties of NS matter, such as minimal core composition (neutrons, protons and leptons) and SFof baryons. Neutron SF in the NS core gives rise to two independent velocity fields: the velocity of SF neutronsand the velocity of remaining components (neutron Bogoliubov thermal excitations, protons, and leptons) . As aresult, SF NSs host specific SF modes in addition to “normal” oscillation modes; the latter are close analogues ofoscillation modes in non-SF NSs. SF modes correspond to counter-motion of SF and normal fluid components, andhence, in contrast to normal modes, dissipate strongly due to powerful mutual-friction mechanism, that tends toequalize the velocities of these components [13–15]. In contrast to normal modes, the eigenfrequencies of SF modesstrongly depend on the stellar temperature (through the temperature dependence of neutron superfluid density). As aconsequence, avoided-crossings of normal and SF modes take place at certain (resonance) stellar temperatures. Nearthe resonances, eigenfunctions of strongly dissipating SF modes admix to those of normal modes and stabilize the Note that proton superconductivity does not lead to an additional independent velocity field, since protons are coupled to other chargedparticles by the electromagnetic forces. latter. Modeling (within the scenario of resonance stabilization of r-modes) shows that an NS in LMXB should spendmost of its life in the vicinity of such avoided-crossing [11, 12, 16, 17].Initially, this scenario was proposed as purely phenomenological one. To put it on a solid ground and to prove thatthe avoided-crossings of the most unstable r-mode with SF modes take place in the parameter range relevant to NSsin LMXBs, one has to calculate the r-mode spectrum for SF NSs at finite temperatures. This goal had been reachedin the previous studies [18, 19] under certain assumptions. Namely, Ref. [18] calculated the r-mode spectrum for aSF neutron star stratified by muons, neglecting the entrainment between neutrons and protons (i.e., assuming thatmotion of one particle species does not induce particle current of another species). Subsequent work [19] accounted forthe entrainment effect, but assumed that NS core consists of neutrons, protons, and electrons only. Here we calculatethe spectrum allowing for both muons and entrainment in the core, and show that together they change the spectrum qualitatively . Note that this paper is an extended version of the Letter [20], where we concentrate more on comparisonof our results with the available observations and constraining the neutron superfluidity model.In our numerical calculations we adopt up-to-date microphysics input. In this sense, although we still work inthe Cowling approximation and in the Newtonian framework when dealing with perturbed oscillation equations, theoscillation spectra calculated in this paper are expected to be realistic at least qualitatively, while the main conclusionswe arrived at in this work, we believe, are robust.In addition to the explanation of controversial observations of NSs in LMXBs, the scenario of resonance r-modestabilization proposes a new method to constrain the properties of superdense matter by finding the resonance tem-peratures in the r-mode spectrum and confronting them with available observations of NSs in LMXBs. Extremesensitivity of the calculated spectra to the model of neutron SF allowed us to constrain the latter.The paper is organized as follows. In Sec. II we provide the equations describing oscillations of rotating SF NSs.Sec. III discusses the expansion of these equations in the limit of slow rotation and weak entrainment. In Sec. IV wedescribe the microphysics input that was used in our numerical calculations. Sec. V presents the results. Finally, inSec. VI we discuss these results and conclude. Appendix contains a detailed analysis of the behavior of SF modes inthe limit of vanishing rotation rate.
II. OSCILLATION EQUATIONS
We consider oscillations of a slowly rotating NS with the spin frequency Ω = 2 πν . Dissipation is assumed to besmall and is neglected when calculating the spectra. We adopt Cowling approximation (do not account for metricperturbations [21]) and work in the Newtonian framework, when considering perturbed hydrodynamic equations. Inwhat follows we allow for muons ( µ ) in the inner layers of NSs, in addition to neutrons ( n ), protons ( p ), and electrons( e ) ( npeµ -composition), and also take into account possible SF of baryons (neutrons and protons) in the core. Let allthe quantities depend on time t as e ıσt in the coordinate frame rotating with the star. Then the linearized equationsgoverning small oscillations of SF NSs in that frame are:(i) Euler equation − σ ξξξ b + 2 ıσ ΩΩΩ × ξξξ b = δww ∇∇∇ P − ∇∇∇ δPw , (1)where w = ( P + ǫ ) /c , P is the pressure, ǫ is the energy density, c is speed of light. Note that Eq. (1) is not a purelyNewtonian one, it respects the fact that P can be comparable to ǫ in NS cores. Here and hereafter, δ stands for theEuler perturbation of some thermodynamic parameter (e.g., δP ). The Lagrangian displacement of baryons (vanishingin equilibrium), ξξξ b , in equation (1) is defined as ξξξ b ≡ jjj b ıσn b , (2)where n b ≡ n n + n p and jjj b ≡ jjj n + jjj p are the baryon number density and baryon current density, respectively; n i and jjj i are the number density and current density of particle species i = n, p, e, µ .(ii) Continuity equations for baryons and leptons δn b + div( n b ξξξ b ) = 0 , (3) δn l + div( n l ξξξ ) = 0 . (4)Here and hereafter, the subscript l = e, µ refers to leptons (electrons and muons); ξξξ ≡ jjj e / ( ıσn e ) is the Lagrangiandisplacement of the normal liquid component [we assume that all the normal-matter constituents (i.e., leptons andbaryon thermal excitations) move with one and the same normal velocity due to efficient particle collisions]. If neutronsare non-SF, then ξξξ = ξξξ b , and hydrodynamic equations become essentially the same as in the normal matter (even ifprotons are SF, see, e.g., Ref. [22]). Bearing this in mind we, for brevity, shall call “normal” (or “non-SF”) the liquidwith non-SF neutrons, irrespective of the actual state of protons.(iii) The “superfluid” equation, analogue of the Euler equation for SF (neutron) liquid component hσ zzz − ıh σ ΩΩΩ × zzz = c n e ∇∇∇ δη e + c n µ ∇∇∇ δη µ , (5)where zzz ≡ ξξξ b − ξξξ characterizes the relative Lagrangian displacement of SF and normal components; and η l ≡ µ n − µ p − µ l is the chemical potential imbalance (in equilibrium η l = 0 [23]). Further, h = µ n n b (cid:20) n b Y pp µ n ( Y nn Y pp − Y np ) − (cid:21) , (6) h = µ n n b (cid:18) n b Y nn µ n + Y np µ p − (cid:19) , (7)where Y ik = Y ki is the relativistic symmetric entrainment matrix [22, 24–26], which is the analogue of the SF mass-density matrix in the non-relativistic theory [27]. SF equation in the form (5) 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 typicalsituation in NSs (see, e.g., Refs. [28, 29]). The above equations should be supplemented with the relation betweenthe thermodynamic quantities, δn i = ∂n i ∂P δP + ∂n i ∂η e δη e + ∂n i ∂η µ δη µ , (8)where again i = n, p, e, µ . In what follows we shall use P , η e and η µ as independent thermodynamic variables.It is convenient to express the non-radial displacements ξ bθ , ξ bφ , z θ , and z φ as a sum of toroidal ( T b , T z ) and poloidal( Q b , Q z ) components [30]: ξ bθ = ∂∂θ Q b ( r, θ ) + ımT b ( r, θ )sin θ , ξ bφ = ımQ b ( r, θ )sin θ − ∂∂θ T b ( r, θ ) , (9) z θ = ∂∂θ Q z ( r, θ ) + ımT z ( r, θ )sin θ , z φ = ımQ z ( r, θ )sin θ − ∂∂θ T z ( r, θ ) , (10)where r and θ are the radial distance and polar angle in spherical coordinate system centered at the stellar center,with the axis z aligned with ΩΩΩ. Then, following the same procedure as for non-SF stars (e.g., [31]), we expand all theunknown functions into associated Legendre polynomials with fixed m : ξ br ( r, θ ) = ı X l ξ br l m ( r ) P ml (cos θ ) , (11) z r ( r, θ ) = ı X l z r l m ( r ) P ml (cos θ ) , (12) Q b ( r, θ ) = X l Q b l m ( r ) P ml (cos θ ) , (13) Q z ( r, θ ) = X l Q z l m ( r ) P ml (cos θ ) , (14) T b ( r, θ ) = X l T b l m ( r ) P ml (cos θ ) , (15) T z ( r, θ ) = X l T z l m ( r ) P ml (cos θ ) , (16) δP ( r, θ ) = X l δP l m ( r ) P ml (cos θ ) , (17) δη e ( r, θ ) = X l δη e l m ( r ) P ml (cos θ ) , (18) δη µ ( r, θ ) = X l δη µ l m ( r ) P ml (cos θ ) , (19)where the summation goes over l = m + 2 k and l = m + 2 k + 1 ( k = 0 , , , . . . ) for “odd” modes, and over l = m + 2 k + 1, l = m + 2 k for “even” modes .Let us consider a slowly rotating NS, and expand all the quantities in a power series in small parameter Ω [belowwe denote by Ω the rotation frequency normalized to the parameter Ω ≡ (cid:0) GM/R (cid:1) / , which is of the order of theKepler frequency; M and R are the stellar mass and radius, respectively]. We are interested in oscillations with theeigenfrequencies σ vanishing at Ω →
0. Thus, σ (normalized to Ω ) in the leading order in Ω can be represented as(e.g., [30, 31, 33]) σ = σ Ω, where σ does not depend on Ω.The equations, describing purely toroidal modes in the leading order in rotation, are given by: ∂∂θ (cid:0) sin θξ bθ (cid:1) + ımξ bφ = 0 , (20) σ ξ bθ + 2 ı cos θξ bφ = − ım ∂∂θ (cid:2) sin θ (cid:0) σ ξ bφ − ı cos θξ bθ (cid:1)(cid:3) , (21) ∂∂θ (cid:0) sin θz θ (cid:1) + ımz φ = 0 , (22) h ( r ) σ z θ + 2 ıh ( r )cos θz φ = − ım ∂∂θ (cid:8) sin θ (cid:2) h ( r ) σ z φ − ıh ( r )cos θz θ (cid:3)(cid:9) . (23)Here the index 0 indicates the leading-order term in the expansion of eigenfunctions in the power series in Ω. Thefirst couple of equations, Eqs. (20) and (21), describe the normal r-modes, analogous to ordinary r-modes in non-SFNSs, while Eqs. (22) and (23) describe SF modes driven by the relative motion (represented by the vector zzz ) of SFand normal (non-SF) liquid components [3, 15, 18, 34]. The solution to these two systems of equations allow us todetermine the eigenfrequencies σ = 2 ml ( l + 1) , (24) σ = 2 ml ( l + 1) h ( r ) h ( r ) (25)and eigenfunctions ξ bθ = ımT b lm ( r ) P ml (cos θ )sin θ , ξ bφ = − T b lm ( r ) ddθ P ml (cos θ ) , (26) z θ = ımT z lm ( r ) P ml (cos θ )sin θ , z φ = − T z lm ( r ) ddθ P ml (cos θ ) (27)of normal and SF modes, respectively.Since the function h ( r ) /h ( r ) in Eq. (25), generally, varies throughout the star, the frequency (25) cannot be aglobal oscillation frequency – each stellar layer has its own different eigenfrequency . However, if we assume that Y np = 0, then h ( r ) = h ( r ) [see equations (6)–(7)] and the eigenfrequency (25) of SF modes reduces to σ = 2 ml ( l + 1) , (28)becoming a global solution, independent of r [3, 15, 18, 34].As discussed in Ref. [18], for vanishing entrainment ( Y np = 0) in the lowest order in rotation, purely toroidal modesare only possible with l = m . For a given m the authors of Ref. [18] found one normal nodeless r -mode and an infiniteset of SF r -modes, all having the same σ = 2 / ( m + 1) [3, 15, 18, 34].However, when neutron and proton SFs co-exist somewhere in an NS, entrainment should be accounted for. Belowwe shall allow for entrainment ( Y np = 0) by considering it as a small perturbing parameter. Equations for odd and even modes completely decouple, thus odd and even modes do not mix with each other [32]. Similarly, oscillationequations (and hence the oscillation modes) completely decouple for different values of m . Except for some special cases when h ( r ) /h ( r ) is constant throughout the core. III. R-MODES IN THE LIMIT OF WEAK ENTRAINMENTA. Failed attempt
Assuming that the entrainment is weak, let us try to develop a perturbation scheme in small parameter ∆ h ≡ h/h − h → Y np →
0) in the leading order in rotation frequency. Our aim is to find the first-order correctionsin ∆ h to the eigenfrequency of r -modes in a SF npeµ NS. This approach is analogous to that of Ref. [19], where itwas shown that in a SF NS with npe -core r -modes can be calculated analytically in the first order in ∆ h .Below all the quantities are taken in the leading order in rotation. We denote the zeroth-order in ∆ h with theindex 0, and the first-order in ∆ h – with the index 1. Using this notation, the eigenfrequency and eigenfunctions canbe expanded in Taylor series in ∆ h as: σ = ( σ + σ )Ω = (cid:18) m + 1 + σ (cid:19) Ω , (29) ξ br = ξ br , T b = T b + T b , Q b = Q b , z r = z r , T z = T z + T z , Q z = Q z , (30) δP = δP Ω , δη e = δη e Ω , δη µ = δη µ Ω . (31)Since in the absence of entrainment the SF and normal r -modes are purely toroidal, the radial and poloidaldisplacements in the zeroth order vanish, ξ br = z r = Q b = Q z = 0. For rotational modes with σ ∝ Ω at Ω →
0, theEuler perturbation of any (scalar) thermodynamic parameter f (e.g., P , µ l , etc.) is proportional to Ω in the leadingorder in rotation and in the absence of entrainment (e.g., [14, 31, 33]). Following Ref. [19], we assume that non-vanishing entrainment does not change this ordering. We will need below only leading-order terms in the expansionsof scalar quantities in the rotation frequency Ω.In the zeroth order in ∆ h (i.e., at vanishing entrainment), as discussed above, the eigenfrequency σ of any r-modeequals σ = 2 m + 1 , (32)and the toroidal displacements are proportional to the l = m associated Legendre polynomial [3, 15, 18, 34], T b = T b mm ( r ) P mm (cos θ ) , T z = T z mm ( r ) P mm (cos θ ) . (33)The functions T b mm ( r ) and T z mm ( r ) cannot be found explicitly in the leading order in entrainment and rotationfrequency. Assuming vanishing entrainment, Ref. [18] proceeded to the next-to-leading order in rotation to calculate T b mm ( r ) and T z mm ( r ). In contrast, in Ref. [19] we worked in the leading order in rotation, and accounted forthe next-to-leading order terms in the entrainment to determine these functions. Here we follow the approach ofRef. [19]. To find the eigenfrequency correction σ and the functions T b mm ( r ) and T z mm ( r ), we shall consider thecontinuity equations (3)–(4), as well as r , φ , and θ -components of the Euler equation (1) and the SF equation (5).The θ -component of the Euler equation (combined with its φ -component) reads, in the first order in ∆ h (ignoringquadratically small terms such as σ ξ bθ ), σ ξ bθ + σ ξ bθ + 2 ı cos θξ bφ = 1 ım ∂∂θ (cid:8) sin θ (cid:2) σ ξ bφ + σ ξ bφ − ı (cid:0) ξ br sin θ + ξ bθ cos θ (cid:1)(cid:3)(cid:9) . (34)Substituting relations (9), (11), (13), (15) into equation (34) divided by sin θ and equating coefficients at the termsproportional to P mm , one can express the function Q b m +1 ,m ( r ) through ξ br m +1 ,m ( r ) and T b mm ( r ). Similarly, usingthe θ - and φ -components of the SF equation, σ z θ + σ z θ + ∆ hσ z θ + 2 ı cos θz φ = 1 ım ∂∂θ (cid:8) sin θ (cid:2) σ z φ + σ z φ + ∆ hσ z φ − ı (cid:0) z r sin θ + z θ cos θ (cid:1)(cid:3)(cid:9) , (35)one can obtain an algebraic relation between Q z m +1 ,m ( r ), z r m +1 ,m ( r ), and T z mm ( r ).Now, taking the coefficient at P mm +1 in the continuity equation for baryons,1 n b r ∂∂r (cid:0) r n b ξ br (cid:1) + 1 r sin θ (cid:20) ∂∂θ (cid:18) sin θ ∂Q b ∂θ (cid:19) − m Q b sin θ (cid:21) = 0 , (36) Note, however, that Ref. [19] defined ∆ h as ∆ h ≡ h /h − and expressing Q b m +1 ,m through T b mm and ξ br,m +1 ,m , we get the first-order ODE for ξ br,m +1 ,m : ddr ξ br,m +1 ,m + (cid:18) n ′ b n b + 3 + mr (cid:19) ξ br,m +1 ,m − (1 + m ) (3 + 2 m ) σ (2 + 4 m ) r T b mm = 0 , (37)where the prime denotes the derivative with respect to r .The continuity equations for electrons and muons,1 n e r ∂∂r (cid:2) r n e (cid:0) ξ br − z r (cid:1)(cid:3) + 1 r sin θ (cid:26) ∂∂θ (cid:20) sin θ (cid:18) ∂Q b ∂θ − ∂Q z ∂θ (cid:19)(cid:21) − m ( Q b − Q z )sin θ (cid:27) = 0 , (38)1 n µ r ∂∂r (cid:2) r n µ (cid:0) ξ br − z r (cid:1)(cid:3) + 1 r sin θ (cid:26) ∂∂θ (cid:20) sin θ (cid:18) ∂Q b ∂θ − ∂Q z ∂θ (cid:19)(cid:21) − m ( Q b − Q z )sin θ (cid:27) = 0 (39)in a stratified star [when d ( n e /n µ ) /dr = 0] imply ξ br − z r = 0 , (40) Q b − Q z = 0 . (41)Expressing Q b m +1 ,m and Q z m +1 ,m in Eq. (41) through, respectively, T b mm , ξ br,m +1 ,m and T z mm , z r,m +1 ,m , andusing Eq. (40), we find T z mm = (1 + m ) σ T b mm (1 + m ) σ + 2∆ h . (42)The φ -component of the Euler equation allows one to express δP m +1 ,m through T b mm , while the φ -component of theSF equation can be used to present δη µ m +1 ,m as a function of δη e m +1 ,m and T z mm .Substituting the obtained expressions for δP m +1 ,m and δη µ m +1 ,m into the r -component of the Euler equation, wederive an ODE of the form: ddr T b mm − mr T b mm + b ( r ) T z mm + c ( r ) δη e m +1 ,m = 0 , (43)while substitution of these expressions into the r -component of the SF equation gives ddr (cid:18) T z mm h n µ (cid:19) − mr T z mm h n µ + c ( r ) δη e m +1 ,m = 0 . (44)Equations (42), (43), and (44) allow us to express δη e m +1 ,m through T b mm , which results in the equation ddr T b mm − mr T b mm − σ a ( r, σ , ∆ h ) T b mm = 0 . (45)The functions b ( r ), c ( r ), c ( r ), and a ( r, σ , ∆ h ) in Eqs. (43)–(45) are known (have been found), but their actualform is not important for us here.Equations (37) and (45) describe the r-mode oscillations in the leading order in rotation frequency in SF npeµ NScore up to the first-to-the-leading-order correction in the entrainment under the assumption that expansions (29)–(31)are valid. This system, Eqs. (37) and (45), should be supplemented with a number of boundary conditions. Regularityin stellar center (at r →
0) requires T b mm ∝ r m , (46) ξ br,m +1 ,m = (1 + m ) σ m ) T b mm . (47)Since the particle number densities in our background model are continuous, the continuity equations imply con-tinuity of baryon and lepton radial displacements, ξ br,m +1 ,m and ξ r,m +1 ,m . This condition leads to the requirementof vanishing z r,m +1 ,m at the SF interface. Moreover, r and φ components of the Euler and SF equations requirecontinuity of the functions T b mm and T z mm throughout the star.At the stellar surface (where we assume that the matter is non-SF and barotropic) we require the Lagrangianperturbation of the pressure to be zero, ∆ P = 0, which means, in the leading order in rotation, ξ br,m +1 ,m ( R ) = 0 . (48)Consider now, for simplicity, a two-layer star composed of SF npeµ core and a barotropic single-fluid crust. To findthe eigenfunctions in the whole star we have to employ oscillation equations in the crust as well: ddr T b mm − mr T b mm = 0 , (49) ddr ξ br,m +1 ,m + (cid:18) n ′ b n b + 3 + mr (cid:19) ξ br,m +1 ,m − (1 + m ) (3 + 2 m ) σ (2 + 4 m ) r T b mm = 0 . (50)Assume first that σ = 0 (vanishing σ was a solution for normal r-mode in npe NSs, see Ref. [19]). Then equationsin the core and in the crust coincide. Moreover, ξ br,m +1 ,m = 0 throughout the star [due to Eqs. (37), (47), and (50)],while in the core z r,m +1 ,m = ξ br,m +1 ,m = 0 and T z mm = 0 due to, respectively, Eqs. (40) and (42). In addition, T b mm ∝ r m . This solution meets all the boundary conditions and describes the normal nodeless r-mode.Now, if σ = 0, then, integrating Eqs. (37), (45), and (49), (50), we have to meet three boundary conditions: in thestellar center (47), at the surface (48), and at the core-crust interface: ξ br,m +1 ,m = z r,m +1 ,m = 0. At the same time,we have only two integration constants (one of which defines the oscillation amplitude) and undefined value of σ . Thisis clearly not enough to meet all the boundary conditions; our system appears to be overdetermined. This happensbecause of restrictions (40) and (42). We come to conclusion that no oscillations at non-zero entrainment are possiblewith σ vanishing at Ω →
0, except for the normal r-mode. However, this conclusion looks to be unphysical, since itis hard to imagine that account for even infinitely small entrainment could eliminate SF modes with σ vanishing atΩ → B. Successful scheme
To demonstrate that such modes do exist, below we do not restrict ourselves to the leading order in rotationfrequency, but instead account for the next-to-the-leading-order corrections in rotation and in the entrainment si-multaneously . Such an approach allows us to relax the restrictions (40) and (42), and find a solution for SF r-modeoscillations at non-vanishing entrainment. In what follows we adopt the following expansions: σ = ( σ + σ )Ω = (cid:18) m + 1 + σ (cid:19) Ω , (51) ξ br = ξ br , T b = T b + T b , Q b = Q b , z r = z r , T z = T z + T z , Q z = Q z . (52)The leading order in both rotation and entrainment of each quantity is labeled with the index 0, while index 1 denotesnext to the leading order corrections, both in entrainment and rotation. For example, as we shall see below, σ behavesas Ω at high rotation frequency, and does not depend on the rotation frequency at small rotation rate and finiteentrainment. Since in the absence of entrainment the r -modes are purely toroidal in the leading order in rotation,the radial and poloidal displacements in the zeroth order vanish, ξ br = z r = Q b = Q z = 0. Equations describing theleading order are the same as in Sec. III A with the solution for the eigenfrequency and eigenfunctions given by Eqs.(32)–(33).To find σ and the functions T b mm ( r ) and T z mm ( r ), we again consider the continuity equations (3)–(4), as well as r , φ , and θ -components of the Euler equation (1) and the SF equation (5). As in Sec. III A, the θ -components of theEuler and SF equations allow us to express Q b m +1 ,m ( r ) through ξ br m +1 ,m ( r ) and T b mm ( r ), and Q z m +1 ,m ( r ) through z r m +1 ,m ( r ) and T z mm ( r ). The relations, however, slightly differ from those in Sec. III A, because here we accountfor the oblateness of an NS (since it is ∝ Ω ), as it is described in Ref. [18]. Oblateness gives rise to additional termsin the θ -components of the Euler and SF equations, as well as in the continuity equations, but does not affect ourapproach qualitatively. The r and φ -components of the Euler and SF equations remain unaffected and give us again[see Eqs. (43)–(44)] ddr T b mm − mr T b mm + b ( r ) T z mm + c ( r )Ω δη e m +1 ,m = 0 , (53) ddr (cid:18) T z mm h n µ (cid:19) − mr T z mm h n µ + c ( r )Ω δη e m +1 ,m = 0 . (54)The main difference from Sec. III A is the continuity equations. Now, in contrast to Sec. III A, they contain thenumber density perturbations. As a result, leptonic continuity equations remain nondegenerate [i.e., the constraints(40) and (41) should not be satisfied]. The coefficients at the polynomial P mm +1 in the three continuity equations (forbaryons, electrons and muons), as well as Eqs. (53)–(54) allow us to express an unknown function δη e m +1 ,m , andeventually arrive at the following system of oscillation equations: ddr T b mm − A ( r ) T b mm − B ( r ) T z mm − C ( r )Ω ( ξ br,m +1 ,m − z r,m +1 ,m ) = 0 , (55) ddr T z mm − A ( r ) T b mm − B ( r ) T z mm − C ( r )Ω ( ξ br,m +1 ,m − z r,m +1 ,m ) = 0 , (56) ddr ξ br,m +1 ,m − (cid:20) σ (1 + m ) (3 + 2 m )(2 + 4 m ) r + Ω A ( r ) (cid:21) T b mm − Ω B ( r ) T z mm − C ( r ) ξ br,m +1 ,m − D ( r ) z r,m +1 ,m = 0 , (57) ddr z r,m +1 ,m − Ω A ( r ) T b mm − (cid:20) (1 + m )(3 + 2 m )( σ + mσ + 2∆ h )(2 + 4 m ) r + Ω B ( r ) (cid:21) T z mm − C ( r ) ξ br,m +1 ,m − D ( r ) z r,m +1 ,m = 0 , (58)where A i ( r ), B i ( r ), C i ( r ), D i ( r ) ( i = 1 , , ,
4) are some known functions of the radial coordinate. Regularity of theseequations at the stellar center ( r →
0) implies T b mm , T z mm , ξ br,m +1 ,m , z r,m +1 ,m ∝ r m , (59) ξ br,m +1 ,m = (1 + m ) m σ T b mm , (60) z r,m +1 ,m = 1 + m m ( σ + mσ + 2∆ h ) T z mm . (61)Vanishing Lagrangian perturbation of the pressure at the surface leads to the condition T b mm ( R ) = (1 + m ) (1 + 2 m ) ξ br,m +1 ,m ( R ) P ′ ( R )4 m Ω w ( R ) R . (62)One also has to require vanishing of z r,m +1 ,m at the SF interface and continuity of the functions ξ br,m +1 ,m , T b mm ,and T z mm throughout the star.Consider again, for simplicity, a two-layer star, consisting of the SF npeµ core and a barotropic single-fluid crust.Having six integration constants (two in the crust and four in the core), we can meet all the necessary boundaryconditions [(60), (61), (62), continuity of ξ br,m +1 ,m and T b mm , and vanishing of z r,m +1 ,m at the core-crust interface),by adjusting the value of σ . Thus, the approach developed in this section allows us to restore the missing SFrotational modes in the spectrum. C. Limit of vanishing rotation rate
Let us now analyze the solution to the system (55)–(58) in the limit Ω →
0, which got us into trouble in Sec.III A. If T z mm = 0, then ∆ h does not enter the equations and we have the standard ordering, σ ∼ Ω , ξ br,m +1 ,m ∼ z r,m +1 ,m ∼ Ω T b mm , δP ∼ δη e ∼ δη µ ∼ Ω . This is the normal r-mode. For SF modes T z mm = 0 and the standardordering is not valid [see the term ∝ T z mm in Eq. (58)]. Instead, as we demonstrate below, σ appears to be finiteand the ordering of eigenfunctions is the following: ξ br,m +1 ,m ∼ z r,m +1 ,m ∼ Ω T b mm ∼ Ω T z mm ∼ Ω dT b mm /dr ∼ Ω dT z mm /dr ∼ Ω dξ br,m +1 ,m /dr ∼ Ω dz r,m +1 ,m /dr ∼ δη e ∼ δη µ ∼ δP/ Ω.To start with, we rewrite Eqs. (55)–(58) assuming the ordering above and neglecting the terms, which are smallaccording to this ordering. Doing this, we step away from the stellar center, in order not to deal with the additionalsmall parameter, r : ddr T b mm − C ( r )Ω ( ξ br,m +1 ,m − z r,m +1 ,m ) = 0 , (63) ddr T z mm − C ( r )Ω ( ξ br,m +1 ,m − z r,m +1 ,m ) = 0 , (64) ddr ξ br,m +1 ,m − σ (1 + m ) (3 + 2 m )(2 + 4 m ) r T b mm = 0 , (65) ddr z r,m +1 ,m − (1 + m )(3 + 2 m )( σ + mσ + 2∆ h )(2 + 4 m ) r T z mm = 0 . (66)The above equations can be combined to (again, up to the terms, leading in Ω) dTdr − K ( r, σ )Ω ξ = 0 , (67) dξdr − T = 0 , (68)where we define ξ ≡ ξ br,m +1 ,m − z r,m +1 ,m , (69) T ≡ σ (1 + m ) (3 + 2 m )(2 + 4 m ) r T b mm − (1 + m )(3 + 2 m )( σ + mσ + 2∆ h )(2 + 4 m ) r T z mm , (70) K ( r, σ ) ≡ σ (1 + m ) (3 + 2 m )(2 + 4 m ) r C ( r ) − (1 + m )(3 + 2 m )( σ + mσ + 2∆ h )(2 + 4 m ) r C ( r ) . (71)Writing down the ratio of Eqs. (67) and (68), T dT = K ( r, σ )Ω ξdξ, (72)we find that ξ ∼ Ω T , or ξ br,m +1 ,m ∼ z r,m +1 ,m ∼ Ω T b mm ∼ Ω T z mm [from Eqs. (63)–(64) it follows that T b mm ∼ T z mm , while Eqs. (65)–(66) imply that ξ br,m +1 ,m ∼ z r,m +1 ,m ]. Notably, this ordering is different from the standardone, which takes place at vanishing entrainment, when we account for the rotational corrections only. In that case∆ h = 0, σ ∝ Ω and one finds from Eqs. (55)–(58) ξ br,m +1 ,m ∼ z r,m +1 ,m ∼ Ω T b mm ∼ Ω T z mm .On the other hand, Eqs. (67)–(68) can be rewritten as d ξdr − K ( r, σ )Ω ξ = 0 . (73)If K ( r, σ ) < r , we have an oscillating solution there (with a vanishing wavelength at Ω → K ( r, σ ) > → K ( r, σ ) < →
0, since they contain an infinite number of nodes. At the same time,we also cannot have K ( r, σ ) > npeµ core and in the remaining star. The solution is only possibleif K ( r, σ ) > → K ( r, σ ) vanishes (seeAppendix for details).As a result, σ for SF modes tends to the finite value at Ω →
0, which can be found from the conditionmin[ K ( r, σ )] = 0 , (74)i.e., the minimum of the function K ( r, σ ) in the SF npeµ -core must vanish. Moreover, all the overtones must have thesame σ [since infinitesimal variation of σ allows one to increase the number of nodes in the infinitely small region, In our analysis we also assume that neutrons are SF,
T < T c n , everywhere in the npeµ core. This assumption allows us to avoid thesingularity related to infinite growth of the function h at the boundary of the SF region, at T → T c n . Note that h → ∞ at T → T c n ,and the excluded terms can be large. However, numerical calculations show that the consideration below is relevant at any temperature. K ( r, σ ) < K ( r, σ ) may have the minimum at any point in SF npeµ core, in particular, inthe centre or at the muon onset density. For example, in the zero-temperature limit the minimum occurs at the outerboundary of the npeµ core for APR EOS (see Sec. IV), while for BSk24 EOS it is located at the stellar center.It is now interesting to discuss why our attempt to find the solution in Sec. III A was not successful. First of all, inSec. III A we assumed the standard ordering of eigenfunctions, which is not valid for SF modes, as we demonstrated inthis section. While the φ -component of the Euler equation implies δP ∼ Ω T b , the imbalances of chemical potentialsscale as δη l ∼ Ω T b due to the scaling ξ br ∼ z r ∼ Ω T b . As a result, perturbations of particle number densities scale as δn i ∼ Ω T b . At first glance, it seems that with this ordering we can skip δn l in Eqs. (38) and (39) in the leading orderin Ω, because dξ br /dr ∼ T b . However, if we consider the difference of Eqs. (38) and (39), dξ br /dr and the functions Q b and Q z will cancel out, leaving us with δn e − δn µ + (cid:0) ξ br − z r (cid:1) (cid:18) n e ∂n e ∂r − n µ ∂n µ ∂r (cid:19) = 0 . (75)All the terms in this equation are of the same order in Ω. Thus, the constraints (40) and (41) of Sec. III A are notapplicable in our situation. IV. PHYSICS INPUT
In our numerical calculations we adopt a model of a three-layer star, consisting of the crust, outer core (composed ofneutrons, protons, and electrons), and the inner core (containing, in addition, muons). Since the crust does not affecteigenfrequencies of global NS oscillations strongly, it is treated within a simplest model of barotropic one-componentfluid. The core is modeled more realistically, adopting modern equations of state (EOSs), which lead to non-barotropicbehavior of the core liquid due to composition gradients, and accounting for possible superfluidity/superconductivityof baryons.We consider two EOSs in the core. The first one is essentially the same as in Ref. [18]. It uses parametrization[35] of Akmal-Pandharipande-Ravenhall (APR) EOS [36], and adopts entrainment matrix from Ref. [37]. The secondEOS is constructed with the BSk24 energy-density functional (BSk24 EOS) [38, 39]. The entrainment matrix for thisEOS is calculated self-consistently, by extracting nucleon Landau parameters from this functional and then followingthe prescription of Refs. [24–26] (see Ref. [40] for details).To illustrate the behavior of the entrainment matrix we plot its elements Y ik ( i, k = n, p ) as functions of densityfor two EOSs in the limit of vanishing temperature (see solid lines in the left and right panels of Fig. 1). In bothpanels density ranges from its value at the core-crust interface up to the central density in the limiting configurationof an NS with the maximum mass. Vertical dots show central densities of an NS with M = 1 . M ⊙ . In the same plotwe also present by dashes the parameter ∆ h . One can see that ∆ h can hardly be considered small at high densities.However, since the r-mode eigenfunctions are mostly localized in the outer core, in what follows we shall treat ∆ h as a small parameter and use expansions discussed in the previous section. Note that the method of expansion inthe entrainment has been proved to be rather accurate in Ref. [19], where the rotational spectrum of an NS withthe superfluid npe core was studied. In that reference we calculated the (temperature-dependent) eigenfrequencies ofsuperfluid r-mode by two different methods, either using the perturbation theory in the entrainment or not using it(exact calculation). Both approaches give very similar temperature-dependent spectra (compare upper solid red lineand dot-dashed line in figure 2 of Ref. [19]).For both EOSs we consider three SF models. In the first one neutron and proton local critical temperatures of SFonset are density independent (model A) and equal T c n = 6 × K and T c p = 5 × K for neutrons and protons,respectively. Although this model, as it is, is not realistic from point of view of microscopic calculations, it can beviewed as a limiting case for realistic wide critical temperature profiles. In our second and third models (models B andC) we adopt density-dependent (local) critical temperature profiles, T c n ( n n ) and T c p ( n p ) (see Fig. 2). These modelscoincide with the models I and II from Ref. [20]. The T c p profile is the same for both models, while T c n profiles differ.Model B has a wide T c n profile, similar to that predicted by modern microscopic calculations (see, e.g., [41, 42]). ModelC, on the contrary, describes a narrow T c n profile, which leaves the outer core non-SF at any reasonable temperature.Similar (purely phenomenological) profiles have been used in a number of works (e.g., [43–46]) to successfully explainthermal properties of isolated neutron stars within the minimal cooling scenario [43, 47]. All the adopted profiles havemaximum critical temperatures that do not contradict both the existing data on cooling NSs [43–50] and microscopiccalculations [41, 42, 51–54].1 Y i k / ( c m − e r g − ) log ρ, g cm − T = 0APR EOS Y np Y pp Y nn ∆ h M = . M ⊙ ∆ h log ρ, g cm − T = 0BSk24 EOS Y np Y pp Y nn ∆ h M = . M ⊙ FIG. 1: Entrainment matrix elements Y ik and the parameter ∆ h versus density in the zero-temperature limit for APR EOS(left panel) and BSk24 EOS (right panel). V. RESULTS
We present the results for the l = m = 2 r-modes, since l = m = 2 normal r-mode without nodes is known to bethe most unstable one [55]. Fig. 3 shows how the correction σ depends on the rotation frequency ν for an NSwith M = 1 . M ⊙ and the stellar red-shifted temperature T ∞ = 10 K (as seen by a distant observer). To plot thefigure, we adopted BSk24 EOS and SF models from Sec. IV. Left panel demonstrates the spectrum for models A andB. Since for these models T ∞ c n ≫ K everywhere in the core (for the chosen stellar configuration), the spectrumfor them is indistinguishable. For these models we find one normal nodeless r-mode ( n ), one SF nodeless r-mode( s ), and an infinite set of SF modes with nodes (only first two overtones, s and s , are presented in the figure).Various oscillation modes are shown by solid lines; in the inset these lines are shown by different colors. Note thatthe modes exhibit avoided-crossings with each other, altering their behavior from normal-like to SF-like and viceversa. Dots (red online) indicate normal nodeless r-mode (at different temperatures different r-modes behave as thenormal one). The normal r-mode is virtually not affected by entrainment, while SF modes deviate strongly from their‘vanishing-entrainment’ behavior (shown by dashes), especially at small rotation frequencies. Diamond at ν = 0 showsthe theoretically predicted limit of σ at ν → σ ∝ Ω , as expected.Right panel shows the spectrum obtained for the model C. For this model, in addition to nodeless r-modes (normaland SF), we find both an infinite set of SF modes with nodes and an infinite set of normal modes with nodes. Theplot shows the main harmonics ( n and s ) and first two overtones ( n , n , and s , s ) of normal ( n ) and SF ( s )modes. The modes exhibit avoided-crossings with each other. Irregularities of the curves at low frequencies are dueto the avoided-crossings with higher-order SF modes. Again, the diamond at ν = 0 shows the theoretically predictedlimit (defined by min [ K ( r, σ )] = 0 , see Eq. 74) for σ at ν → σ for the normal modes tends to zero at ν →
0, as In the left panel the rotation frequency ranges up to unphysical values of the order of ν ≈ For BSk24 EOS and T ∞ = 10 K, K ( r ) has a minimum at the stellar center for all SF models. T c n / K n n , fm − m o d e l B m o d e l C T c p / K n p , fm − FIG. 2: Left panel: Neutron (local) critical temperatures T c n versus n n for models B and C. Right panel: Proton (local) criticaltemperature T c p ( n p ) (the same for models B and C). Vertical lines (dashes for APR EOS, dots for BSk24 EOS) indicate centralnumber densities of (from left to right) M = 1 . M ⊙ , M = 1 . M ⊙ , and M = 1 . M ⊙ NSs. expected (see Sec. III).One can see that the spectra in both panels differ qualitatively. The SF models A and B do not support normalmodes with nodes, while model C supports them. This happens because for models A and B at T ∞ = 10 K there isno non-SF non-barotropic region in the star: the crust is assumed to be barotropic and the whole core is SF. At thesame time, in model C the outer core is normal, and, since we use the non-barotropic EOS in the core, the outer coremay support normal r-modes with nodes [32, 33]. In contrast to normal modes, we have an infinite set of SF modeswith nodes for all the three models, because for all of them there is a SF region stratified by muons at T ∞ = 10 K.In Fig. 3 one can see avoided-crossings of normal nodeless r-mode with SF r-modes. At the rotation frequenciescorresponding to the avoided-crossings, ν n ,s α ( α = 0 , , , . . . ), one should expect stabilization of normal r-modeby the resonance interaction with SF r-modes [11, 12]. At different temperatures the position ν n ,s α of the avoidedcrossings will be different. Fig. 4 illustrates how ν n ,s α depend on T ∞ for APR EOS (left panel) and BSk24 EOS(right panel). Solid lines correspond to model C, dots represent model B, dashes correspond to the model A. Thicklines (both solid lines, dashes, and dots) in both panels show the results for an NS with M = 1 . M ⊙ , thin lines – for M = 1 . M ⊙ . Each curve is marked with the sign n , s or n , s , which correspond to the avoided-crossings of normalnodeless r-mode with, respectively, the main harmonic or first overtone of SF r-modes. For comparison, points witherror bars in Fig. 4 show the positions of the observed NSs in LMXBs taken from Refs. [12, 56]. Here we do notindicate the names of these sources to avoid cluttering of the figure, one can find the names in Fig. 5.Generally (at not too high temperatures, see below), n , s avoided crossing takes place at unphysically highrotation rates. Only when T ∞ approaches the maximum value of T ∞ c n in the region of the core, where neutron andproton SFs co-exist (denoted by T ∞ c n max in what follows), ν n ,s ( T ∞ ) decreases rapidly with increasing T ∞ and fallsto zero at T ∞ = T ∞ cn max . While the local value of T c n max equals T c n max = 6 × K for all our SF models, thered-shifted value, T ∞ c n max , depends on the NS mass and EOS through the redshift parameter and varies in the range T ∞ c n max ∼ (2 . − × K. As a result, we have almost vertical drop of ν n ,s at T ∞ ∼ (2 − × K (see Fig.4) . The exception is the low-mass configurations in the model C, which have low values of T ∞ c n max (see Fig. 2).Dot-dashed line in the right panel illustrates this point, corresponding to the n , s avoided-crossing for M = 1 . M ⊙ NS in the model C. We do not plot the curves for an NS with M = 1 . M ⊙ to avoid cluttering of the figure; they are very similar to those plotted for an NSwith M = 1 . M ⊙ . σ ν, Hz BSk24 EOS M = 1 . M ⊙ T ∞ = 10 Kmodels A, B s s s n s s s s s n n ν, Hz M = 1 . M ⊙ T ∞ = 10 KBSk24 EOSmodel C n s n n s s s s FIG. 3: σ versus ν at T ∞ = 10 K for models A and B (left panel) and model C (right panel). The spectrum is plotted for anNS with M = 1 . M ⊙ assuming BSk24 EOS. See text for details. The n , s avoided-crossing occurs at lower rotation frequencies than the n , s one. For the models A and B(representing wide critical temperature profiles) the corresponding n , s curves pass through the sources in theinstability window, and some stellar models (e.g., high-mass configurations for APR EOS), allow one to explainnot too hot sources. At the same time, in the model C, n , s avoided-crossing lies at much lower rotation rates.This happens because in this model T c n profile drops sharply in the outer core, shrinking the SF region even at lowtemperatures. As we checked for various SF models, such shrinking of SF region in the outer core leads to a dramaticdecrease of ν n ,s at a given T ∞ . Analogous shrinking of SF region due to drop of T c n at its higher-density slope (inthe stellar center) leads to the same effect, which is, however, not so dramatic.In the vicinity of the avoided-crossings normal r-mode exhibits stabilizing resonance interaction with the SF r-modes[11, 12]. This leads to formation of “stability strips” along the ν n ,s α ( T ∞ ) curves in the ν − T ∞ plane (termed stabilitypeaks in the initial scenario of Refs. [11, 12]). Fig. 5 shows two examples of the instability windows with such stripsfor an NS with M = 1 . M ⊙ . The figure is plotted for the models A (left panel) and B (right panel), assuming APREOS. To calculate these windows we accounted for the shear viscosity and mutual friction dissipation, as described inRef. [18]. Note that our results imply that stability peaks are not vertical, as simple model of Refs. [11, 12] suggested.Nevertheless, according to the scenario of Refs. [11, 12], the star in the course of its evolution in LMXB should stillspend most of its life climbing up the left edge of the “peak”. Thus, all the observed sources should be located withinthe stability region, which, however, varies with the stellar mass. VI. DISCUSSION AND CONCLUSIONS
In order to confirm the viability of the phenomenological scenario of r-mode resonance stabilization suggested inRefs. [11, 12] and to put it on a solid ground, we have calculated temperature-dependent r-mode spectrum for a slowlyrotating SF NS. In contrast to the previous studies [18, 19], we, for the first time, accounted for both entrainmentbetween neutrons and protons ( Y np = 0) and the presence of muons in the NS core. Both muons and entrainmentaffect the spectrum qualitatively. Accounting for these effects together leads to elimination of rotational SF modesfrom the oscillation spectrum, if one works in the leading order in rotation and assumes the standard ordering ofeigenfunctions in Ω. This unphysical elimination takes place because non-vanishing entrainment requires non-vanishingradial displacement of the normal fluid component, while stratification by muons forbids such displacements in theleading order in rotation. To restore SF r-modes and calculate the spectrum we developed an original perturbation4 ν [ H z ] log T ∞ model B1 . M ⊙ model A1 . M ⊙ model A1 . M ⊙ model B1 . M ⊙ model B1 . M ⊙ model C1 . M ⊙ model C1 . M ⊙ APR EOS n , s n , s n , s n , s n , s n , s n , s log T ∞ model C1 . M ⊙ model A1 . M ⊙ model A1 . M ⊙ model B1 . M ⊙ model B1 . M ⊙ model C1 . M ⊙ model C1 . M ⊙ BSk24 EOS n , s n , s n , s n , s n , s n , s n , s FIG. 4: Values of T ∞ and ν at which normal r-mode experiences avoided-crossing with the main harmonic (marked with n , s )and first overtone (marked with n , s ) of SF r-modes. Thick lines (dashes for the model A, dots for the model B, and solidlines for the model C) show results for an NS with M = 1 . M ⊙ , while thin lines (again dashes for the model A, dots for themodel B, and solid lines for the model C) correspond to an NS with M = 1 . M ⊙ . Left panel: APR EOS, right panel: BSk24EOS. In the vicinity of these curves normal r-mode experiences stabilizing interaction with corresponding SF r-mode (see Fig.5). Dot-dashed line in the left panel shows ν n ,s ( T ∞ ) for the limiting configuration of an NS ( M = 1 . M ⊙ ) in the model B.In the right panel dot-dashed line shows ν n ,s ( T ∞ ) for M = 1 . M ⊙ NS in the model C. Region filled gray is classical stableregion defined by the shear viscosity only (calculated as described in Ref. [11] for an NS with M = 1 . M ⊙ ). Points with errorbars describe available observational data on NSs in LMXBs [12, 56]. scheme, whose leading order corresponds to the leading order in rotation and vanishing entrainment, while the next-to-the-leading order includes corrections due to non-vanishing entrainment (which is treated as a small parameter)and rotation simultaneously .Using this perturbation scheme we proposed a simple method to calculate the eigenfrequency of SF r-modes [namely,the first correction to the known leading-order value σ = 2 / ( m + 1)] in the limit of vanishing rotation rate. Wedemonstrated that it can be calculated from the simple algebraic equation (74).Then, applying this scheme to the more general situation when the rotation frequency is small but non-vanishing,we have calculated temperature-dependent r-mode spectrum of an NS for two EOSs and three SF models (A, B, andC). We found that the normal nodeless l = m = 2 r-mode (that is, the most unstable one [55]) exhibits avoided-crossings with the SF r-modes. Near the avoided-crossings normal r-mode experiences stabilizing interaction with SFmodes, which suppresses the r-mode instability. When the avoided-crossing takes place in the region of the classicalinstability window, it results in the stability strip in the “stellar temperature – rotation frequency” plane (termed“stability peak” in the initial phenomenological scenario [11, 12]). Although the strips, calculated here, are not verticalas suggested in Refs. [11, 12], an NS in LMXB will spend the majority of time climbing up the left edge of such strip[11, 12, 16, 17], according to the resonance stabilization scenario. Our calculations demonstrate that for certain SFmodels avoided-crossings take place in the range of parameters relevant to NSs in LMXBs, falling within the classicalinstability window (see Fig. 4). This puts the scenario of Refs. [11, 12] on a solid ground, making it a quantitativetheory, which can explain observations of hot, rapidly rotating NSs in LMXBs.Our calculations open up a possibility to constrain parameters of neutron SF. First, one can estimate T c n max bynoticing that the hottest sources can be stabilized by the resonance with the main harmonic of SF r-mode ( n , s ).For both considered EOSs, all the stellar masses, and all SF models considered by us (not only the models A, B, Cdiscussed in the paper) the curve ν n ,s ( T ∞ ) rapidly drops to zero at T ∞ = T ∞ c n max . In order for the curve to passthrough the hottest sources in Fig. 4, one should require T c n max ∼ (3 − × K (see Ref. [20] for details). Notethat this estimate is the lower limit for the maximal T c n , since T c n max denotes the maximum of T c n in the region,5
4U 1608-522SAX J1750.8-2900IGR J00291-5934MXB 1659-298 EXO 0748-676Aql X-1KS 1731-260SWIFT J1749.4-2807SAX J1748.9-2021XTE J1751-305SAX J1808.4-3658 IGR J17498-2921HETE J1900.1-2455XTE J1814-338IGR J17191-2821IGR J17511-3057NGC 6440 X-2XTE J1807-294XTE J0929-314Swift J1756-2508
4U 1608-522SAX J1750.8-2900IGR J00291-5934MXB 1659-298 EXO 0748-676Aql X-1KS 1731-260SWIFT J1749.4-2807SAX J1748.9-2021XTE J1751-305SAX J1808.4-3658 IGR J17498-2921HETE J1900.1-2455XTE J1814-338IGR J17191-2821IGR J17511-3057NGC 6440 X-2XTE J1807-294XTE J0929-314Swift J1756-2508 ν [ H z ] log T ∞ M = 1 . M ⊙ APR EOSmodel A log T ∞ M = 1 . M ⊙ APR EOSmodel B
FIG. 5: Instability windows for l = m = 2 normal nodeless r-mode calculated for M = 1 . M ⊙ NS with APR EOS and the SFmodels A and B (left and right panels, respectively). In the region filled gray NS is stable. At low frequencies high-temperaturestability peak in the model A is not calculated due to numerical problems. Dotted lines show ν n ,s ( T ∞ ) (the same lines as inFig. 4). We do not plot the curves ν n ,s ( T ∞ ) here, but they follow the corresponding stability peaks. See text and Fig. 4 forfurther details. where both neutrons and protons are SF. Our constraint is consistent with microscopic calculations [41, 42, 51–54],as well as with observations of cooling NSs [43–50].Second, the resonance with the first overtone of SF r-modes ( n , s ) can stabilize less hot rapidly rotating NSs forour models A and B (wide neutron critical temperature profiles, most of the stellar core is SF), but not for the modelC (narrow profile). In the case of narrow T c n profiles, resulting in much lower values of ν n ,s , these NSs can bestabilized only by the ( n , s ) resonance, if they have sufficiently low masses (see dot-dashed line in the right panel ofFig. 4). However, this explanation is less plausible since NSs in LMXBs are generally believed to have high masses[57, 58]. It is also worth mentioning that if T c n max > K and T c n profile is wide, i.e., T c n remains high in the wholeNS core, then hottest rapidly rotating sources, could be stabilized by the resonance with the first overtone of SFr-mode. Unfortunately, in this case we do not see a possibility to stabilize other, moderately heated, sources withinour minimal scenario.Resuming, to explain all wealth of observational data on NSs in LMXBs within the scenario of resonance r-modestabilization [11, 12] (which we consider as a minimal extension of the classical scenario), one needs a wide neutroncritical temperature profile (see above), such that in the region of co-existence of neutron and proton SFs the maximumvalue of T c n is T c n max ∼ (3 − × K. The real maximum of T c n throughout the whole density range can be largerthan T c n max . Note that proton SF model cannot be constrained in our scenario since proton pairing only weaklyaffects the oscillation modes [22, 59, 60].In principle, the positions of stability strips are sensitive not only to the neutron SF model, but also to EOS: as onecan see from Fig. 4, APR and BSk24 EOSs yield the results that differ by tens of percent. However, in order to obtainconstraints on EOS based on the proposed mechanism, one has to calculate the r-mode spectrum more accurately, byallowing for gravitational field perturbations, higher-order terms in the expansions (51)–(52), as well as the GeneralRelativity effects.Definitely, one should also keep in mind that the real instability window may also be affected by additional r-modestabilization mechanisms such as Ekman layer dissipation [61, 62], bulk viscosity in hyperon/quark matter [63–65],enhanced mutual friction dissipation [66] etc., which are not considered in the present paper. Accounting for all theseeffects can, in principle, modify our constraints on the parameters of neutron SF.6 Acknowledgments
This work is supported by the Russian Science Foundation (grant number 19-12-00133).
Appendix A: Constructing a solution for superfluid modes in the limit of vanishing rotation rate
To start with, we again, for simplicity, consider a two-layer star consisting of the SF npeµ core and the crust. Morerealistic structure with npe -layer in between is discussed in the end of this Appendix, and analyzed by us numericallyas well. It shows the same qualitative behavior as the simplified two-layer stellar model.Let us try to build a solution for SF modes in the limit Ω →
0. Integrating oscillation equations in the crust, we findeigenfunctions ξ br,m +1 ,m and T b mm , which have a standard ordering in the crust: ξ br,m +1 ,m ∼ Ω T b mm . Generally,these functions do not vanish at the core-crust interface ( r = R cc ). This means that, due to the boundary conditions atthe interface (continuity of ξ br,m +1 ,m and T b mm ), ξ br,m +1 ,m and T b mm should also be ordered as ξ br,m +1 ,m ∼ Ω T b mm in the core at r = R cc . However, as it was shown in Sec. III C, ξ ∼ Ω T . Generally this means that, in the SF npeµ core, ξ br,m +1 ,m ∼ Ω T b mm [see Eqs. (69) and (70)], and as a result, the boundary conditions at the core-crust interfacecannot be satisfied.There is, however, a possibility to solve this problem. The idea is to construct such a solution in the core that ξ and thus ξ br,m +1 ,m (note that ξ = ξ br,m +1 ,m at r = R cc ) is suppressed by a factor of Ω (or, in other words, almostvanish) at r = R cc , while T is not suppressed. Such solution satisfies the desired ordering, ξ br,m +1 ,m ∼ Ω T b mm .Continuity of T b mm can be ensured by an appropriate choice of the normalization constant, while slightly varying σ one can choose ξ br,m +1 ,m at r = R cc in such a way to satisfy the continuity of ξ br,m +1 ,m at the core-crust interface.To verify if we indeed can construct such a solution, we should solve Eq. (73): d ξdr − K ( r, σ )Ω ξ = 0 . (A1)Strictly speaking, Eq. (A1) is not valid in the very vicinity of the stellar center (because when deriving it, we assumedthat r is not small). On the other hand, the asymptotic formulas, Eqs. (59)–(61), are valid in the very vicinity of thecenter. Let us denote by r c a coordinate such that the asymptotes (59)–(61) are valid at r ≤ r c , while Eq. (A1) isvalid at r c ≤ r ≤ R cc . The value of r c can, in principle, be determined from Eqs. (55)–(58) by analyzing behavior ofvarious terms at Ω → r →
0. In the limit of vanishing rotation rate r c →
0. In what follows, we make use ofthe asymptotic formulas, Eqs. (59)–(61) at 0 ≤ r ≤ r c , and solve Eq. (A1) at r c ≤ r ≤ R cc .Note that Eq. (A1) has the same form as the Schr¨odinger equation. The limit Ω → K ( r, σ ) vanishes] can be writtenas [67] ξ = A K ( r, σ ) / exp " R rr c p K (˜ r, σ ) d ˜ r Ω + A K ( r, σ ) / exp " − R rr c p K (˜ r, σ ) d ˜ r Ω , (A2)where A and A are some constants.Consider first a situation when K ( r, σ ) > npeµ core (at r c ≤ r ≤ R cc ). Then we have no“turning points” and the solution (A2) is valid throughout the core. Let us check if we can meet all the boundaryconditions in this case.The asymptotic solution (59)–(61) implies that [see also Eqs. (69) and (70)] T = 3 + 2 mr c ξ (A3)at 0 ≤ r ≤ r c . Keeping in mind that T = dξ/dr [see Eq. (68)] at r ≥ r c , and making use of the continuity ofeigenfunctions ξ br,m +1 ,m , z r,m +1 ,m , T b mm , and T z mm (and thus the continuity of T and ξ ) at r = r c , we find dξdr = 3 + 2 mr c ξ (A4)at r = r c , or in view of Eq. (A2):3 + 2 mr c ( A + A ) = ( A − A ) p K ( r c , σ )Ω −
14 ( A + A ) K ′ ( r c , σ ) K ( r c , σ ) ≈ ( A − A ) p K ( r c , σ )Ω , (A5)7where the prime denotes the derivative with respect to r . One can see that the constants A and A are of the sameorder of magnitude. Thus, at the core-crust interface the decreasing exponent in the solution (A2) will be negligiblysmall and we find ξ = A K ( R cc , σ ) / exp " R R cc r c p K (˜ r, σ ) d ˜ r Ω , (A6) T = A K ( R cc , σ ) / Ω exp " R R cc r c p K (˜ r, σ ) d ˜ r Ω (A7)at r = R cc with ξ ∼ Ω T and no possibility to match eigenfunctions at the core-crust interface (see above).Assume now that K ( r, σ ) vanishes at some radius, then the solution (A2) is invalid in the vicinity of this radius[67] and consideration above is not applicable. Consider a situation when K ( r, σ ) has a negative minimum at someradius r such that r c < r < R cc , K ( r , σ ) = − K . The function K ( r, σ ) can be expanded near r as K ( r, σ ) = − K + β ( r − r ) , (A8)with β ≡ K ′′ ( r, σ ) | r / >
0, and Eq. (A1) takes the form d ξdx − (cid:18) a + 14 x (cid:19) ξ = 0 , (A9)where a ≡ − K / (2Ω β / ), x ≡ (4 β ) / ( r − r ) / √ Ω. Note that, in the limit Ω → x changes from −∞ to + ∞ . Thesolution to Eq. (A9) can be expressed in terms of parabolic cylinder functions [68]. General solution at x ≥ U ( a, x ) and V ( a, x ) introduced in Ref. [68], ξ = B U ( a, x ) + B V ( a, x ),and exhibiting the following asymptotic behavior at x ≫ | a | U ( a, x ) ≈ exp (cid:18) − x (cid:19) x − a − , (A10) V ( a, x ) ≈ r π exp (cid:18) x (cid:19) x a − . (A11)Notice that, since Eq. (A9) is symmetric with respect to the transformation x → − x , the functions U ( a, − x ) and V ( a, − x ) are also solutions to Eq. (A9). Since we are interested in real ξ , we shall present the solution to Eq. (A9)in the region x ≤ U ( a, − x ) and V ( a, − x ), ξ = B U ( a, − x ) + B V ( a, − x ), withthe following asymptotic behavior at − x ≫ | a | U ( a, − x ) ≈ exp (cid:18) − x (cid:19) ( − x ) − a − , (A12) V ( a, − x ) ≈ r π exp (cid:18) x (cid:19) ( − x ) a − . (A13)The constants B , B , B , and B should be found by matching the functions ξ and T in different regions describedbelow.Let us consider three regions. In the region I, r c ≤ r ≤ r I , K ( r, σ ) is positive and r I is sufficiently far fromthe turning point, so that quasi-classical approximation and hence the solution (A2) are valid. In the region III, r III ≤ r ≤ R cc , K ( r, σ ) is also positive and the point r III is again chosen far enough away from the turning point, sothat the solution in this region reads ξ = C K ( r, σ ) / exp " R R cc r p K (˜ r, σ ) d ˜ r Ω + C K ( r, σ ) / exp " − R R cc r p K (˜ r, σ ) d ˜ r Ω . (A14)In region II, r I ≤ r ≤ r III , we assume that the expansion (A8) is valid, and, moreover, | x | ≫ a at r = r I and r = r III .Boundary condition at r = r c prescribes the constants A and A to be of the same order of magnitude (see Eq.A5). As a result, the eigenfunctions in the region I at r = r I read (the decreasing exponents are neglected) ξ = A K ( r I , σ ) / exp " R r I r c p K (˜ r, σ ) d ˜ r Ω , (A15) T = dξdr = A K ( r I , σ ) / Ω exp " R r I r c p K (˜ r, σ ) d ˜ r Ω . (A16)8Similarly, the boundary condition at the core-crust interface, ξ ( R cc ) = 0 (see above), requires C = − C . Thus, theeigenfunctions in the region III at r = r III equal (we again omit the small exponent ∝ C ) ξ = C K ( r III , σ ) / exp " R R cc r III p K (˜ r, σ ) d ˜ r Ω , (A17) T = dξdr = C K ( r III , σ ) / Ω exp " R R cc r III p K (˜ r, σ ) d ˜ r Ω . (A18)Further, using the asymptotic behavior (A10)–(A13) of the eigenfunctions in the region II and continuity of eigen-functions ξ and T at r = r I and r = r III we find that B ≪ B and B ≪ B .At x = 0 eigenfunctions also have to be continuous, which means B U ( a, x ) = B U ( a, − x ) , (A19) B dU ( a, x ) dx = − B dU ( a, − x ) dx (A20)at x →
0. The above system has two solutions. First, U ( a,
0) = 0, B = − B . Second, dU ( a, x ) /dx | x =0 = 0, B = B .Since (e.g., [68]) U ( a,
0) = √ π a + Γ (cid:0) + a (cid:1) , (A21) dU ( a, x ) dx (cid:12)(cid:12)(cid:12)(cid:12) x =0 = − √ π a − Γ (cid:0) + a (cid:1) , (A22)the first and the second solutions correspond to, respectively, + a = − n and + a = − n , where n is a naturalnumber. Merging these constraints, we find the “quantization rule” for the parameter a : a = − − n . Recalling thedefinition of a , we get K = 2Ω p β (cid:18)
12 + n (cid:19) . (A23)This indirect condition on the frequency σ can also be rewritten in the form of the Bohr-Sommerfeld quantizationrule as Z r r p − K ( r, σ ) dr = π Ω (cid:18)
12 + n (cid:19) , (A24)where r and r are the turning points, K ( r , σ ) = K ( r , σ ) = 0. It is interesting to note that the condition(A23) could be immediately obtained from the analogy to the quantum harmonic oscillator problem [67]. This is notsurprising, since Eq. (A9) describes harmonic oscillations, while the boundary conditions imposed in the region IIrequire vanishing of the “wave-function” ξ at x → ±∞ . Therefore, this problem is indeed completely equivalent tothe classical problem of quantum mechanics, and leads to the same spectrum (A23). As in quantum mechanics, n inEq. (A23) determines the number of zeros of the function ξ and its parity.Now the asymptotic solution in the region II at | x | ≫ | a | can be represented as ξ = B U ( a, x ) , x ≫ | a | ,ξ = ( − n B U ( a, − x ) , − x ≫ | a | . (A25)This asymptotic solution should be matched with the solution (A15)–(A16) in the region I at r = r I (where − x ≫ | a | )and with the solution (A17)–(A18) in the region III at r = r III (where x ≫ | a | ). In this way one can relate thecoefficients A , B , and C , which completes our solution. We leave this exercise for the reader.Above we demonstrated how to build a solution for superfluid r-modes in the limit Ω →
0. Our results imply thatone needs to have a region with K ( r, σ ) < npeµ core. This region has to be small, r − r ∼ √ Ω, in orderfor the solution to possess a finite number of nodes, n . In other words, in the limit Ω → σ for SF r-modes is setby the condition min[ K ( r, σ )] = 0 , (A26)9i.e., the minimal value of K ( r, σ ) throughout the SF npeµ core must vanish. Obviously, σ for various harmonics (withfinite number of nodes) must coincide in this limit. The consideration above demonstrates that r-mode eigenfunctionscannot be expanded in Taylor series in Ω at Ω → non-analytic in this limit.We assumed above that minimum value of K ( r, σ ) corresponds to a minimum of the function K ( r, σ ), that is K ′ ( r , σ ) = 0. However, K ( r, σ ) can reach its minimum value at the stellar center or at the core-crust interface,where K ′ ( r, σ ), generally, does not vanish. Considering this situation in a way similar to that described above, onecan easily check that our main conclusions, in particular, Eq. (A26), remain valid also in this case.Let us now turn to the more realistic case of a three-layer star, containing npe layer in the outer core, between theinner npeµ core and crust. We limit ourselves to considering two possibilities. First, neutrons are SF all the way fromthe stellar center to some point r = r s inside the npeµ core. Then at r > r s neutrons are normal and eigenfunctionsobey the standard ordering at r s ≤ r ≤ R , ξ br,m +1 ,m ∼ Ω T b mm . Hence all the reasoning presented above remainsunchanged in this case, with the only difference, that now one should vanish ξ at r = r s : ξ ( r s ) = 0 + O (Ω).The second possibility is that neutron SF extends from the stellar center up to some point r = r s outside the npeµ core, r s > R µ . This situation differs slightly from that analyzed above. Now, while in the crust and in the normalpart of the npe core ξ br,m +1 ,m ∼ Ω T b mm , in the SF npe matter ξ br,m +1 ,m ∼ z r,m +1 ,m ∼ T b mm ∼ T z mm [19]. Thus,in this situation (and in the limit Ω →
0) not ξ , but the function T should vanish in the SF npeµ core interface, at r = R µ , to guarantee the continuity of eigenfunctions at r = R µ . Nevertheless, even in this case the conclusions ofthis Appendix [in particular, Eq. (A26)] remain unaffected. [1] N. Andersson, Astrophys. J. , 708 (1998), arXiv:gr-qc/9706075.[2] J. L. Friedman and S. M. Morsink, Astrophys. J. , 714 (1998), arXiv:gr-qc/9706073.[3] N. Andersson and G. L. Comer, Mon. Not. R. Astron. Soc. , 1129 (2001), astro-ph/0101193.[4] Y. Levin, Astrophys. J. , 328 (1999), astro-ph/9810471.[5] B. Haskell, International Journal of Modern Physics E , 1541007 (2015), 1509.04370.[6] K. Glampedakis and L. Gualtieri, Gravitational Waves from Single Neutron Stars: An Advanced Detector Era Survey(2018), vol. 457 of Astrophysics and Space Science Library, p. 673.[7] R. Bondarescu, S. A. Teukolsky, and I. Wasserman, Phys. Rev. D , 064019 (2007), 0704.0799.[8] M. G. Alford, S. Mahmoodifar, and K. Schwenzer, Phys. Rev. D , 044051 (2012), 1103.3521.[9] R. Bondarescu and I. Wasserman, Astrophys. J. , 9 (2013), 1305.2335.[10] B. Haskell, K. Glampedakis, and N. Andersson, Mon. Not. R. Astron. Soc. , 1662 (2014), 1307.0985.[11] M. E. Gusakov, A. I. Chugunov, and E. M. Kantor, Phys. Rev. D , 063001 (2014), 1305.3825.[12] M. E. Gusakov, A. I. Chugunov, and E. M. Kantor, Physical Review Letters , 151101 (2014), 1310.8103.[13] M. A. Alpar, S. A. Langer, and J. A. Sauls, Astrophys. J. , 533 (1984).[14] L. Lindblom and G. Mendell, Phys. Rev. D , 104003 (2000), gr-qc/9909084.[15] U. Lee and S. Yoshida, Astrophys. J. , 403 (2003), astro-ph/0211580.[16] E. M. Kantor, M. E. Gusakov, and A. I. Chugunov, Mon. Not. R. Astron. Soc. , 739 (2016), 1512.02428.[17] A. I. Chugunov, M. E. Gusakov, and E. M. Kantor, Mon. Not. R. Astron. Soc. , 291 (2017), 1610.06380.[18] E. M. Kantor and M. E. Gusakov, Mon. Not. R. Astron. Soc. , 3928 (2017), 1705.06027.[19] V. A. Dommes, E. M. Kantor, and M. E. Gusakov, Mon. Not. R. Astron. Soc. , 2573 (2019), 1810.08005.[20] E. M. Kantor, M. E. Gusakov, and V. A. Dommes, Phys. Rev. Lett. (2020).[21] T. G. Cowling, Mon. Not. R. Astron. Soc. , 367 (1941).[22] M. E. Gusakov and N. Andersson, Mon. Not. R. Astron. Soc. , 1776 (2006).[23] P. Haensel, A. Y. Potekhin, and D. G. Yakovlev, eds., Neutron Stars 1 : Equation of State and Structure, vol. 326 ofAstrophysics and Space Science Library (2007).[24] M. E. Gusakov, E. M. Kantor, and P. Haensel, Phys. Rev. C , 055806 (2009).[25] M. E. Gusakov, E. M. Kantor, and P. Haensel, Phys. Rev. C , 015803 (2009).[26] M. E. Gusakov, P. Haensel, and E. M. Kantor, Mon. Not. R. Astron. Soc. , 318 (2014), 1401.2827.[27] A. F. Andreev and E. P. Bashkin, Soviet Journal of Experimental and Theoretical Physics , 164 (1976).[28] G. Mendell, Astrophys. J. , 515 (1991).[29] N. Andersson, T. Sidery, and G. L. Comer, Mon. Not. R. Astron. Soc. , 162 (2006), astro-ph/0510057.[30] H. Saio, Astrophys. J. , 717 (1982).[31] K. H. Lockitch and J. L. Friedman, Astrophys. J. , 764 (1999), gr-qc/9812019.[32] S. Yoshida and U. Lee, Astrophys. J. Suppl. Ser. , 353 (2000), astro-ph/0002300.[33] J. Provost, G. Berthomieu, and A. Rocca, Astron. Astrophys. , 126 (1981).[34] N. Andersson, K. Glampedakis, and B. Haskell, Phys. Rev. D , 103009 (2009), 0812.3023.[35] H. Heiselberg and M. Hjorth-Jensen, Astrophys. J. Lett. , L45 (1999), astro-ph/9904214.[36] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C , 1804 (1998).[37] M. E. Gusakov and P. Haensel, Nuclear Physics A , 333 (2005), astro-ph/0508104. [38] S. Goriely, N. Chamel, and J. M. Pearson, Phys. Rev. C , 024308 (2013).[39] J. M. Pearson, N. Chamel, A. Y. Potekhin, A. F. Fantina, C. Ducoin, A. K. Dutta, and S. Goriely, Mon. Not. R. Astron.Soc. , 2994 (2018), 1903.04981.[40] E. M. Kantor and M. E. Gusakov, in Electromagnetic Radiation from Pulsars and Magnetars (2020), Journal of Physics:Conference Series.[41] D. Ding, A. Rios, H. Dussan, W. H. Dickhoff, S. J. Witte, A. Carbone, and A. Polls, Phys. Rev. C , 025802 (2016),URL https://link.aps.org/doi/10.1103/PhysRevC.94.025802 .[42] A. Sedrakian and J. W. Clark, European Physical Journal A , 167 (2019), 1802.00017.[43] M. E. Gusakov, A. D. Kaminker, D. G. Yakovlev, and O. Y. Gnedin, Astron. Astrophys. , 1063 (2004), astro-ph/0404002.[44] M. E. Gusakov, A. D. Kaminker, D. G. Yakovlev, and O. Y. Gnedin, Mon. Not. R. Astron. Soc. , 555 (2005), astro-ph/0507560.[45] P. S. Shternin, D. G. Yakovlev, C. O. Heinke, W. C. G. Ho, and D. J. Patnaude, MNRAS , L108 (2011).[46] K. G. Elshamouty, C. O. Heinke, G. R. Sivakoff, W. C. G. Ho, P. S. Shternin, D. G. Yakovlev, D. J. Patnaude, andL. David, Astrophys. J. , 22 (2013), 1306.3387.[47] D. Page, J. M. Lattimer, M. Prakash, and A. W. Steiner, Astrophys. J. Suppl. Ser. , 623 (2004), astro-ph/0403657.[48] D. Page, M. Prakash, J. M. Lattimer, and A. W. Steiner, Phys. Rev. Lett. , 081101 (2011).[49] W. C. G. Ho, K. G. Elshamouty, C. O. Heinke, and A. Y. Potekhin, Phys. Rev. C , 015806 (2015), 1412.7759.[50] S. Beloin, S. Han, A. W. Steiner, and D. Page, Phys. Rev. C , 015804 (2018).[51] U. Lombardo and H.-J. Schulze, in Physics of Neutron Star Interiors, edited by D. Blaschke, N. K. Glendenning, andA. Sedrakian (2001), vol. 578 of Lecture Notes in Physics, Berlin Springer Verlag, p. 30.[52] D. G. Yakovlev, K. P. Levenfish, and Y. A. Shibanov, Sov. Phys.— Usp. , 737 (1999).[53] A. Gezerlis, C. J. Pethick, and A. Schwenk, ArXiv e-prints (2014), 1406.6109.[54] J. M. Dong, U. Lombardo, and W. Zuo, Physics of Atomic Nuclei , 1057 (2014).[55] N. Andersson and K. D. Kokkotas, International Journal of Modern Physics D , 381 (2001), gr-qc/0010102.[56] A. S. Parikh and R. Wijnands, Mon. Not. R. Astron. Soc. , 2742 (2017), 1707.05606.[57] F. ¨Ozel, D. Psaltis, R. Narayan, and A. Santos Villarreal, Astrophys. J. , 55 (2012), 1201.1006.[58] J. Antoniadis, T. M. Tauris, F. Ozel, E. Barr, D. J. Champion, and P. C. C. Freire, arXiv e-prints arXiv:1605.01665 (2016),1605.01665.[59] M. E. Gusakov, E. M. Kantor, A. I. Chugunov, and L. Gualtieri, Mon. Not. R. Astron. Soc. , 1518 (2013), 1211.2452.[60] L. Gualtieri, E. M. Kantor, M. E. Gusakov, and A. I. Chugunov, Phys. Rev. D , 024010 (2014), 1404.7512.[61] Y. Levin and G. Ushomirsky, Mon. Not. R. Astron. Soc. , 917 (2001), astro-ph/0006028.[62] K. Glampedakis and N. Andersson, Phys. Rev. D , 044040 (2006), astro-ph/0411750.[63] M. Nayyar and B. J. Owen, Phys. Rev. D , 084001 (2006), astro-ph/0512041.[64] M. G. Alford, S. Mahmoodifar, and K. Schwenzer, Phys. Rev. D , 024007 (2012), 1012.4883.[65] D. D. Ofengeim, M. E. Gusakov, P. Haensel, and M. Fortin, Phys. Rev. D , 103017 (2019), 1911.08407.[66] B. Haskell, N. Andersson, and A. Passamonti, Mon. Not. R. Astron. Soc.397