Nonanalytic Relativistic r -Modes of Slowly Rotating Nonbarotropic Neutron Stars

: We show that the r -modes of slowly rotating nonbarotropic neutron stars are described by nonanalytic functions of stellar angular velocity, which makes the perturbation techniques, used so far in the r -mode theoretical studies, inapplicable. In contrast to those studies and in accordance with numerical calculations beyond the slow rotation approximation, the obtained r -mode spectrum is discrete, which resolves the continuous spectrum problem, lasting since 1997. Our ﬁndings imply that the relativistic r -modes in slowly rotating neutron stars dramatically differ from their Newtonian cousins, which may have important implications for the detectability of r -mode signatures in observations, in particular for the r -mode excitation efﬁciency during the neutron star inspirals.

A standard approach in theoretical studies of r-modes in slowly rotating stars (i.e., stars with angular velocity, Ω, significantly lower than the Keplerian velocity, Ω K ) is to assume that the r-mode eigenfunctions and eigenfrequencies can be Taylor-expanded in small parameter Ω/Ω K and the r-mode equations can be obtained following the perturbation theory with respect to this parameter. While this traditional approach turns out to be very fruitful in the Newtonian theory (see [19,20] and the references therein), in relativistic theory, it leads to a set of contradictory results, known as "the problem of the continuous r-mode spectrum". In General Relativity (GR), the frame-dragging effect at position r leads to the entrainment of local inertial reference frames with angular velocity Ωω(r), as measured by the distant observer. As a result, in barotropic stars (where the pressure is a function of one thermodynamic variable), r-modes cease to exist [21], while in nonbarotropic stars (the pressure depends on several variables), the traditional approach predicts the r-mode frequencies, σ, to fill the continuous band [21][22][23][24][25][26][27]: where l and m are the r-mode quantum numbers and R is the stellar radius. The continuous spectrum is not affected by accounting for the gravitational radiation in the problem [28,29]; the discrete r-modes [21,30] with frequencies beyond the band (1) very likely do not exist in typical neutron stars [30], and the isolated r-modes [27] with hidden discrete frequencies within the band (1) have divergent velocity perturbations (see also a discussion of different spectrum regularization methods in [31][32][33]). Numerical calculations [34,35] beyond the slow rotation approximation, in turn, predict a discrete r-mode spectrum, as in the Newtonian theory, thus indicating the failure of the traditional approach in application to relativistic r-modes. In this letter, we develop a new original approach to the r-mode problem, which reveals the reason for that failure and uncovers quite unexpected properties of relativistic r-modes.

General Oscillation Equations
We considered a neutron star as a non-magnetized quasineutral mixture of normal (i.e., nonsuperfluid and nonsuperconducting) degenerate particle species (further labeled by Latin indices i, j, k, . . . ). Thermodynamic functions in this case-pressure p, energy density ε, enthalpy density w = p + ε, and chemical potentials µ k -depend only on the set of number densities n k of different particle species. For simplicity, we ignore perturbations of the gravitational field (Cowling approximation), which leads to a relative error in the r-mode oscillation frequency of only a few percent [36]. Then, gravity plays the role of the stationary background, described in the x µ = (ct, r, θ, ϕ)-coordinates (c is the speed of light) by Hartle's metric tensor g µν of a slowly rotating massive body [37]: Under these assumptions, small stellar oscillations are governed by the linearized Euler equation, linearized continuity equations for different particle species, thermodynamic relations in degenerate matter, and the equation of state (EOS), assumed in this study to be nonbarotropic (summation over all repeated indices is implied): Here, u µ is the four-velocity of the macroscopic material flows inside the star (u µ = u µ = (δ µ t + Ω/c δ µ ϕ )u t in equilibrium), and covariant derivative ∇ ρ is associated with the interval (2). Due to the axial symmetry of a rotating star, the periodic Eulerian perturbation δ f of any quantity f can be sought in the form δ f ∝ e iσt+imϕ .
Below, we express δu µ through the Lagrangian displacement ξ ρ , showing the variation of the world lines of fluid elements in the star due to an oscillation (see, e.g., [38]): The gauge freedom ξ ρ → ξ ρ + f u ρ in the definition (4), where f is some arbitrary function, allows us to impose the additional algebraic condition u µ ξ µ = 0.

Deriving r-Mode Equations
Based on the results of the Newtonian theory, one expects relativistic r-modes to be oscillations with frequency σ ∝ Ω and quasi-toroidal displacement (by definition, f (1) f (0) ): where we introduced the toroidal function T to be found below. Although for σ ∝ Ω, only even powers of Ω enter the equations (3) and it seems natural to assume that f (1) ∝ Ω 2 (for brevity, below, we omit the 1/Ω K factor in Ω/Ω K ), such an approach leads to the continuous spectrum problem. Here, to facilitate the analysis of the equations, we assumed that the frame-dragging effect, which is the cause of the continuous spectrum, is weak, ω(r) = ω(r), where = ω(0) 1. We also ignored the oblateness of the star (O(Ω 2 ) terms in Equation (2)), which leaves the mathematical structure of the problem unaffected. Then, we look for the solution of the oscillation Equation (3) in the form (5), assuming that the terms f (0) correspond to the limit → 0 and Ω 2 / → 0 (i.e., we account for the frame-dragging effect for arbitrarily slow rotation rates), while the terms f (1) , whose dependence on and Ω is not postulated, should be small due to small Ω/Ω K , a weak frame-dragging effect, or both.
To derive general r-mode equations, we followed three subsequent steps [33]. First, we excluded scalar perturbations δn k , δp, and δw from the system (3). Second, we used decompositions (5) in the thus-reduced system and in the → 0 and Ω 2 / → 0 limit found: where P m l (cos θ) is the associated Legendre polynomial. Third, we simplified the reduced system using (7) and discarded small terms, according to the following selection rule: if there is a term f in a chosen equation, then the terms Ω 2 f and f in this equation are ignored. An analysis of the thus-obtained equations reveals that the r-mode radial displacement has the form and that the problem reduces to the solution of the closed system of equations for T l (r) and ξ r l±1 (r), which for the l = m case (most CFS unstable mode) can be written as The explicit form [33] of the coefficients C 1,2,3 (r) and G 1,2 (r) is not important for further discussion. Near the stellar center, the solution should smoothly match its asymptote T l (r) ∝ r l [33]. At the stellar surface, the total pressure should vanish, which in our case is equivalent to [33] where the prime denotes d/dr. Those values of σ (1) , for which these boundary conditions can be simultaneously met, form the relativistic r-mode spectrum. In the Newtonian limit, Equation (9) and boundary conditions correctly reduce to those of traditional Newtonian theory.

R-Modes in the Ω → 0 Limit
It is easy to see that, if ξ r l+1 ∝ Ω 2 , as in the traditional ordering, then the first equation of the system (9) leads to a continuous spectrum in the Ω → 0 limit. At the same time, by numerically solving the system, we obtain regular r-mode eigenfunctions and discrete eigenfrequencies with no signatures of the continuous spectrum (see below). To find the ordering, corresponding to this solution, we note that for extremely slow rotation rates, the leading contribution to σ (1) is defined by the weak frame-dragging effect (rotational corrections become negligible). Since σ (1) andω(r) enter the system (9) only in combination σ (1) + 2 ω(r)/(l + 1), this implies that in this limit, σ (1) ∼ is a quantity of order . Then, we assume that, in the Ω → 0 limit, the leading contributions to the r-mode frequency correction σ (1) and eigenfunctions behave as where σ (10) ∼ 1. The latter condition in (11) means that the operator d/dr, when acting on the eigenfunctions, should be considered as a "quantity" of order d 1 Ω d 2 . The inequalities d 1 = 0 or d 2 = 0 indicate that the functions under consideration are nonanalytic in or Ω (see below for a further explanation). Under these assumptions, some terms in the system (9) can be small compared to others. Since this system should allow finding the eigenfrequency correction, the terms with σ (10) andω(r) should be of the same order as at least one of the other terms in the first equation, while the remaining terms are allowed to be comparatively small. In the second equation, at least two terms should be of the same order and not smaller than the remaining term; otherwise, the solution is trivial. The analysis shows that these conditions are met only for k 1 = d 1 = 1/2 and k 2 = −d 2 = 1; therefore, Using this result, we discard small terms in the system (9) and eventually reduce it to In the vicinity of the stellar center, 0 ≤ r ≤ r c , some of the discarded terms are not small, and, therefore the solution to Equation (13) does not satisfy (9) (the value of r c can be estimated from (9) and tends to zero for small κ). Further, from (9) and the asymptote T l ∝ r l of the toroidal function, it follows that, near the center, ξ r l+1 ∼ κ 2 T l . The leading contribution to the radial displacement ξ r l+1 ∼ κT l cannot satisfy this condition and, therefore, should vanish when approaching r c . Similarly, it should also vanish at the stellar surface. As a result, the boundary conditions for the leading contribution to the radial displacement are ξ r l+1 (r c ) = ξ r l+1 (R) = 0. In order to satisfy the "correct" boundary conditions, discussed after (9), one should account for the higher-order rotational corrections to the radial displacement.
For a simple toy model with q σ (r) = q = const, we obtain T l ∝ e ± √ qr/κ for q > 0 and T l ∝ sin( |q|r/κ + φ) for q < 0. These nonanalytic functions of κ produce, under the action of the operator d/dr, the 1/κ factor, which reveals the idea behind the notation d/dr ∼ 1/κ in (12). In reality, however, the function q σ (r) is not constant and can have turning points, where it changes sign. Since the explicit form of the coefficients G 2 (r) and C 1 (r) implies that G 2 (r)/C 1 (r) > 0 [33] andω(r) is a positive monotonically decreasing function [37], we can have only one turning point r t , related (if exists) to the eigenfrequency as (l + 1)σ (10) + 2ω(r t ) = 0.
Based on the toy model, we expect the toroidal function will behave exponentially at r c < r < r t and oscillate at r t < r < R. The r-mode eigenfunctions can then be classified by the number of nodes (zeros) n of the toroidal function. For any n, the distance between the nodes tends to zero as κ → 0; therefore, the turning point r t,n → R and eigenfrequencies σ (10) n → −2ω(R)/(l + 1). To find the approximate solution of (13) in the κ → 0 limit, we used the WKB method. Far from the turning point, where q σ (r) noticeably differs from zero (region I with r c < r < r t and region III with r t < r < R), we looked for the solution in the form T l (r) = e h(κ,r)/κ , where h(κ, r) is some analytic function of κ. In the vicinity of the turning point (region II), we used the Taylor expansion q σ (r) ≈ α 2 (r t − r), introduced a new variable z = (r t − r)(α/κ) 2/3 , and reduced the second-order Equation (13) to the Airy equation. In the regions, where q σ (r) noticeably differs from zero, but its Taylor expansion is still accurate, the solutions in regions I and III should, with high accuracy, coincide with that found in the region II. This condition and the equality ξ r l+1 (r c ) = 0 allow us to find the r-mode eigenfunctions: The remaining boundary condition at the surface ξ r l+1 (R) = 0 then leads to the quantization rule: which, being combined with (15), defines the discrete r-mode spectrum in the κ → 0 limit. For those κ, for which the turning point is so close to the stellar surface that the Taylor expansion of q σ (r) can be applied up to r = R, one can obtain the following explicit expressions for the turning point and for the r-mode spectrum:

Gravitational Radiation Timescale
The instability with respect to the emission of gravitational waves causes the r-mode energy E to slowly increase with average rateĖ GW , leading to the mode growth on a typical timescale τ GW = 2E/Ė GW . A reasonable and commonly used way to estimate the efficiency of this process is to adopt the general post-Newtonian expressions for E andĖ GW [10,39]. In our case, τ GW is mainly defined by the current multipole moment and equals where γ is some (weakly depending on Ω) combination of the fundamental constants G and c, stellar mass M and radius R, the mode order l, and the total frequency σ.
Since Ω nonanalytically enters the toroidal function, the dependence τ GW (Ω) becomes highly non-trivial and significantly differs from the traditional τ GW ∝ 1/Ω 2l+2 . In particular, for those rotation rates, for which the Taylor expansion q σ (r) ≈ α 2 (r t − r) is valid up to r = R, one can replace the toroidal function T with T II (see Equation (16)) and expand in Taylor series other, slowly changing subintegral functions near the stellar surface. Then, one can show that, in this limit, τ GW ∝ 1/Ω 2l+2+4/3 .

Numerical Results
To verify the predictions of the theory in the Ω → 0 limit, we calculated the l = m = 2 r-mode spectrum and eigenfunctions from (9) for different Ω. We employed the BSK24 equation of state [40] and considered a neutron star model with mass M = 1.4M (M is the solar mass).
According to our calculations, for Ω/Ω K 0.05, relativistic toroidal eigenfunctions weakly depend on Ω and resemble their Ω-independent Newtonian counterparts, while for Ω/Ω K < 0.05, as Ω decreases, they start to be localized near the stellar surface. Figure 1, in which we plot the relativistic (GR, blue line) and Newtonian (Newt, red dashes) toroidal functions with four nodes, shows the example of such behavior at Ω/Ω K = 0.005. In accordance with (16), the eigenfunction is exponentially suppressed at r < r t and oscillates at r > r t . In our calculations, performed for Ω/Ω K ∈ [0.1; 0.001], we found no signatures of the continuous part in the r-mode oscillation spectrum. The spectrum is discrete and, as one can see from Figure 2, coincides for extremely slow rotation rates with that predicted by the theory in the Ω → 0 limit. Here, the eigenfrequencies, obtained numerically via solving Equation (9), are shown by filled circles, while curves show the spectrum, obtained from the quantization rule (17). Different colors correspond to different numbers of nodes n of the toroidal function.

Conclusions
Summarizing, we applied a new original approach to study relativistic r-modes of slowly rotating nonbarotropic neutron stars. Unlike the traditional one, it does not require any Taylor expansions in Ω/Ω K of the r-mode eigenfunctions and oscillation spectrum, and the only essential assumptions made are the stellar rotation being slow and the framedragging effect being weak. Within this approach, we explicitly showed that the relativistic generalization of the Newtonian r-modes in nonbarotropic stars is described by nonanalytic functions of the stellar angular velocity. Because of this nonanalyticity, these functions are undetermined at Ω = 0 (see Equation (16)), and relativistic r-modes cannot be considered as small rotational corrections to the zero-frequency modes of the non-rotating star. They cannot be Taylor-expanded in small parameter Ω/Ω K and, therefore, cannot be studied within the traditional approach used so far.
Nonanalyticity affects the r-mode ordering: instead of the "classical" ξ r l+1 ∼ Ω 2 T l and σ (1) ∝ Ω 2 , we have ξ r l+1 ∼ √ ΩT l , σ (1) ∼ , and moreover, the operator d/dr ∼ √ /Ω, when acting on the relativistic eigenfunctions. These results may affect the interpretation of the observed periodic oscillations in electromagnetic radiation of rotating NSs [41,42]. Further, analysis of the general oscillation equations (Equation (3)) shows [33] that, while the pressure perturbation δp ∼ Ω 2 T l retains its usual order, other scalar perturbations δε ∼ δn k ∼ δµ k ∼ √ ΩT l instead of being of order Ω 2 , as in the Newtonian case. Such nonanalyticity and peculiar ordering should significantly influence the typical mode dissipation timescales due to shear viscosity τ η , bulk viscosity τ ζ , and diffusion τ diff [43,44]. We expect the effect to be especially pronounced for τ ζ and τ diff , defined by non-toroidal components of the motion. Moreover, we expect that such ordering of scalar functions should modify the mode coupling to the tidal potential and affect the r-mode excitation during NS inspiral, which may imprint in the gravitational wave signal (e.g., [45]). Finally, according to our post-Newtonian estimate, because of the nonanalyticity, the mode growth timescale τ GW due to the CFS mechanism has non-standard dependence τ GW (Ω), reducing for extremely slow rotation rates to τ GW ∝ 1/Ω 2l+2+4/3 instead of the traditional τ GW ∝ 1/Ω 2l+2 .
The performed numerical calculations show that the nonanalytic r-mode spectrum is discrete and, for extremely slow rotation rates, is accurately described by the quantization rule (17) (which, notably, has no counterparts in the Newtonian theory). Relativistic rmode eigenfunctions in this limit drastically differ from the Newtonian ones: they are exponentially suppressed everywhere in the star except for the narrow region in the vicinity of the stellar surface, where they oscillate. Such behavior makes the employed Cowling approximation even more justified: the r-mode dynamics can be significantly influenced only by the gravitational field perturbations of the internal stellar layers, where the modes are suppressed. Such behavior also justifies the assumption of the weak frame-dragging effect: near the surface, where the r-modes are localized, we have ω(R) ≈ 0.15, and it is, indeed, weak. We note that, according to our preliminary results, the r-mode nonanalyticity takes place beyond the weak framedrag approximation. However, the r-mode equations are more complicated in this case and deserve separate consideration. Finally, accounting for the oblateness of the rotating star produces additional terms depending on Ω 2 in the general system (9). Generally, these terms should be accounted for, but at Ω → 0, they become small and do not affect the results, obtained in this limit.
Although we have managed to obtain discrete relativistic r-modes, our study, strictly speaking, did not completely exclude the existence of the r-modes with the continuous oscillation spectrum discussed in the literature [21][22][23][24][25][26][27][28][29]. We believe, however, that, since the continuous spectrum does not manifest itself in our numerical calculations, as well as in the numerical calculations of [34,35], such r-modes possibly do not exist as a result of some internal inconsistency of the theory, based on the traditional approach.
The mathematical structure of the equations, describing relativistic r-modes, resembles that of the Newtonian r-mode equations in differentially rotating stars. The ideas, developed in this letter, can be applied to study these equations and may lead to the discovery of a new class of nonanalytic solutions, similar to the nonanalytic relativistic r-modes. In particular, this may affect the modeling of the recently observed solar Rossby waves [46]. Keeping in mind that solar rotation is strongly differential, while Ω/Ω K is extremely small, one could expect the peculiar behavior of solar Rossby waves analogous to those discussed in this letter.